Designing a New Shuttle Service to Meet Large-Scale Instantaneous Peak Demands for Passenger Transportation in a Metropolitan Context: A Green, Low-Cost Mass Transport Option

Currently, the green, sustainable development of metropolises is hindered by problems caused by Large-scale Instantaneous Peak-demands for Passenger-transportation (LIPP), such as traffic congestion and air pollution. To mitigate these problems, we propose a new type of demand-responsive service as an alternative to inefficient “door-to-door” service. The proposed service is based on service units designed to aggregate passengers for shuttle service. To guarantee service quality and efficiency, a maximum passenger walking time constraint, a request rejection mechanism and a scheme for ensuring solution feasibility are considered. Through numerical experiments, we prove the following: (i) the proposed transport option exhibits better performance (by 40.37% for passengers and by 35.79% for operators) than the door-to-door transport option for solving real cases. (ii) By testing different datasets, we prove that the proposed service is more suitable for the request distributions that are spatiotemporally concentrated. (iii) Regarding the individual components of the proposed clustering-first, routing-second solution framework, the proposed soft clustering algorithm exhibits better performance than the classical hard clustering method (by 8%), and the proposed routing algorithm is 1.5 times more efficient than the commercial solution software GAMS.


Background
As very common phenomena in the development of metropolises, Large-scale Instantaneous Peak-demands for Passenger-transportation (LIPP) have led to a number of problems related to mobility management and environmental concerns, such as traffic congestion, air pollution, and traffic accidents [1][2][3]. These phenomena are of great concern to the managers of metropolises because of the characteristics of the corresponding demands. These characteristics are markedly different from those of the transportation demands in general cities in the following ways.
(i) The peak demands in metropolises are unpredictable and instantaneous. Multiple sources (e.g., commuters and tourism traffic) and many factors (e.g., weather [4], festivals and special events) can cause demands to fluctuate irregularly and without warning.
(ii) The high population densities and clear boundaries of urban functional areas result in peak demands with certain spatiotemporal similarities in terms of origins and destinations (e.g., demands of commuters who live in the same residential areas and work in the same business districts). Regarding methodologies, combinations of request clustering and routing methods based on geographical similarities have been the focus of existing research. Cullen et al. [28] proposed a clustering-first, routing-second heuristic model. Bodin and Sexton [29] also constructed a clustering-first, routing-second heuristic model for the pickup and delivery problem with time windows (PDPTW). In this model, passengers are first partitioned into clusters. Then, the clusters are assigned to routes using the algorithm developed by Sexton and Bodin [30,31]. Finally, the algorithm swaps passengers between routes for route re-optimization. Dumas et al. [6] introduced the concept of mini-clusters, in which passengers with high spatiotemporal proximity are clustered together. In this algorithm, a heuristic method is used to produce a set of mini-clusters. Then, a column generation (CG) algorithm is used to optimally combine the mini-clusters into vehicle routes. Finally, the algorithm performs route re-optimization to obtain the optimal schedule for each vehicle. Desrosiers et al. [32] proposed another approach to mini-cluster construction by means of a parallel insertion method based on the spatiotemporal proximity of the requests. Ioachim et al. [33] applied an optimization-based technique instead of a heuristic method for the construction of mini-clusters. In addition to these studies, Pankratz [34], Bard and Jarrah [35], Qu and Bard [36], and Masson et al. [37] also applied clustering algorithms for the PDPTW. The results of these studies show that high-quality vehicle plans can be obtained only through the integrated consideration of clustering and routing.
In the field of operations research, a routing problem for demand-responsive service is typically formulated, as either a vehicle routing problem with time windows (VRPTW) [38,39] or a PDPTW [40,41]. In most VRPTWs or PDPTWs, the number of passengers served, the operating cost, and the service quality are the three major objectives to be optimized, either separately or simultaneously [42]. Typical vehicle-to-passenger matching and routing constraints include time window constraints, capacity constraints, coupling constraints, and maximum ride time constraints, which typically give rise to an NP-hard problem [43,44]. Many researchers have focused on studying corresponding solution techniques. CG algorithms [45,46], branch-and-price algorithms [15], branch-and-cut algorithms [47], branch-and-cut-and-price algorithms [48], tabu search algorithms [49], and dynamic-programming-based metaheuristics [14] can be used to solve problems of this type. The proposed LIPDRS problem can be regarded as an offline multivehicle routing problem because all the requests are assumed to be known beforehand.
Although significant improvements have been achieved in solving demand-responsive service design problems, LIPDRS design still poses several challenges, as discussed in Section 1.1. To overcome these challenges, we model the LIPDRS decision process using a mixed-integer programming model with dual objectives to optimize the service from the perspectives of both passengers and operators. This model includes a carefully designed passenger walking time constraint and request rejection mechanism to ensure that the resulting service plan can guarantee not only good quality of service, but also improved operational efficiency. Additionally, this paper proposes a clustering-first, routing-second framework for solving the proposed model. A modified version of the fuzzy c-means clustering algorithm and a CG algorithm are utilized in this framework to determine the service units and perform vehicle assignment and routing, respectively. In addition, to overcome the challenges caused by relaxing the model variables, a branching scheme is proposed.
The rest of the paper is organized as follows. Section 2, presents details of the methodology for solving the LIPDRS design model, including the (first) clustering phase, i.e., the determination of service units (DSU), and the (second) routing phase, i.e., fleet assignment and routing (FAR). Section 3 presents a numerical experiment conducted to analyze the performances of the proposed services on different request datasets. In addition, we present comparisons between the proposed LIPDRS option and the existing door-to-door transport option based on various artificial datasets to evaluate the application scope of our proposed method.

Methodology
In this section, the methodology of the present study is presented. We first introduce the whole architecture for LIPDRS design, including the overall service decision process and the construction of the LIPDRS design model (LIPPD). Second, we introduce the corresponding clustering-first, routing-second algorithm. Finally, we introduce the criteria for evaluating the results.

Overall Service Decision Process
To enable an effective and feasible shuttle service, a seven-step process in which the requesters can be actively involved is designed, as shown in Figure 1.

1.
Request Collection. Potential passengers submit their travel requests (each consisting of a preferred location tuple with boarding and alighting locations and a time window tuple with the expected boarding and alighting times) through mobile phones or other devices. The data collection and planning period is designated as one day; that is, in one day, the demands for the next day are collected and planned. Any request that is submitted before the deadline (either a specified moment, such as 12:00 p.m., or some specified number of hours ahead of the service time) will be pushed into the request pool to wait for stage 2. Other requests will be marked as overdue requests and will wait for stage 5.

2.
Determination of Service Units (DSU). Requests in the request pool are aggregated into a set of service units (which specify the locations at which and time windows in which the requesters will be served) using the approach described in Section 2.3. The set of service units will be sent to the service pool to wait for stage 3. Isolated requests (i.e., requests that do not belong to any service unit) from the request pool are also detected in this stage. These isolated requests will wait for stage 4.

3.
Fleet Assignment and Routing (FAR). Vehicles are assigned to specific service units, and the fleet is routed, subject to capacity constraints, spatiotemporal constraints, procedure constraints, etc. Through this process, the operators obtain the service routes and scheduled service plan (timetable) for the LIPDRS. 4.
Economic Evaluation. The operators use evaluation metrics to compare the generated plan with other transport options and evaluate whether the plan will be economically acceptable to passengers. These evaluation metrics are discussed in Section 2.5.

5.
Request Rejection. To ensure the quality and efficiency of demand-responsive transport, reasonable rejection of certain transport requests is an integral part of designing a transport plan. Three types of requests need to be rejected according to different rejection rules: overdue requests, isolated requests, and uneconomical requests. These three types of rejection occur in different parts of the planning process, as shown in Figure 1. The specific rejection mechanism will be discussed in Section 2. 6.
Information Feedback. This stage consists of two parts: feedback to requesters who have been assigned to service units, with scheduled boarding and alighting information, and feedback to rejected requesters. 7.
Service Implementation. After receiving their feedback information, passengers will walk or bike to their specified boarding locations within their specified time windows for service.

LIPDRS Design Model (LIPPD)
The main idea on which the formulation of the LIPPD is based is that it is necessary to first determine service units based on the received requests and then assign and route the fleet in accordance with those service units. We formally formulate the LIPPD as follows: Objectives:

LIPDRS Design Model (LIPPD)
The main idea on which the formulation of the LIPPD is based is that it is necessary to first determine service units based on the received requests and then assign and route the fleet in accordance with those service units. We formally formulate the LIPPD as follows: Objectives: x i∈V,j∈V The related objective functions, constraints, variables, sets and parameters are defined as follows. • Objective Functions Equation (1a) represents the cost of operation, including the fixed and variable costs of using vehicles and the penalties for the cancelation (rejection) of service units. A smaller objective function value is better. Equation (1b) guides the DSU process. Unlike in the case of taxi and door-to-door services, LIPDRS requesters will be served at specified boarding/alighting locations (i.e., the locations of the service units) instead of the requesters' home/work locations [50]. The number of service units impacts the efficiency, effectiveness and equity of passenger service and vehicle operation. If the service units are too densely distributed (meaning that the clusters are small in size), the speed of service will not be effectively accelerated; if the service units are too sparsely distributed, more requesters will have to walk farther. Therefore, we select the ratio between compactness and divergence as one of the objective functions of the LIPPD, where (a) the divergence drives the distances between service units to be as large as possible and (b) the compactness drives the walking distances (times) between requesters and service units to be as small as possible. In this case, a larger objective function value is better. To extend the problem considered in this paper to other fields, different objective functions can be chosen from existing studies [7,51,52].

• Constraints
Equation (2) is a spatiotemporal constraint, namely, a maximum walking time constraint, where mwt denotes the maximum walking time; the norm · At W will be explained in Section 2.2.2. For example, in Figure 2, there are 17 requests and 6 candidate service units. As a result of the DSU process, 2 of the candidate service units are selected, and the 17 requests are grouped into service unit 2 and service unit 5. Within each service unit, the time windows for boarding and alighting satisfy the intracluster operation time constraint. Service unit 2 consists of 14 close requests, and the walking time between each requester and the service location is less than the predefined maximum walking time mwt. Equation (3) ensures that the service units are selected from among the candidate service units because the vehicles should not park at arbitrary locations. Equation (4) ensures that each service unit is assigned to exactly one vehicle or is canceled. Equation (5) enforces that a vehicle must visit both the boarding and alighting points of every service unit assigned to it. Equations (6)-(8) ensure that the vehicles either perform assigned tasks or remain idle. Equations (9)-(12) together form the precedence constraints. Equations (13)-(15) together form the capacity constraints. Equations (16) and (17) represent the maximum distance constraint and the maximum driving time constraint, respectively.

• Decision Variables
The variables of this model are as follows. x k ij (i, j) ∈ (V × V)· (De + , j) j ∈ V · ( j, De − ) j ∈ V , k ∈ M is a binary variable that is equal to 1 if vehicle k is scheduled to travel from location i to location j and is equal to 0 otherwise. z k i (i ∈ N, k ∈ M) is a binary variable that is equal to 1 if service unit i is assigned to vehicle k and is otherwise equal to 0. c, a positive integer variable, represents the number of clusters (i.e., the number of service units). The u ji ( c−1 j=1 u ji = 1, 0 ≤ u ji ≤ 1) are positive variables. U n×c is the matrix of membership values, and u ji ∈ U n×c is the membership value of the j th request with respect to service unit i. D i (i ∈ V ∪ De) is an integer variable that specifies the departure time at vertex i. y i (i ∈ V ∪ De) is an integer variable that specifies the load of the vehicle when arriving at vertex i. We define q De + = 0 for all k ∈ M.
• Sets In this model, N denotes the set of service units. For each service unit i ∈ N, a request group of size q i ∈ N needs to be transported from boarding point i + to alighting point i − . It is known that q i = q i + = −q i − , where positive quantities correspond to pickups and negative quantities correspond to dropoffs. We define N + i + i ∈ N as the set of boarding points and N − i − i ∈ N as the set of alighting points. Let V N + ∪ N − . Furthermore, M represents the set of vehicles. Each vehicle k ∈ M has a capacity Q ∈ N, a starting location De + , and an ending location De − . We set De De + ∪ De − .

• Parameters
For all i, j ∈ V ∪ De, c ij denotes the travel cost. F represents the fixed cost of a vehicle, and P represents the penalty for cancelling serving a service unit. Then, m is a fuzzifier with a value greater than 1, where a higher value indicates a higher degree of ambiguity [51]. n is the number of requests. r represents the average of the spatiotemporal features of all requests. r j represents the spatiotemporal features of the jth request, and v i represents the spatiotemporal features of service unit i. These spatiotemporal features consist of location and time window information. We represent a request as an eight-dimensional vector r j = lo j+ , tw j+ , lo j− , tw j− = lng j+ , lat j+ , e j+ , l j+ , lng j− , lat j− , e j− , l j− , where '+' indicates features corresponding to the origin and '-' indicates features corresponding to the destination. Similarly, each service unit is associated with an eight-dimensional feature vector of the same form, v i = (lo i+ , tw i+ , lo i− , tw i− ) = (lng i+ , lat i+ , e i+ , l i+ , lng i− , lat i− , e i− , l i− ), where '+' indicates boarding information and '-' indicates alighting information. Each candidate service unit is similarly represented in the form v i = (lo i+ , tw i+ , lo i− , tw i− ) = (lng i+ , lat i+ , e i+ , l i+ , lng i− , lat i− , e i− , l i− ). The mwt denotes the maximum walking time. Then, we denote the set of requests by R = r j 0 ≤ j < n , the set of service units by S = {v i |0 ≤ i < c} (where c is the number of service units) and the set of candidate service units by Cad = {v i |0 ≤ i < c }, where n, c and c are the numbers of elements of S, R and Cad, respectively. [e i , l i ] denotes the time window, and op i denotes the operation time. For all i, j ∈ V ∪ De, d ij denotes the travel distance, and t ij denotes the travel time. For each service unit i ∈ N, the number of requests q i ∈ N is given. Q is the capacity of vehicle. Additionally, each vehicle has an operation distance range [0, d max ] and an operation time range [0, t max ] (representing the driver's working time).
Based on the description given in Figure 2, we further illustrate the model by means of a simple example of the LIPDRS design process. Objective function (1.2) and constraints (2) and (3) are used to select service units from among the candidate service units and assign requests to them in accordance with the spatiotemporal features of the requests. The bottom of Figure 2 shows the process of selecting service units. In this example, the requesters submit requests that contain spatiotemporal origin and destination information, represented by oblique-textured and cross-textured circles, respectively. The candidate service units are represented by untextured squares. Ultimately, two service units are selected, represented by textured squares; they are identified as service unit 2 and service unit 5. These two service units serve two branches of requests with different destinations. Each service unit has two components: the boarding component (oblique-textured squares) and the alighting component (cross-textured squares). Both components contain the corresponding locations and time windows In this example, five service units are selected, and the operator assigns them to vehicles in the fleet and designs routes for the fleet. Objective function (1.1) and constraints (3)- (20) are used for assignment and routing. Once the service routes have been determined, the fleet departs from the depot(s), visits each service unit within the corresponding service time window, and then returns to the depot. Generally, the LIPDRS solution comprises several service routes, such as the two example routes (0-1-1 -2-2 -0 and 0-4-5-4 -5 -3-3 -0) shown in Figure 2.  Based on the description given in Figure 2, we further illustrate the model by means of a simple example of the LIPDRS design process. Objective function (1.2) and constraints (2) and (3) are used to select service units from among the candidate service units and assign requests to them in accordance with the spatiotemporal features of the requests. The bottom of Figure 2 shows the  figure, where the cycles represent requests (the diagonal striped circles represent origins, and the striped circles represent destination); the squares represent boarding parts and alighting parts of service units (the diagonal striped squares represent selected boarding parts, the striped squares represent selected alighting parts, the filled squares represent the candidate boarding parts, and the blank squares represent the candidate alighting parts). Both the boarding part and the alighting part include information about location and time window. The figure shows the process of design service units and service routes from requests that exists in the optimization model. In the middle and the bottom of the figure, requests are assigned to service units according to spatiotemporal features of the requests and constraints. Then, in the top of the figure, the boarding parts and alighting parts of service units are assigned to service routes and to be served by the vehicles in the fleet.

Clustering-First, Routing-Second Algorithm for Solving the LIPPD
The solution process for the LIPPD is the process of finding the optimal service units and service routes in accordance with the spatiotemporal features of the requests. Ideally, we wish to find an optimal solution to the LIPPD that considers the trade-off between passenger costs and vehicle costs [53]. However, from the model, we can see that it is not easy to obtain solutions to this problem: the objective functions are multitarget functions with nonlinear structures. Therefore, instead of looking for the globally optimal solution, we use a two-stage algorithm to find a near-optimal solution. This section introduces the clustering-first (DSU), routing-second (FAR) solution approach. This solution framework implies the following assumption about the priority of the targets: the burden placed on the passengers by this new service mode is minimized first, and the operational cost is then minimized. In other words, we pay more attention to the requesters' costs than to the operators' costs.
To facilitate clearer description throughout the rest of the paper, we use RM to denote the submodel comprising Equations (1a), (4)- (17), (19) and (20) and use CM to denote the submodel comprising Equations (1b), (2), (3) and (18). The DSU process focuses on solving CM, and the FAR process focuses on solving RM.

Preprocessing
The preprocessing phase consists of two tasks: (i) the detection and deletion of self-contradictory requests and (ii) the detection and deletion of abnormal requests.
Since the service considered here needs to be feasible in an actual scenario, it is necessary to exclude user-submitted requests with self-contradictions. In general, such requests are not spatiotemporally feasible; for example, the separation between the expected time windows for boarding and alighting may be too small, meaning that a vehicle cannot complete the specified task, while obeying the speed limit in the corresponding road network. Such self-contradictory requests should be deleted to avoid infeasible solutions to the model. For this purpose, we use formula (21) to detect self-contradictory requests: where ||v i || Am ∀v i ∈ Cad is a norm defined as ||v i || Am ||l i+ , l i− || lt R − (e i− − e i+ ), with ς representing the operation time redundancy for a vehicle and a, b lt R representing the vehicle travel time (obtained from a digital map) between locations a and b. The definition of the norm ||a, b|| lt R will be discussed in Section 2.2.2.
In addition to self-contradictory requests, abnormal requests should also be detected and removed from the request pool. There are two main types of abnormal requests: (i) requests caused by misoperation, such as requests whose starting and ending locations are too close or repeated requests with the same ID number, and (ii) spatially invalid requests, such as requests corresponding to locations that are significantly beyond the scope of the service. The screening for such abnormal requests relies mainly on simple heuristic rules.

Distance Measurement
As a service operated in a metropolitan context, LIPDRS requires a reasonably accurate distance calculation method. The temporal and spatial distances between two points cannot be directly calculated using the Euclidean distance but should instead be calculated in accordance with the conditions of the road network, speed limits, etc. Measuring the distance by means of the path planning function of an electronic map is a viable option. Here, we define two temporal distance norms, In addition, we define ||x i − v j || At W as a measure of the temporal distance between a request and a service unit, as expressed by Equation (3). This measure is vitally important to the process of assigning requests to service units: where ||(a, b)|| lt W represents the passenger walking time (obtained from a digital map) between locations a and b and ||(a, b)|| twt represents the time difference between time windows a and b.

Main DSU Flows
We introduce the DSU solution process in the form of the flow chart shown in Figure 3. During the DSU solution process, on one hand, the locations and time windows of the service units must be selected to allow the requests to be assigned to service units, and on the other hand, as in other clustering-based problems, determining the number of clusters c is also an important task [7]. Therefore, a cluster number determination method is required. To this end, the structure of this flow chart is composed of two loops and corresponding initialization modules. The outer loop is designed to determine the optimal number of service units, and the inner loop applies a clustering method to determine the service units in accordance with a given number of clusters. calculation method. The temporal and spatial distances between two points cannot be directly calculated using the Euclidean distance but should instead be calculated in accordance with the conditions of the road network, speed limits, etc. Measuring the distance by means of the path planning function of an electronic map is a viable option. Here, we define two temporal distance norms, , and , , and one spatial distance norm, , , where , and , are obtained from a digital map (e.g., AMAP).
In addition, we define − as a measure of the temporal distance between a request and a service unit, as expressed by Equation (3). This measure is vitally important to the process of assigning requests to service units: where ‖ , ‖ represents the passenger walking time (obtained from a digital map) between locations a and b and ‖ , ‖ represents the time difference between time windows and .

Main DSU Flows
We introduce the DSU solution process in the form of the flow chart shown in Figure 3. During the DSU solution process, on one hand, the locations and time windows of the service units must be selected to allow the requests to be assigned to service units, and on the other hand, as in other clustering-based problems, determining the number of clusters is also an important task [7]. Therefore, a cluster number determination method is required. To this end, the structure of this flow chart is composed of two loops and corresponding initialization modules. The outer loop is designed to determine the optimal number of service units, and the inner loop applies a clustering method to determine the service units in accordance with a given number of clusters.  In Figure 3, modules (3)-(7) compose the inner loop for the clustering task. The purpose of clustering, an unsupervised machine learning method, is to aggregate the requesters by generating centroids (which can be used as the basis for selecting service units) based on similarities among the spatiotemporal features of the requests. Among the many clustering methods available, we choose fuzzy c-means clustering (FCM) as the basic method of identifying service units and assigning requests. FCM is a soft clustering method. Unlike hard clustering methods (e.g., k-means), FCM simultaneously considers the possibilities that a data point may belong to multiple clusters by means of membership degrees, resulting in a much finer degree of data modeling that yields numerical results that can be used for further analysis [7]. However, the original FCM method cannot accommodate the walking time constraint, which is vitally important to our research. We therefore propose a new similarity measure to modify the original FCM to consider a limit on the walking time. In Figure 3, modules (3)-(7) compose the inner loop for the clustering task. The purpose of clustering, an unsupervised machine learning method, is to aggregate the requesters by generating centroids (which can be used as the basis for selecting service units) based on similarities among the spatiotemporal features of the requests. Among the many clustering methods available, we choose fuzzy c-means clustering (FCM) as the basic method of identifying service units and assigning requests. FCM is a soft clustering method. Unlike hard clustering methods (e.g., k-means), FCM simultaneously considers the possibilities that a data point may belong to multiple clusters by means of membership degrees, resulting in a much finer degree of data modeling that yields numerical results that can be used for further analysis [7]. However, the original FCM method cannot accommodate the walking time constraint, which is vitally important to our research. We therefore propose a new similarity measure to modify the original FCM to consider a limit on the walking time.
The outer loop, which is composed of the inner loop and modules (1), (2), (8) and (9), selects a number of candidate clusters from a set and passes this cluster number to the inner loop. This is an enumeration process. At the end of the loop, the number of clusters that results in the best objective function value is selected as the final number of service units. The number of service units impacts the efficiency, effectiveness and equity of passenger service and vehicle operation. If the service units are too densely distributed (meaning that the clusters are small in size), the speed of service will not be effectively accelerated; if the service units are too sparsely distributed, more requesters will have to walk farther. Thus, the number of service units is a key factor in striking this balance. We will introduce each module in detail throughout the rest of Section 2.3.

Initialization
There are two initialization modules in Figure 3. Module (1) is for the initialization of the set of the candidate numbers of clusters, cc, wherein we construct a set-in increments of one within a given range. The other initialization module is module (3). The elements in the matrix of membership values U n×c for FCM are randomly initialized in accordance with the following constraints: the elements u ji must satisfy c−1 j=1 u ji = 1 and 0 ≤ u ji ≤ 1.

Measurement of Similarity between a Request and a Service Unit
The norm for measuring the similarity of the jth request to the centroid of cluster i is vitally important because it directly influences the objective function (1.2). In the flow of the DSU solution process, the similarities are calculated in module (5). A typical norm is the A-norm proposed by Bezdek et al. [54]. However, the A-norm can measure only geographical proximity and thus is not suitable for measuring spatiotemporal proximity, as is necessary for the walking time constraint. Therefore, we propose a Gaussian-kernel-based norm that allows the maximum walking time to be considered, as presented in (23): In (23), we use the temporal distance ||r j − v x || At W to represent the distance between two locations instead of calculating the Euclidean distance. The definition of ||r j − v x || At W can be found in Section 2.2.2. A Gaussian kernel, a kernel type that is widely used in machine learning [55], can make objective function (1.2) valid only in a local area (i.e., if a requester's walking time is greater than the maximum passenger walking time mwt, then the similarity to the corresponding service unit will be very low). However, as a result of the mechanism of generating the initial solution for FCM [7,56], the initially selected service units might be far away from any request. In other words, a service unit might have minimal similarity to all requests. Therefore, we design (23) as a piecewise function to avoid empty clusters (i.e., when there is at least one request in a service unit, the Gaussian kernel is used). In (23), σ is a width parameter that controls the range in the radial direction, whereby we can control mwt.

Iteration Method and Termination Criteria
In modules (4) and (6), the service units v i are reselected, and the membership values u ji are updated. For this purpose, we propose Equations (24) and (25). The norm a mCad can be used to enumerate the candidate service units and find the nearest one to the centroid a i of the requests, where the centroid a i is calculated as a i = n−1 j=0 u ji m r j / n−1 j=0 u ji m .
Following Bezdek [56], m = 2 is adopted in this study. Two criteria for terminating the iterative process are designed for module (7). One is a threshold ε. If U p+1 − U p < ε, then the iterative process terminates, where U p is the membership matrix in the pth iteration. The other is a maximum number of iterations, mi. After either termination criterion is triggered, any remaining requests with very low similarities (close to 0) to any service unit become isolated requests.

Routing (Second Stage): Fleet Assigning and Routing (FAR)
We develop a model based on the pickup and delivery problem with time windows (PDPTW) for formulating the FAR problem. Considering the complexity of the search space in a large-scale network, we use a column generation (CG) algorithm in the framework for solving the model. In addition, to address the influence of variable relaxation, a branching scheme is used to generate a binary solution tree to enable the selection of a feasible and optimal solution.

Main FAR Flows
The FAR process converts the service units determined via the DSU process into service routes, which represent the specific times and order of the service units served by the fleet. We introduce the FAR solution process in the form of the flow chart shown in Figure 4. The structure of this flow chart is composed of two loops and corresponding initialization modules. These two loops are designed to calculate the model and to guarantee that there is no infeasible solution, namely, no noninteger solution, which can occur as a side effect of the relaxation of the model.
In Figure 4, the inner loop (denoted by A) applies the CG algorithm, iteratively solving the master problem and the subproblem. This loop is composed of modules (4)- (9). The outer loop (denoted by B) applies a branching scheme to avoid noninteger solutions. It is composed of loop A and modules (10)- (15). Modules (1)-(3) are the initialization modules for the two loops. Focusing on these two loops and the modules they include, we will elaborate on the details of the algorithm in Sections 2.4.2-2.4.6.
Following Bezdek [56], = 2 is adopted in this study. Two criteria for terminating the iterative process are designed for module (7). One is a threshold . If − < , then the iterative process terminates, where is the membership matrix in the iteration. The other is a maximum number of iterations, . After either termination criterion is triggered, any remaining requests with very low similarities (close to 0) to any service unit become isolated requests.

Routing (Second Stage): Fleet Assigning and Routing (FAR)
We develop a model based on the pickup and delivery problem with time windows (PDPTW) for formulating the FAR problem. Considering the complexity of the search space in a large-scale network, we use a column generation (CG) algorithm in the framework for solving the model. In addition, to address the influence of variable relaxation, a branching scheme is used to generate a binary solution tree to enable the selection of a feasible and optimal solution.

Main FAR Flows
The FAR process converts the service units determined via the DSU process into service routes, which represent the specific times and order of the service units served by the fleet. We introduce the FAR solution process in the form of the flow chart shown in Figure 4. The structure of this flow chart is composed of two loops and corresponding initialization modules. These two loops are designed to calculate the model and to guarantee that there is no infeasible solution, namely, no noninteger solution, which can occur as a side effect of the relaxation of the model.
In Figure 4, the inner loop (denoted by A) applies the CG algorithm, iteratively solving the master problem and the subproblem. This loop is composed of modules (4)- (9). The outer loop (denoted by B) applies a branching scheme to avoid noninteger solutions. It is composed of loop A and modules (10)- (15). Modules (1)-(3) are the initialization modules for the two loops. Focusing on these two loops and the modules they include, we will elaborate on the details of the algorithm in B A Figure 4. FAR Flows.

Initialization
(1)-(2): Dantzig-Wolfe Decomposition and Relaxation of the Master Problem To circumvent the difficulties facing such large-scale problems [57], we decompose the model introduced in Section 2.1.2 via Dantzig-Wolfe (DW) decomposition. After DW decomposition, the original model has been decomposed into an equivalent combination of two smaller problems: a master problem (the fleet assignment problem) and a subproblem (the routing and pricing problem).
The master problem is the set coverage problem, and the subproblem is the route planning and pricing problem. The decomposed model can be solved via a CG algorithm. To obtain the marginal cost of the master problem, we relax the integer decision variables of the master problem to positive variables.
(3): Solution Tree A side effect of the relaxation of the master problem is that the route selection process may yield noninteger solutions. A mature approach for overcoming this problem of noninteger solutions is the branching approach [57]. To obtain information about when to start and end each branch, we design a solution tree with index L and a feasibility judgment method to store solutions and judge whether the process should terminate or continue.
As shown in Figure 4, the solution tree is a layered binary tree. The original solution is stored as the root. The feasibility judgment method is used only to judge the solutions in the Lth layer. If the feasibility judgment method detects that a noninteger route has been selected in the solution, then the branching method described in [57] is applied. Otherwise, no noninteger route has been selected in the solution, and the algorithm ends. Finally, we choose the best solution in the solution tree as the output.

Master Problem (Relaxed Set Partitioning Problem, MP)
The master problem MP is a set partitioning problem, which can be formulated as follows. Objective: Minimize r k ∈Ω c k θ k + i∈N P 1 − r k ∈Ω δ k i θ k = i∈N P + r k ∈Ω θ k c k − i∈N Pδ k i (26) subject to Here, M denotes the set of available vehicles, Ω denotes the set of all feasible pickup and drop-off routes (r k ) served by vehicle k, and c k denotes the cost of route r k .
The primary objective of MP is to minimize the number of vehicles needed, and the secondary objective is to minimize the total cost. Therefore, c k is calculated as c k = FC + VCL k , where L k denotes the length of route r k , VC denotes the unit cost of a vehicle, and FC denotes the fixed cost of a vehicle, which should be relatively large. This cost structure can be achieved by defining the travel costs c ij = d ij (i De) and c Dej = FC + d Dej ( j De).
Furthermore, let parameters δ k i be defined as follows: Let the θ k be positive variables of the master problem: If θ k = 1, then route r k ∈ Ω is selected, and if θ k = 0, then route r k ∈ Ω is not selected; otherwise, if θ k ∈ (0, 1), the solution is infeasible. The process of replacing the θ k with positive variables rather than integer variables is called relaxation, which is performed in module (2). Equation (27) ensures that each service unit will be either served or canceled, and the penalty for canceling service units is represented by i∈N P 1 − r k ∈Ω δ k i θ k in objective function (26). Equation (28) ensures that the number of selected routes cannot exceed the scale of the fleet.

Subproblem (Routing and Pricing Problem, S k )
Let the graph of the subproblem be denoted by G = (N, E), and let the cost of each arc (i, j) ∈ E be denoted by c G ij . Then, the subproblem S k is formulated as follows: Objective: Minimize i∈V j∈V c G ij x ij (32) subject to i∈V,j∈V x ij d ij ≤ d max (43) i∈V,j∈V The target of S k is to find routes that minimize the reduced cost for the master problem. In objective function (32) and w belong to → u , w , which is the marginal cost of the master problem MP. In S k , the following variables are significant: z i (i ∈ N), a binary variable, is equal to 1 if service unit i ∈ N is served by the vehicle and is equal to 0 otherwise. x ij ((i, j) ∈ V × V), a binary variable, is equal to 1 if the vehicle travels from location i ∈ V to location j ∈ V and is equal to 0 otherwise. D i (i ∈ V), an integer variable, denotes the departure time at location i ∈ V. y i (i ∈ V), an integer variable, denotes the load of the vehicle when arriving at location i ∈ V.
Regarding the constraints, Equations (33) and (34) ensure that the vehicle visits either both the boarding and alighting points of a service unit or neither of them. Equations (35) and (36) ensure that the vehicle starts and ends at its depot. Equations (37)- (39) and (45) together form the precedence constraints. Equations (40)- (42) and (46) together form the capacity constraints. Equations (43) and (44) ensure that the vehicle will return within a limited travel distance and duration.

Iteration Termination Criteria
The solution process for MP and S k is embedded in loop A. After the input data and initial routes are obtained, MP is solved first, yielding a marginal cost → u , w for S k . Then, S k is solved using a labeling method based on dynamic programming [58], and the solution is a new route for MP. The new route is the one that incurs the minimal reduced cost for MP. This process does not stop until the solution to S k meets the following termination criterion: the reduced cost of the returned route is equal to or greater than 0.
The individual termination criteria of the two loops are also important. In module (8), if the reduced cost obtained from the subproblem is greater than 0, then the loop stops. Otherwise, the loop continues. In module (11), if the solutions in the Lth layer are feasible, then the loop stops. Otherwise, the loop continues.

Branching Scheme
When a solution appears infeasible, namely, noninteger in nature, a branching scheme should be applied, as in modules (12)- (14). The rules of the branching scheme are as follows: (i) For each arc (i, j), calculate the values ) and select an arc that satisfies 0 < f ij < 1.
(ii) Derive a branch on which the selected arc cannot belong to the solution (enforce f ij = 0 and delete the corresponding edges and routes).
(iii) Derive a branch on which the arc must be included in the solution (enforce f ij = 1 and delete the corresponding edges and routes).
After branching, the index L is updated, and the algorithm flow continues until the termination criterion is satisfied.

Evaluation Metrics
Evaluation metrics are needed to (i) provide a basis for the comparison of different transport options with different objective functions and (ii) evaluate whether a plan is economically feasible, i.e., whether passengers will be economically willing to accept this travel option.
To evaluate the proposed transport option, the difference between it and other existing methods must be quantified. However, the objective functions of different models are different; therefore, it is not possible to directly compare the quality of different service methods based on their objective function values. Here, we use several indicators to evaluate different modes of transport. Economic indicators are often the best basis for comparison. This article proposes indicators for evaluating services from two perspectives: (i) fleet operating costs and (ii) passenger travel costs. The former perspective is the focus of the operators, and the latter is the focus of both the operators and the passengers and can be considered to directly determine whether a particular transport option is acceptable to passengers. These indicators are related to the objective function but provide more detailed information about the solution.
Here, we define two sets of notations for evaluation metrics. From the perspective of the passenger travel costs, we define (i) C PW : the cost that passengers expect to incur due to additional walking time (CNY/min); (ii) T TW : the total walking time (min); (iii) D TPV : the total distance traveled by all passengers on vehicles (km); (iv) C G : the unified generalized cost (in terms of time and comfort) of passengers using the service; (v) T TPV : the total time spend by all passengers on vehicles; (vi) P F : the fixed price for using the service; (vii) P V : the unit variable price for using the service; and (viii) N S : the number of requests served. From the perspective of the fleet operating costs, we define (ix) N: the number of vehicles used; (x) D TTV : the total distance traveled by vehicles (km); (xi) C F : the unit fixed cost of using a vehicle; (xii) C V : the unit variable cost of using a vehicle; (xiii) P RR : the penalty cost for a rejecting request; and (xiv) N RR : the number of rejected requests. Then, we calculate two costs for evaluation, as shown below.
α. The passenger cost, which is composed of (i) the walking cost c tw , c tw = C PW ·T TW ; (ii) the ticket cost c tt , c tt = P F ·N S + P V ·D TPV ; and (iii) the general cost c gt incurred due to travel time, c gt = C G ·T TPV . The formula for α is given in Equation (49): β. The fleet cost, which is composed of (i) the fixed cost c f f of the fleet, c f f = N·C F ; (ii) the variable cost c v f of the fleet, c v f = VC V ·D TTV ; and (iii) the rejection penalty p r , p r = N RR ·P RR . The formula for β is given in Equation (50): Intuitively, by comparing α and β for different transport options, we can judge which option is economically superior.

Data Preparation
Consistent with the process discussed in Section 2.1.1, we used records from a booking system for collecting real-world requests as the basic dataset (RRD) for our numerical experiment. We focused on requests located in northwest Beijing, that is, in the Changping and Haidian districts, because obvious LIPP phenomena arise in these districts as a result of their functional differences: Changping consists of large-scale residential areas, whereas Haidian consists of many business, educational and scientific areas. Many people work in Haidian and live in Changping; therefore, they require transportation between these two districts. The dataset used in this experiment contains 712 requests scattered throughout northwest Beijing during the period of 6:00-9:00 (the morning peak period) on 2019/04/21.
To further explore the performances of the proposed option and other options, we would like to consider different types of datasets (different request distributions) to test the application scopes of the different models. However, it would be difficult to obtain arbitrary request distributions without forcing users to make particular choices. Therefore, we artificially generated a series of datasets based on RRD (i.e., we simulated the booking process by randomly generating requests based on RRD in accordance with a series of rules). By using these datasets, we could test the model performances in various contexts. We defined the following artificial datasets: (i) ARR, a dataset of artificial random requests, contains 700 artificially generated requests that are evenly and randomly distributed in time and space. (ii) ACR, a dataset of artificial clustered requests, contains 700 artificially generated requests separated into 50 clusters that are relatively concentrated in time and space. (iii) ARCR25%, a dataset consisting of a mix of artificial requests with random and clustered structures, contains 700 artificially generated requests, 25% of which are relatively concentrated in time and space and the remaining 75% of which are not. (iv) ARCR50%, another dataset consisting of a mix of artificial requests with random and clustered structures, contains 700 artificially generated requests, 50% of which are relatively concentrated in time and space and the remaining 50% of which are not. (v) ARCR75%, the final dataset consisting of a mix of artificial requests with random and clustered structures, contains 700 artificially generated requests, 75% of which are relatively concentrated in time and space and the remaining 25% of which are not. The preprocessing method discussed in Section 2.2.1 was applied to obtain requests that were suitable for scheduling. Then, we conducted our experiment, as described in Section 3.2.

Experimental Scheme
The numerical experiment consisted of four steps. In step (i), we used the proposed clustering algorithm to determine the optimal service units for the real-world dataset RRD. Then, we compared the performance of the proposed clustering method with that of the classical k-means method. In step (ii), we compared the performance of the proposed routing algorithm with that of GAMS (a commercial modeling and solution platform) based on cases with 13-25 service units. In step (iii), we evaluated the performance of the proposed LIPDRS on real requests and analyzed the request rejection situation as well as the economic feasibility on the real-world dataset RRD in comparison to the existing door-to-door transport option. In addition, by combining the results of this step with those of step (i), we selected the optimal plan for RRD. In step (iv), we extended the experiment to an evaluation of the performances of the proposed transport option and the existing door-to-door option on the various artificial request datasets to determine the scopes of application of the different options.
The parameters that needed to be prespecified for this experiment are listed in Table 2. In addition, we used information on the real road network in Beijing, as obtained from 'AMAP', to obtain distance measurements using the method discussed in Section 2.1.3.  2 We used a currently operating bus hub as the depot. 3 Expected speed in the road network (ESRN), obtained from 'AMAP', a digital map. 4 In the context of our study, we regard the travel distance as being proportional to the travel cost. However, the proportionality between travel distance and cost is not a key element of our study. Therefore, we use a conceptual unit to bypass the proportion determination problem, following reference [57]. According to survey data, the average fuel consumption of a large bus is approximately 0.3-0.35 L/km. Based on the price of diesel in Beijing, one unit corresponds to approximately 2 CNY.

Results and Discussion of Step (i)
The objective function values L(c) (Equation (1b)) corresponding to different numbers c of service units were calculated in this step via enumeration for two clustering methods (the proposed fuzzy clustering method and the classical k-means hard clustering method), and the results are shown in Figure 5. For both the proposed method and the classical method, the L(c) values are high for c values of {36, 42, 44, 51, 53} and reach a peak when c is set to 53.
Additionally, the performance of k-means is lower than that of the proposed fuzzy clustering method. For example, when c is set to 53, the L(c) value (Equation (1b)) obtained with the proposed clustering method is 8% higher than that obtained with k-means. The difference is even larger (approximately 2%) when c is set to 40. These findings suggest that the proposed soft clustering method is superior to the hard-clustering method.
We selected the number of clusters resulting in the highest L(c), as an example to illustrate the distribution of the service units, as shown in Figure 6. However, this is not necessarily the most appropriate number of service units. After calculating all plans in their entirety, including the service units and service routes, we can evaluate the statistical characteristics of each plan based on the metrics discussed in Section 2.5 to select the optimal number of service units. In our discussion of step (iii), we will further consider the optimal plan for RRD.

Results and Discussion of Step (i)
The objective function values (c) (Equation (1.2)) corresponding to different numbers of service units were calculated in this step via enumeration for two clustering methods (the proposed fuzzy clustering method and the classical k-means hard clustering method), and the results are shown in Figure 5. For both the proposed method and the classical method, the ( ) values are high for values of {36, 42, 44, 51, 53} and reach a peak when is set to 53.
Additionally, the performance of k-means is lower than that of the proposed fuzzy clustering method. For example, when is set to 53 , the (c) value (Equation (1.2)) obtained with the proposed clustering method is 8% higher than that obtained with k-means. The difference is even larger (approximately 2%) when is set to 40. These findings suggest that the proposed soft clustering method is superior to the hard-clustering method.
We selected the number of clusters resulting in the highest (c), as an example to illustrate the distribution of the service units, as shown in Figure 6. However, this is not necessarily the most appropriate number of service units. After calculating all plans in their entirety, including the service units and service routes, we can evaluate the statistical characteristics of each plan based on the metrics discussed in Section 2.5.2 to select the optimal number of service units. In our discussion of step (iii), we will further consider the optimal plan for RRD.

Results and Discussion of Step (ii)
To validate the proposed routing method, problem instances with 13-55 service units were solved to obtain two types of solutions, SOGs and SCGs, where the SOG for a particular problem instance is the solution to RM (discussed in Section 2.1.3) as obtained by GAMS (a commercial solution software), whereas the SCG is the solution to the models MP and S k obtained using the CG algorithm. We regard the SOG as the standard solution to the problem. By comparing the SOGs and SCGs in terms of their objective values and calculation times, we can evaluate whether the proposed method is feasible and efficient. Figure 7 presents a comparison of the SOGs and SCGs from the perspectives of the calculation times and objective values. The proposed routing method obtains the same results as GAMS does; thus, the proposed method is feasible. Regarding the calculation time, as the number of service units increases, the calculation times for the SCGs exhibit acceptable linear growth when the number of service units is less than 35. As the number of service units increases greater than 35, the calculation time increases significantly. When the number of service units reaches 55, the calculation time is 21,612 s (approximately 6 h), which is essentially the maximum acceptable limit. In contrast, the calculation times for the SOGs indicate that GAMS cannot address large-scale problems (involving approximately 50 service units). With more than 21 service units, the SOG calculation times exhibit much more rapid growth than the SCG calculation times do; when the number of service units is 39, the calculation time

Results and Discussion of Step (ii)
To validate the proposed routing method, problem instances with 13-55 service units were solved to obtain two types of solutions, SOGs and SCGs, where the SOG for a particular problem instance is the solution to (discussed in Section 2.1.3) as obtained by GAMS (a commercial solution software), whereas the SCG is the solution to the models and obtained using the CG algorithm. We regard the SOG as the standard solution to the problem. By comparing the SOGs and SCGs in terms of their objective values and calculation times, we can evaluate whether the proposed method is feasible and efficient. Figure 7 presents a comparison of the SOGs and SCGs from the perspectives of the calculation times and objective values. The proposed routing method obtains the same results as GAMS does; thus, the proposed method is feasible. Regarding the calculation time, as the number of service units increases, the calculation times for the SCGs exhibit acceptable linear growth when the number of service units is less than 35. As the number of service units increases greater than 35, the calculation time increases significantly. When the number of service units reaches 55, the calculation time is 21,612 s (approximately 6 h), which is essentially the maximum acceptable limit. In contrast, the calculation times for the SOGs indicate that GAMS cannot address large-scale problems (involving approximately 50 service units). With more than 21 service units, the SOG calculation times exhibit much more rapid growth than the SCG calculation times do; when the number of service units is 39,   Now, let us consider a case with 13 service units as an example to illustrate how the branching scheme works during the FAR solution process. The edges in this case are indexed by numbers from 0 to 156. The relationship between the objective value and the iteration number throughout the entire solution process is shown in Figure 8. In the 15th iteration, the termination criterion is met. The objective value of the solution is 1097.4995 (with 10 vehicles), and the calculation time is 14,230 ms. Now, let us consider a case with 13 service units as an example to illustrate how the branching scheme works during the FAR solution process. The edges in this case are indexed by numbers from 0 to 156. The relationship between the objective value and the iteration number throughout the entire solution process is shown in Figure 8. In the 15 th iteration, the termination criterion is met. The objective value of the solution is 1097.4995 (with 10 vehicles), and the calculation time is 14,230 ms.  Table 3 presents the details of the original solution obtained with the proposed algorithm. Routes 14, 17, and 22 are fractional; therefore, the solution is not feasible. Consequently, the branching scheme [57] must be applied to guarantee the feasibility of the final solution.   Table 3 presents the details of the original solution obtained with the proposed algorithm. Routes 14, 17, and 22 are fractional; therefore, the solution is not feasible. Consequently, the branching scheme [57] must be applied to guarantee the feasibility of the final solution. We choose edge 101 as the branching point. After branching, we obtain a left branch and a right branch. The left branch eliminates edges 162, 5, 23, 56, 70, 80, 162, and 193 as well as routes including those edges. By contrast, the right branch forces edge 101 to be selected. From the objective values, we find that the left branch (1139.103) is superior to the right branch (1157.919). The objective value of the left branch is equal to that of the SOG (1139.103).

Results and Discussion of Step (iii)
Based on the various service units generated in step (i), we further planned the routes and evaluated the results using the metrics proposed in Section 2.5 to choose the optimal service plan. The results based on different numbers of service units (c = {50, 51, 52, 53, 54}) are shown in Table 4.
From the perspective of the total passenger cost, as the number of service units increases, the walking time of passengers (and the corresponding generalized cost, c tw ) gradually decreases. This trend is also roughly reflected in the passengers' travel time (and the corresponding ticket cost, c tt ) and ride time (and the corresponding general cost due to travel time, c gt ), but the following difference is found: when the number of service units is 53, the passengers spend the least time and distance traveling (the corresponding c tt and c gt values are the smallest). The reason for this phenomenon is that these two indicators are jointly determined by the distribution of the service units and the designed routes. In a real-world road network, even a slight change in the service units can result in a change in the routes and, thus, the costs.
From the perspective of the total fleet cost, when the number of service units is in the range of 50-54, the number of rejected requests is not high; in fact, only one request is rejected (p r = 180). When the number of service units is 50-53, the number of vehicles used remains almost constant, being either 19 or 20 (c f f = 1710 or 1800). When the number of clusters is 54, however, the number of vehicles used increases to 23 (c f f = 2070), which is much greater than the numbers of vehicles used in the other cases. A similar trend is seen in terms of the variable cost (c v f ), and when the number of service units is 53, the cost is optimal. The reason for this finding may be that the routes corresponding to the specific distribution of the service units in the real-world road network are relatively simple when the number of service units is 53.
In summary, for RRD, the best plan is obtained when the number of service units is 53. In this plan, requesters must walk for approximately 6.8 min (c tw = 34.00) on average, and approximately 13.43 requests are served per service unit. 1 c tw Walking cost, c tw = C PW ·T TW , where C PW represents the cost that passengers expect to incur for additional walking time (CNY/min) and T TW represents the total walking time (min). 2 c tt : Ticket cost, c tt = P F ·N S + P V ·D TPV , where P V is the unit variable price for using the service, P F is the fixed price for using the service, N S is the number of requests served, and D TPV represents the total distance traveled by all passengers in vehicles (km). 3 c gt : General cost due to travel time, c gt = C G ·T TPV , where C G represents the unified generalized cost (in terms of time and comfort) of passengers using the service and T TPV represents the total time spent by all passengers in vehicles. 4 c f f : Fixed cost of the fleet, c f f = N·C F , where N is the number of vehicles used and C F is the unit fixed cost of using a vehicle. 5 c v f : Variable cost of the fleet, c v f = C V ·D TTV , where C V is the unit variable cost of using a vehicle and D TTV represents the total distance traveled by vehicles (km). 6 p r : Rejection penalty, p r = N RR ·P RR , where N RR is the number of rejected requests and P RR is the penalty cost for rejecting a request.

Results and Discussion of Step (iv)
In this step, we compared the proposed LIPDRS with the existing door-to-door service (based on the model presented in [57]) on various datasets (discussed at the beginning of Section 3). The results are reported in Tables 5 and 6.
where P V is the unit variable price for using the service, P F is the fixed price for using the service, N S is the number of requests served, and D TPV represents the total distance traveled by all passengers in vehicles (km). 3 c gt : General cost due to travel time, c gt = C G ·T TPV , where C G represents the unified generalized cost (in terms of time and comfort) of passengers using the service and T TPV represents the total time spent by all passengers in vehicles. 4 c f f : Fixed cost of the fleet, c f f = N·C F , where N is the number of vehicles used and C F is the unit fixed cost of using a vehicle. 5 c v f : Variable cost of the fleet, c v f = C V ·D TTV , where C V is the unit variable cost of using a vehicle and D TTV represents the total distance traveled by vehicles (km). 6 p r : Rejection penalty, p r = N RR ·P RR , where N RR is the number of rejected requests and P RR is the penalty cost for rejecting a request.
where P V is the unit variable price for using the service, P F is the fixed price for using the service, N S is the number of requests served, and D TPV represents the total distance traveled by all passengers in vehicles (km). 3 c gt : General cost due to travel time, c gt = C G ·T TPV , where C G represents the unified generalized cost (in terms of time and comfort) of passengers using the service and T TPV represents the total time spent by all passengers in vehicles. 4 c f f : Fixed cost of the fleet, c f f = N·C F , where N is the number of vehicles used and C F is the unit fixed cost of using a vehicle. 5 c v f : Variable cost of the fleet, c v f = C V ·D TTV , where C V is the unit variable cost of using a vehicle and D TTV represents the total distance traveled by vehicles (km). 6 p r : Rejection penalty, p r = N RR ·P RR , where N RR is the number of rejected requests and P RR is the penalty cost for rejecting a request.
First, building on the results of steps (i) and (iii), we compared the plans designed for the LIPDRS and door-to-door (DtD) options. Regarding the passenger cost, LIPDRS results in a cost of approximately 160,723.46 (= 17409.55 + 4763.66 + 138550.25) CNY, and DtD results in a cost of approximately 269,556.67 (= 0 + 8959.04 + 260597.63) CNY. Obviously, from the passenger's point of view, LIPDRS costs 40.37% less than DtD. The difference in the general cost due to travel time is the main reason. To meet the LIPP in the real-world road network, the DtD service will require passengers to spend more time on vehicles after boarding, and the corresponding cumulative waiting time will constitute a significant difference in the passengers' travel expenses. Regarding the fleet cost, LIPDRS costs approximately 4369.43 units, and DtD costs approximately 6804.38 units. Thus, LIPDRS costs approximately 35.79% less than DtD. This difference is due to the establishment of the service units, which greatly reduces the transportation burden placed on the vehicles and thus reduces operating costs. Therefore, in this real-world case, the proposed transport option has an economic advantage.
Next, we analyzed the scope of application of the proposed solution by comparing the two different transport options on the artificial datasets. From ARR to ACR (including ARCR25-75%), the spatiotemporal concentration of requests gradually increases. For LIPDRS, the walking time of passengers (and the corresponding generalized cost, c tw ) is relatively high on ARR. In addition, the c tt and c gt values are high, almost the same as those for the DtD option. In contrast, DtD does not require passengers to walk. Consequently, the total passenger cost (259,834.19 CNY) is lower than that for LIPDRS (314,414.32 CNY). Moreover, LIPDRS also costs more than DtD from the perspective of the fleet cost. These results indicate that LIPDRS lacks adaptability to random and decentralized distributed requests such as those in ARR.
As the level of demand centralization increases (i.e., as the dataset changes from ARR to ARCR25%, ARCR50% and ARCR75%), the c tw value of LIPDRS decreases, and the c tt and c gt values show the same trend. In terms of the fleet cost, LIPDRS begins to cost less than DtD with ARCR50%. At this point, the advantage of LIPDRS begins to emerge. When the ratio of requests that are spatiotemporally concentrated is above 50%, LIPDRS is the more suitable option. Therefore, we believe that the proposed LIPDRS is much more suitable for accommodating the LIPP in a metropolitan context. For example, on ARCR75%, LIPDRS can reduce the cost for passengers by approximately 17.94% and the cost for operators by approximately 32.14%.

Conclusions
To alleviate the problems caused by LIPP in metropolises, this paper proposes an innovative type of demand-responsive public transport service called LIPDRS. By providing an attractive transportation option for passengers who typically drive or use conventional public transport, this type of demand-responsive service can potentially reduce car usage and expand passenger accessibility in congested metropolitan areas.
To obtain the optimal operation plan for LIPDRS, we model the LIPDRS decision process using a mixed-integer programming model with dual objectives to optimize the service from the perspectives of both passengers and operators. This model includes a carefully designed passenger walking time constraint and request rejection mechanism to ensure that the resulting service plan can guarantee not only good quality of service but also improved operational efficiency. Additionally, this paper proposes a clustering-first, routing-second framework for solving the proposed model. A modified version of the fuzzy c-means clustering algorithm and a column generation algorithm are utilized in this framework to determine the service units and perform vehicle assignment and routing, respectively. In addition, to overcome the problems caused by relaxing the model variables, a branching scheme is proposed.
Based on the designed methodology, a four-step experiment was conducted to investigate different transport options on various types of datasets. The results demonstrate the advantages of the proposed method and illustrate its application scope. The contributions of this paper can be summarized as follows. (i) We propose a new transport option along with a complete system of modeling and solution methods. A maximum passenger walking time constraint, a request rejection mechanism and a scheme for ensuring solution feasibility are all considered in this scheduling system. (ii) The proposed transport option has been proven to achieve better performance (by 40.37% for passengers and by 35.79% for operators) than the existing "door-to-door" transport option does when faced with LIPP. (iii) By testing on different datasets, we prove that the proposed service is more suitable for the request distributions that are spatiotemporally concentrated. We have found that when the ratio of requests that are spatiotemporally concentrated is above 50%, LIPDRS is the more suitable option. (iv) The performance advantages of the proposed clustering-first, routing-second solving framework have been proven through experiments. Regarding the individual components of the framework, the proposed soft clustering algorithm has been proven to be superior to the k-means hard clustering method (by 8%), and the proposed routing algorithm has been proven to be 1.5 times more efficient than the commercial solution software GAMS.
Overall, the proposed methodology has the following attractive features for solving the issues faced due to LIPP in metropolitan contexts: (i) The service considers the trade-off between operator profit and requester cost. (ii) At the expense of increasing the passengers' walking times, the proposed service can effectively reduce the passengers' ride times and the operating costs of the fleet relative to existing door-to-door services. This service is suitable for use in cases of insufficient capacity. The LIPDRS approach can also be used for the scheduling of liquefied natural gas (LNG) buses or pure electric buses. Additionally, it can be used for transport between residential areas and central business districts and for acyclic transportation related to holidays and large events. The use of this mode of transport can help to support green and sustainable metropolitan development.
In this study, only one type of vehicle is considered for LIPDRS planning. However, for a metropolis with diverse demands, this assumption is not appropriate. In future research, we will consider hybrid fleets with different types of vehicles (e.g., with different speeds and capacities) to improve both operator profits and passenger costs. In addition, the whole LIPDRS design process requires some preliminary preparation time for finding the optimal c value, which limits the capability of online implementation. Therefore, we will seek improvements in computational efficiency in future research. For example, some intelligent optimization methods [61][62][63] or cloud computing techniques [64] can be used to obtain solutions efficiently.