Data-Driven Robust Integrated Car-Sharing Scheduling with Heterogeneous Vehicle Considering Uncertain Dispatcher and Relocation Times

1. Introduction

Car-sharing transportation, as a typical example of shared transportation (see Figure 1), has become an important travel mode for urban residents [1]. The global car-sharing market is estimated at USD 11.52 billion in 2025 and is projected to reach USD 28.67 billion by 2030, corresponding to a CAGR of 20.0% over 2025–2030 [2]. Car-sharing offers novel solutions to urban transportation challenges like traffic congestion, taxi availability issues, and parking difficulties [3]. Additionally, it ensures a comparable level of travel comfort to private vehicles [4]. Car-sharing services come in two main types: station-based and free-floating. In densely populated cities, there will be little parking space and high parking costs, which will limit the service capacity of free-floating systems [5]. Within the station-based model, services can be further classified as either one-way or round-way. Consumers are increasingly prioritizing comfort and convenience, causing the traditional round-way model to lose its attractiveness gradually [6].

However, one-way station-based car-sharing services may face a spatial mismatch between vehicle availability and customer demand [7]. (i) How to address this imbalance through the development of effective scheduling plans? In practical operations, uncertainty in both dispatcher movement time and vehicle relocation time may prevent scheduling plans from being executed accurately, thereby affecting the stability of order fulfillment. In the long run, ensuring order fulfillment rates and subsequently improving customer satisfaction are key to establishing a competitive advantage for the car-sharing platform [8]. So, (ii) how can scheduling plans maintain good performance in the face of uncertainty? Additionally, (iii) how can a car-sharing system be designed to meet diverse customer demands and provide personalized services?

To develop effective scheduling plans that address supply-demand imbalances, vehicle relocation becomes a key factor influencing the performance of these plans. Vehicle relocation determines which vehicle should be moved from its supply location to its demand location. This process is relatively complex because, unlike other small-scale shared mobility systems like bicycles, vehicles cannot be moved in bulk on a single carrier [9]. Therefore, the process of relocating vehicles necessitates dispatchers driving them to the demand station. Dispatchers can move around the station where the vehicle is located by folding bicycles, scooters [10], public transport [11] or shuttle vehicles [9]. In order to address supply-demand imbalances, integrated scheduling plans are required to simultaneously coordinate order allocation, vehicle relocation and dispatcher movement.

In practical operations, both dispatcher movement time and vehicle relocation time can be influenced by various subjective and objective factors, such as individual riding/driving proficiency, weather conditions, road congestion, and excessive traffic lights, resulting in variability in execution time [12]. Utilizing deterministic models may not offer effective scheduling plans to meet customer demands. To address uncertainties, researchers often turn to methods such as stochastic programming. However, stochastic programming has its limitations, particularly in determining the data distribution of uncertain parameters. This becomes especially challenging when dealing with large-scale problems and high-dimensional historical data. Robust optimization handles time uncertainty by optimizing decisions against the worst-case realization in a prescribed uncertainty set, while a small set of protection parameters controls the level of conservatism [13]. With the growing availability of operational data, uncertainty sets can be calibrated directly from historical observations within a data-driven robust optimization framework. In this paper, we develop a data-driven interval uncertainty set from historical records of the Taxi & Limousine Commission (TLC) in Queens, New York, by estimating bounds for dispatcher movement time and vehicle relocation time, and use a budget parameter to cap how many time components can simultaneously take their worst-case values under this set. Therefore, data-driven robust optimization produces integrated schedules that are less sensitive to travel-time disruptions, leading to more stable execution and higher order-fulfillment rates in practice.

The growing diversity of travel demands has encouraged car-sharing platforms to offer heterogeneous vehicle services [14]. While most studies assume a homogeneous fleet, real-world fleets are mixed, and ignoring heterogeneity can distort assignment and relocation decisions. In this paper, we model vehicle heterogeneity by categorizing vehicles based on seating capacity, considering three types: 2-seat, 5-seat, and 7-seat vehicles, which serve solo/couple, small-group, and larger family/group trips, respectively. Accordingly, each order must be served by a vehicle whose capacity is no smaller than the requested number of passengers [15]. By determining the appropriate fleet size, the quantity of different vehicle types and the allocation of vehicle tasks, the overall cost can be minimized [16]. In summary, this paper investigates an integrated robust scheduling system for one-way station-based car-sharing with heterogeneous vehicles to meet customer demands. The main contributions of this paper are as follows:

(1) To the best of our knowledge, in the domain of one-way station-based car-sharing, we first integrate the issues of order allocation, vehicle relocation, dispatcher movement, uncertain dispatcher movement time and vehicle relocation time into the problem;

(2) We consider a heterogeneous vehicle fleet, where 2-seat vehicles are for solo/couple, 5-seat for small groups, and 7-seat for larger family/group trips, to better meet diverse customer demand within the integrated scheduling framework;

(3) Based on analysis of historical data from TLC that is similar to this problem, uncertainty sets for dispatcher movement time and vehicle relocation time are defined in interval form by a data-driven method. Specifically, a robust optimization model that can be solved quickly and effectively is proposed to improve order fulfillment rate and the stability of integrated scheduling plans.

The remainder of the paper is organized as follows. Section 2 reviews the related literature. Section 3 describes the one-way station-based car-sharing setting and formulates the integrated scheduling problem. Section 4 presents the deterministic MILP and its data-driven robust counterpart under uncertain dispatcher movement and vehicle relocation time. Section 5 designs a hybrid algorithm for robust models, denoted as ALNS/GA. Section 6 reports numerical experiments, including scenario comparisons and sensitivity analyses. Section 7 concludes the paper and discusses limitations and future research directions.

2. Literature Review

2.1. Vehicle Relocation and Dispatcher Activity

Relocation tasks can be performed either directly by customers [3] or by the service provider [17]. In customer-based relocation, the platform offers incentives to induce users to relocate vehicles toward stations that alleviate distribution imbalances [18]. Zhang et al. [3] implement this idea through price discounts that encourage customers to return vehicles to desired locations. In practice, however, the realized effect of customer-based relocation may vary with users’ willingness to participate and respond to incentives [19]. Consequently, many operators also rely on provider-based relocation, where the platform transfers vehicles from surplus locations to shortage locations to restore balance [17]. Li et al. [20] illustrate that operator-based relocation is frequently embedded in daily operations to mitigate imbalances in vehicle availability. Xu and Wu [21] further refine this operator-based perspective by explicitly modeling how relocation actions reshape vehicle availability over space and time.

Early service-provider studies often model relocation mainly from the vehicle side and focus on deciding where vehicles should be moved, while leaving dispatcher activity implicit [11]. In practice, however, operator-based relocation must be carried out by dispatchers: a dispatcher needs to travel to the pickup station, take the vehicle, and drive it to a shortage station. These steps consume time and dispatcher resources, so a relocation plan that ignores dispatcher activity may be difficult to implement [22]. Yu et al. [23] include dispatcher-based relocation as part of the rebalancing design, indicating that platform relocation still relies on dispatchers to move vehicles in daily operations. Guo et al. [6] then point out a basic requirement: to relocate vehicles, dispatchers must also travel between stations, and this travel affects whether a relocation plan is feasible and how good it is in practice. Bruck et al. [24] state the same idea in a bike-sharing setting and draw a direct implication for modeling: a rebalancing plan should specify dispatcher routes and station-to-station travel, not only the target vehicle distribution. Therefore, vehicle relocation and dispatcher activity should be considered jointly in daily operations.

Some studies integrate vehicle relocation and dispatcher activity, achieving rebalancing of car-sharing system by dispatchers to drive vehicles from surplus stations to shortage stations [6]. Wu et al. [25] consider four types of travel arcs for dispatchers, including waiting, maintenance, movement, and relocation. Yang et al. [26] categorize vehicle travel arcs into three types: parking, user driving, and dispatcher relocating. Tian et al. [27] formulate vehicle relocation and dispatcher rebalancing as a joint operational problem in a one-way free-floating car-sharing system and construct relocation plans by assigning vehicle relocation tasks to dispatchers in a way that respects dispatcher availability and station-to-station travel. Guo et al. [28] integrate vehicle rebalancing with operational relocation effort by jointly optimizing charging schedules and inter-station relocation decisions. From a routing perspective, Gleditsch et al. [29] generate executable operator tours for dynamic rebalancing, linking relocation quantities to feasible routes and the resulting visiting sequence of locations. In recent studies, dispatcher movement between stations is typically realized by small personal devices such as folding bicycles, scooters, public transport or shuttle vehicles. Chang et al. [10] model dispatcher movement by allowing dispatchers to use a folding bicycle or an electric scooter to travel to the next relocation task. Eilertsen et al. [11] allow dispatchers to move between tasks by either a folding bicycle or public transport. Liu et al. [9] model dispatcher movement through dedicated shuttles that transport dispatchers between stations.

2.2. Uncertain Factor

There is relatively limited research on uncertainties in car-sharing systems, with most studies focusing on demand uncertainty and employing stochastic or multi-stage methods. Specifically, Oliveira et al. [30] incorporate a choice-based demand model into operator decisions, providing a structured way to represent demand variability when optimizing relocation actions. Zhang et al. [3] formulate a two-stage stochastic programming model under uncertain demand and derive operator relocation plans to mitigate supply–demand imbalances. Basciftci et al. [14] develop a two-stage stochastic mixed-integer program that integrates multi-period planning decisions with operational relocations under demand uncertainty. In summary, some studies mainly consider demand uncertainty, while less attention has been paid to other factors that may affect the operational efficiency of car-sharing systems, such as dispatcher movement time and relocation travel time.

Robust optimization has been increasingly applied to address uncertainty in scheduling problems [31]. Compared to traditional methods that rely on inputs with random distributions, robust optimization minimizes input requirements and enhances model versatility to tackle challenges like the curse of dimensionality and limited historical data [32]. In shared-mobility rebalancing, Bruck et al. [24] formulate rebalancing as a robust optimization problem and show that robust counterparts can still be solved exactly with tailored algorithms. Chen et al. [33] incorporate distributional robustness into free-floating bike rebalancing, so that the resulting plan is less sensitive to demand perturbations. Zhang et al. [34] push this direction further by calibrating the uncertainty description from historical observations and embedding it in a multi-period robust repositioning model. Beyond demand-side uncertainty, execution-time uncertainty is also crucial in shared mobility, and it has been increasingly incorporated into routing and operational optimization models [35]. Time uncertainty has also been treated explicitly in robust routing models, where uncertain travel times are embedded into the feasibility requirements of time-dependent constraints [36]. In a real-time dispatch-and-relocation setting, Chen et al. [37] handle travel-time uncertainty through a distributionally robust chance-constrained formulation while making vehicle relocation decisions online. Overall, although time uncertainty has been increasingly considered in shared-mobility operations, robust car-sharing scheduling still rarely models uncertainty in vehicle relocation time or in dispatcher movement time, and studies that consider both simultaneously are even rarer.

2.3. Research Gap

Table 1 summarizes the features of previous research efforts. The distinctions of this study from previous research are as follows: Firstly, References [8,21,23,24,25,28,29,32,33] are similar to this study as they also address the dual scheduling of one-way station-based car-sharing system. However, these studies all adopt deterministic models, without considering the uncertainties in real-world scenarios. Secondly, References [3,22,26,34,39] address the uncertainty of customer demands but do not tackle integrated scheduling problem involving both vehicles and dispatchers. To our knowledge, there is currently no research that simultaneously considers the uncertainty in dispatcher movement time and vehicle relocation time, or integrates such uncertainties into the integrated scheduling problem. Thirdly, although the importance of heterogeneous fleets has been studied in the fields of VRP and LRP, it has only received some attention in the domain of car-sharing in recent years. Smet [38] investigates a car-sharing system with heterogeneous fleets, where customers are assigned to alternative vehicles when the requested vehicle type is unavailable. However, this study is limited to a round-way scenario and does not consider vehicle relocation. Based on the above comparison, few studies in one-way station-based car-sharing simultaneously integrate scheduling, a heterogeneous fleet, and uncertainty in both dispatcher movement time and vehicle relocation time.

3. Problem Description

3.1. System Setting

Consider a station-based car-sharing system for one-way rentals where customers provide booking details such as start time, end time, departure station and destination station. We define the start time, end time, departure station and destination station for each demand as $d_c,{\overline{d}}_c,u_c,{\overline{u}}_c$, respectively. The car-sharing platform offers heterogeneous vehicle services with known types and quantities. The location and number of stations are predetermined, each station with limited capacity. Each service-providing vehicle may not necessarily start its next task from the destination station it arrived at after completing the previous one. Failure to meet customer demands incurs penalty costs. This heterogeneity in service leads to significant spatial imbalance between vehicle availability and customer demands. Therefore, dispatchers are required to ride scooters between stations, and relocate vehicles to meet demands and address imbalance issues. Dispatchers begin executing tasks from their initial stations at the start of scheduling and return to those stations by the end of scheduling. We name the problem as car-sharing integrated scheduling with heterogeneous fleet, abbreviated as CISH for later references.

The topology of the one-way station-based car-sharing system can be described using a time-space network, where nodes correspond to stations at different points in time, and arcs denote the interactions involving customers, vehicles, and dispatchers. There are three types of activities in the system: travel, relocation and movement. Travel refers to customers using vehicles through the system. Relocation involves dispatchers driving vehicles between stations. Movement refers to dispatchers riding scooters from one station to the station where relocation tasks need to be executed.

Figure 2 illustrates the integrated scheduling process of the entire car-sharing system. Initially, orders are allocated to vehicles based on customer request information. As imbalances arise between customer demand and vehicle availability, the algorithm integrates customer departure times, pick-up station, drop-off station, the current location of dispatchers and vehicles to calculate potential relocation arcs. Integrated scheduling involves order allocation, vehicle relocation and dispatcher movement. Vehicle relocation routes contain 2-3, 4-5, 7-8, 9-10, 11-6, 13-10 and 14-15. The initial stations with three dispatchers are 1, 6, and 12. The movement routes are 3-4, 8-9, 10-11 and 10-14. The dummy arcs for dispatcher departures and returns from stations are as follows: 1-2, 5-1, 6-7, 12-13, and 15-12. Taking the first dispatcher as an example, he departs from station 1 to supply station 2 and drives the required vehicle to demand station 3. Then, he rides a scooter to supply station 4 and delivers the required vehicles to station 5. Finally, he completes all tasks and returns to his starting point. Each dispatcher is assigned a unique initial station, which serves as the dispatcher’s departure and return location during the scheduling horizon. These initial stations do not undertake any demand-related tasks to maintain the structural clarity of the scheduling model. Figure 3 complements Figure 2 by showing how integrated scheduling decisions translate into executable time sequences for individual dispatchers. Each panel decomposes the execution time into dispatcher movement, vehicle relocation, and the customer-service interval. For task 1, dispatcher 1 performs a direct vehicle relocation from station 4 to station 1, followed by the customer-service segment. For task 2, dispatcher 2 first moves from station 3 to station 4 and then relocates a vehicle from station 4 to station 1; the movement segment, relocation segment, and customer-service segment are displayed separately to reflect their temporal order. By accounting for these execution-time components, the integrated schedule better matches vehicles to bookings, improving order fulfillment and service reliability of the car-sharing system.

3.2. Time Discretization for Data-Driven Uncertainty Sets

To formalize the above operations on the time–space network and ensure temporal consistency in subsequent order allocation, vehicle relocation, and dispatcher movement, we discretize the planning horizon as follows [39]. We consider a scheduling horizon $[T_0,T_{\mathrm{end}}]$ and discretize it into equal-length periods of $\Delta$ minutes. Let $T=\{1,2,\dots ,|T|\}$ denote the set of discrete time indices, where

$$\left|T\right|=⌊{\displaystyle \frac{T_{\mathrm{end}}-T_0}{\Delta }}⌋+1,\tau (t)=T_0+(t-1)\Delta ,\forall t\in T.$$

Here $\tau (t)$ is the real clock time corresponding to index t. The discrete time indices are used to map customer trips, construct candidate scheduling arcs, and update station inventories over the planning horizon. The exact start times of dispatcher tasks remain continuous and are measured in minutes to preserve scheduling accuracy. To avoid constructing arcs that leave the planning horizon, we only allow new customer trips, vehicle relocations and dispatcher movements to depart from periods $t\in T\setminus \{|T|\}$. The last period $\left|T\right|$ is used solely to complete ongoing tasks and record the final inventories at stations.

For each demand $c\in C$, the booking data provide a continuous start time ${\widehat{d}}_c$, a continuous end time ${\widehat{\overline{d}}}_c$, a departure station $u_c$, and a destination station ${\overline{u}}_c$. The continuous times are mapped to discrete period indices by

$$\begin{array}{c}\hfill d_c=1+⌊{\displaystyle \frac{{\widehat{d}}_c-T_0}{\Delta }}⌋,{\overline{d}}_c=1+⌊{\displaystyle \frac{{\widehat{\overline{d}}}_c-T_0}{\Delta }}⌋,\forall c\in C.\end{array}$$

In this way, $d_c$ represents the earliest period in which demand c can start using a vehicle, and ${\overline{d}}_c$ represents the latest period by which the vehicle must be returned. Therefore, the continuous time interval $[{\widehat{d}}_c,{\widehat{\overline{d}}}_c]$ is fully contained in the discrete interval $[\tau (d_c),\tau ({\overline{d}}_c)]$. Hence, each customer travel task is represented by a fixed travel arc $(u_c,d_c,{\overline{u}}_c,{\overline{d}}_c)$ in the time–space network.

Given the nominal relocation time $t_{ij}$ and movement time ${\pi }_{ij}$ between stations i and j, their discrete counterparts in periods are defined as $t_{ij}^{\Delta }=\lceil t_{ij}/\Delta \rceil$ and ${\pi }_{ij}^{\Delta }=\lceil {\pi }_{ij}/\Delta \rceil$, respectively. A time–space arc $a=(i,t,j,k)$ is constructed if departing from station i at period t reaches station j at

$$\begin{array}{cc}\hfill k& =t+t_{ij}^{\Delta }(\mathrm{vehicle}\mathrm{relocation}),\hfill \\ \hfill k& =t+{\pi }_{ij}^{\Delta }(\mathrm{dispatcher}\mathrm{movement}),\hfill \end{array}$$

which guarantees that all arcs in the network are aligned with the discrete time grid.

Note that t and k are indices of coarse time periods, while $t_{ij}$ and ${\pi }_{ij}$ are continuous-time durations (in minutes). Dividing by $\Delta$ converts a duration into the number of periods it spans. The ceiling operator guarantees temporal feasibility on the coarse time axis, i.e., the actual arrival time $\tau (t)+t_{ij}$ (or $\tau (t)+{\pi }_{ij}$) always falls no later than the end of period k.

It is worth noting that the proposed model adopts a two-scale representation of time. Coarse-grained time indices $t\in T$ are used to construct the time–space network and enforce vehicle inventory balance at stations. In contrast, the continuous scheduling variables $t_a^s$, $t_a^{\min }$ and $t_a^{\max }$ are defined in minutes to capture the finer dispatching sequence of relocation tasks within each period. This two-layer time structure allows us to keep the network size tractable while retaining sufficient temporal resolution for dispatcher scheduling.

For the budget uncertainty sets introduced later, we adopt a rolling window with physical length W minutes, corresponding to $H^{\mathrm{win}}=⌈{\displaystyle \frac{W}{\Delta }}⌉$ consecutive periods. For each $t\in T$, the associated window is defined as

$$\begin{array}{c}\hfill {\mathbb{N}}_t^{H^{\mathrm{win}}}=\{r\in T:t\le r\le \min (t+H^{H^{\mathrm{win}}}-1,\left|T\right|)\}.\end{array}$$

Hence, the budget parameter $\Gamma$ controls the maximum number of periods that may simultaneously attain worst-case deviations within each window $N_t^{H^{\mathrm{win}}}$.

4. Model

In this section, we first formulate a deterministic mixed-integer linear programming (MILP) model for the CISH. We then construct data-driven budgeted interval uncertainty sets for vehicle relocation time and dispatcher movement time. Based on these sets, the deterministic model is reformulated into its robust counterpart by robustifying the key temporal-feasibility constraints and the generalized operating-cost objective via standard dual reformulation.

4.1. Deterministic Model of CISH

4.1.1. Notations

Sets

C	Set of customer demands
S	Set of dispatchers
T	Set of discrete time periods, $T=\{1,\dots ,t,\dots ,|T|\}$
I	Set of stations
H	Set of vehicle types
$H^c$	Set of vehicle types suitable for customer demand $c\in C$, $H^c\subseteq H$
A	Set of all potential vehicle-relocation arcs
$A_{\mathrm{start}}^s$	Set of arcs through which dispatcher s departs from the initial station
$A_{\mathrm{end}}^s$	Set of arcs through which dispatcher s returns to the initial station
$A^s$	Set of all relocation arcs associated with dispatcher s, where $A^s=A\cup A_{\mathrm{start}}^s\cup A_{\mathrm{end}}^s$
$A_1^s$	Set of potential relocation arcs and initial departure arcs for dispatcher s, where $A_1^s=A\cup A_{\mathrm{start}}^s$
$A_2^s$	Set of potential relocation arcs and terminal return arcs for dispatcher s, where $A_2^s=A\cup A_{\mathrm{end}}^s$
a	Time–space arc $a=(i,t,j,k)$, representing a route that departs from station i in period t and arrives at station j in period k. Similarly, $a^{\prime }=(i^{\prime },t^{\prime },j^{\prime },k^{\prime })$, $a^{″}=(j,t^{\prime },j^{\prime },k^{\prime })$, and $a^{‴}=(j,k,i,t)$

Indices and parameters

$i,j,i^{\prime },j^{\prime }$	Indices for stations
$t,k,t^{\prime },k^{\prime }$	Indices for time
c	Index for customer demands, $c\in C$
s	Index for dispatchers, $s\in S$
h	Index for vehicle types, $h\in H$
$q_i$	Capacity of station i
$t_{ij}$	Nominal relocation time (minutes) for moving a vehicle from station i to station j
${\pi }_{ij}$	Nominal movement time (minutes) for a dispatcher moving from station i to station j
$t_{ij}^{\Delta }$	Discrete vehicle relocation duration (in periods) from station i to station j, defined as $t_{ij}^{\Delta }=⌈t_{ij}/\Delta ⌉$
${\pi }_{ij}^{\Delta }$	Discrete dispatcher movement duration (in periods) from station i to station j, defined as ${\pi }_{ij}^{\Delta }=⌈{\pi }_{ij}/\Delta ⌉$
$\alpha$	Penalty for abandoning a customer demand
$\beta$	Vehicle relocation cost per unit time
$\gamma$	Cost of vehicle relocation by a dispatcher per unit time
$\delta$	Dispatcher movement cost per unit time
$n_{(i,1)}^h$	Initial number of vehicles of type h at station i
M	A sufficiently large number
$F^h$	Fixed cost of a vehicle of type h
$d_c$	Discrete start-period index of demand c, mapped from ${\widehat{d}}_c$
${\overline{d}}_c$	Discrete end-period index of demand c, mapped from ${\widehat{\overline{d}}}_c$
$u_c$	Departure station of demand c
${\overline{u}}_c$	Destination station of demand c
$\Delta$	Length of each discrete time period (minutes) used to discretize the planning horizon
$\tau (t)$	Clock time corresponding to discrete index t, i.e., $\tau (t)=T_0+(t-1)\Delta$, $\forall t\in T$
W	Physical length (minutes) of the rolling time window in the budgeted uncertainty sets, which induces $H^{\mathrm{win}}=⌈W/\Delta ⌉$ consecutive periods
$n_{(i,t)}^h$	Number of vehicles of type h at station $i\in I$ in period $t\in \{2,\dots ,|T|\}$

Decision variables

$x_c^h$	Binary variable equal to 1 if customer demand c is served by a vehicle of type h, and 0 otherwise
$y_a^{sh}$	Binary variable equal to 1 if dispatcher s relocates a vehicle of type h along relocation arc $a\in A^s$, and 0 otherwise
$z_{a{a}^{\prime }}^{sh{h}^{\prime }}$	Binary variable equal to 1 if dispatcher s relocates a vehicle of type h along arc $a\in A_1^s$ and subsequently relocates a vehicle of type $h^{\prime }$ along arc $a^{\prime }\in A_2^s$, and 0 otherwise
$e_a^h$	Number of vehicles of type h relocated along arc $a\in A$
$t_a^s$	Starting time (minutes) at which dispatcher s begins the relocation task along arc $a\in A^s$

4.1.2. Model Composition

The deterministic model of CISH is formulated as follows.

Objective function: The minimization objective function includes the penalty cost of unmet demands, vehicles and dispatchers relocation cost, dispatchers movement cost, and the fixed cost of vehicles h used to meet demand in operation hours, which are collectively referred to as generalized daily operating costs.

$$\begin{array}{cc}\hfill P& =\min \alpha \sum _{c\in C}\sum _{h\in H^c}(1-x_c^h)+\beta \sum _{h\in H}\sum _{a\in A}{t}_{ij}{e}_a^h+\gamma \sum _{s\in S}\sum _{h\in H}\sum _{a\in A}{t}_{ij}{y}_a^{sh}\hfill \\ \hfill & +\delta (\sum _{s\in S}\sum _{h\in H}\sum _{a\in A_1^s}\sum _{a^{\prime }\in A_2^s}{\pi }_{j{i}^{\prime }}{z}_{a{a}^{\prime }}^{sh{h}^{\prime }}+\sum _{s\in S}\sum _{h\in H}\sum _{a\in A_{start}^s\cup A_{end}^s}{\pi }_{ij}{y}_a^{sh})+\sum _{c\in C}\sum _{h\in H^c}{F}^h{x}_c^h,\hfill \end{array}$$

(1)

where $\sum _{s\in S}\sum _{h\in H}\sum _{a\in A_{start}^s\cup A_{end}^s}{\pi }_{ij}{y}_a^{sh}$ in the dispatchers movement cost is the dummy relocation of dispatchers from and to the depot.

Flow conservation: The available vehicles at station i at time $t+1$ is calculated as follows: it equals the number of available vehicles at time t plus the vehicles that arrived at station i during the interval $[t,t+1)$, minus the vehicles that departed from station i during the interval $[t,t+1)$, plus the vehicles relocated to station i and minus those relocated from station i during the same interval. We assume that no travels, relocations, or movements occur at the final time $\left|T\right|$.

$$\begin{array}{cc}\hfill & n_{(i,t+1)}^h=n_{(i,t)}^h-\sum _{c\in C:u_c=i,d_c=t}{x}_c^h+\sum _{c\in C:{\overline{u}}_c=i,{\overline{d}}_c=t}{x}_c^h-\sum _{\begin{array}{c}j\in I,k\in T:\\ a\in A\end{array}}{e}_a^h+\sum _{\begin{array}{c}j\in I,k\in T:\\ a^{‴}\in A\end{array}}{e}_{a^{‴}}^h,\forall i\in I,\forall h\in H,\forall t\in T\setminus \{|T|\}.\hfill \end{array}$$

(2)

Vehicles and stations: Constraint (3) ensures that, for each vehicle type h, the number of type-h vehicles available at station i at time t, together with the type-h vehicles relocated to station i at time t, is no less than the total number of type-h vehicles that are assigned to customer trips or relocation tasks departing from station i at time t. Constraint (4) ensures that each demand is allocated to at most one type of vehicle. Constraint (5) states that at any given time, the number of vehicles at each station must not exceed the station’s capacity.

$$\begin{array}{cc}\hfill & \sum _{c\in C:u_c=i,d_c=t}{x}_c^h+\sum _{\begin{array}{c}j\in I,k\in T:\\ a\in A\end{array}}{e}_a^h\le n_{(i,t)}^h+\sum _{\begin{array}{c}j\in I,k\in T:\\ a^{‴}\in A\end{array}}{e}_{a^{‴}}^h,\forall i\in I,\forall t\in T,\forall h\in H,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill & \sum _{h\in H^c}{x}_c^h\le 1\forall c\in C,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & n_{(i,t)}^h\le q_i\forall i\in I,\forall t\in T,\forall h\in H.\hfill \end{array}$$

(5)

Scheduling: Constraint (6) states that relocation tasks are performed exclusively by dispatchers, hence the number of times vehicles are relocated equals the number of relocation tasks executed by dispatchers. Constraints (7) and (8) ensure that if dispatcher s does not perform a relocation task on arc a, then any relocation task on arcs connected to this arc, either preceding or succeeding it, cannot be executed. Constraint (9) ensures that for dispatcher s, performing a relocation task on arc $a^{\prime }$ must be preceded by execution of a connected arc a; in other words, each relocation task cannot occur independently and must be linked to a preceding arc. Constraint (10) ensures the flow balance of dispatchers, requiring that for each dispatcher, the number of relocation flows departing from arc a equals the number of relocation flows entering it. Constraint (11) specifies that the start time for relocating vehicles on arc $a^{\prime }$ by a dispatcher should be greater than or equal to the start time for relocating vehicles on arc a, plus the relocation time from i to j, and the time it takes the dispatcher to move from j to $i^{\prime }$. Constraint (12) is similar to constraint (11) but applies to dummy relocation arcs. Since the end station of a dummy relocation arc coincides with the start station of the subsequent relocation task, only the dispatcher’s movement time between arcs is considered, and the vehicle relocation time is excluded. Constraints (13) and (14) ensure that each dispatcher departs from their initial station at the start of scheduling and returns to the initial station at the end of scheduling.

$$\begin{array}{cc}\hfill & e_a^h=\sum _s{y}_a^{sh}\forall a\in A,\forall h\in H,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & z_{a{a}^{\prime }}^{sh{h}^{\prime }}\le y_a^{sh}\forall s\in S,\forall a\in A_1^s,\forall a^{\prime }\in A_2^s,\forall h\in H,\forall h^{\prime }\in H,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & z_{a{a}^{\prime }}^{sh{h}^{\prime }}\le y_{a^{\prime }}^{sh}\forall s\in S,\forall a\in A_1^s,\forall a^{\prime }\in A_2^s,\forall h\in H,\forall h^{\prime }\in H,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & \sum _{a\in A_1^s}{z}_{a{a}^{\prime }}^{sh{h}^{\prime }}=y_{a^{\prime }}^{s{h}^{\prime }}\forall s\in S,\forall a^{\prime }\in A_2^s,\forall h\in H,\forall h^{\prime }\in H,\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & \sum _{a^{\prime }\in A_2^s}{z}_{a{a}^{\prime }}^{sh{h}^{\prime }}=\sum _{a^{\prime }\in A_1^s}{z}_{a^{\prime }a}^{sh{h}^{\prime }}\forall s\in S,\forall a^{\prime }\in A,\forall h\in H,\forall h^{\prime }\in H,\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & t_a^s+t_{ij}+{\pi }_{j{i}^{\prime }}\le t_{a^{\prime }}^s+M\left(1-z_{a{a}^{\prime }}^{sh{h}^{\prime }}\right)\forall s\in S,\forall a\in A,\forall a^{\prime }\in A_2^s,\forall h\in H,\forall h^{\prime }\in H,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & t_a^s+{\pi }_{ij}\le t_{a^{″}}^s+M(1-z_{a^{″}}^{sh{h}^{\prime }})\forall s\in S,\forall a\in A_{start}^s,\forall a^{″}\in A_2^s,\forall h\in H,\forall h^{\prime }\in H,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & \sum _{a\in A_{start}^s}{y}_a^{sh}=1\forall s\in S,\forall h\in H,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & \sum _{a\in A_{end}^s}{y}_a^{sh}=1\forall s\in S,\forall h\in H.\hfill \end{array}$$

(14)

Time window: Constraint (15) is the time window constraint for dispatcher s to perform the relocation task on arc a departing from station i. It ensures that the start time $t_a^s$ must lie within the bounds defined by the earliest and latest feasible start times, $t_a^{min}$ and $t_a^{max}$, respectively. Specifically, $t_a^{min}$ denotes the earliest possible start time (in minutes) of the relocation task on arc a, determined by the dynamic scheduling requirements, whereas $t_a^{max}$ denotes the latest possible start time.

$$\begin{array}{c}\hfill t_a^{min}-M(1-y_a^{sh})\le t_a^s\le t_a^{max}+M(1-y_a^{sh})\forall s\in S,\forall a\in A^s,\forall h\in H.\end{array}$$

(15)

Variables: Constraints (16)-(18) restrict that $x_c^h,y_a^{sh}$ and $z_{a{a}^{\prime }}^{sh{h}^{\prime }}$ are binary variables. Constraints (19)-(22) define that $e_a^h,n_{(i,t)}^h,t_a^s,t_a^{min}$ and $t_a^{max}$ are integer variables.

$$\begin{array}{cc}\hfill & x_c^h\in \{0,1\}\forall c\in C,\forall h\in H^c,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & y_a^{sh}\in \{0,1\}\forall s\in S,\forall a\in A_s,\forall h\in H,\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & z_{a{a}^{\prime }}^{sh{h}^{\prime }}\in \{0,1\}\forall s\in S,\forall a\in A_1^s,\forall a^{\prime }\in A_2^s,\forall h\in H,\forall h^{\prime }\in H,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & e_a^h\in \mathbb{N}\forall a\in A,\forall h\in H,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & n_{(i,t)}^h\in \mathbb{N}\forall i\in I,\forall t\in \{2,\dots ,|T|\},\forall h\in H,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & t_a^s\in \mathbb{N}\forall s\in S,\forall a\in A_s,\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & t_a^{min},t_a^{max}\in \mathbb{N}\forall a\in {\cup }_{s\in S}{A}^s.\hfill \end{array}$$

(22)

4.2. Robust Counterpart for CISH

4.2.1. Uncertainty Set Construction

In real-world operations, the relocation time of dispatchers is subject to various sources of uncertainty, including traffic congestion, dynamic route changes, unexpected detours, and heterogeneous driving proficiency. As a result, the actual relocation time may deviate from its nominal estimate and exhibit time-varying fluctuations. Following the data-driven budgeted-uncertainty-set framework of literature [40], we construct a budget-based uncertainty set over a rolling time window $N_t^{H^{win}}$, which consists of $H^{win}$ consecutive periods from t to $(t+H^{win}-1)$.

Specifically, for each window starting at t, the interval bounds and the associated budget parameters are calibrated from historical observations and can be updated in a rolling-horizon manner, while the robust optimization problem solved at t treats these calibrated quantities as given inputs. The uncertainty set is formulated as follows:

$$\begin{array}{c}\hfill {\mathcal{T}}_t^{ij}=\left\{{({\tilde{t}}_{t,r}^{ij})}_{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}\in {\mathbb{R}}^{|{\mathbb{N}}_t^{H^{\mathrm{win}}}|}|\begin{array}{c}{\displaystyle {\tilde{t}}_{t,r}^{ij}\in [l_{t,r}^t,u_{t,r}^t],{\displaystyle \sum _{k=t}^r}\left|\frac{{\tilde{t}}_{t,k}^{ij}-{\overline{m}}_{t,k}^t}{{\widehat{m}}_{t,k}^t}\right|\le {\Gamma }_{t,r}^t,\forall r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}\hfill \end{array}\right\}.\end{array}$$

(23)

(1) For each period $r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}$, $[l_{t,r}^t,u_{t,r}^t]$ is the interval support of the uncertain relocation time ${\tilde{t}}_{t,r}^{ij}$.

(2) ${\overline{m}}_{t,r}^t=(u_{t,r}^t+l_{t,r}^t)/2$ and ${\widehat{m}}_{t,r}^t=(u_{t,r}^t-l_{t,r}^t)/2$ are the mean and the half-interval length of ${\tilde{t}}_{t,r}^{ij}$, respectively.

(3) ${\Gamma }_{t,r}^t\in [0,r-t+1]$ is used to control the conservatism of the uncertainty set, representing the cumulative budget of uncertainty.

In real-world scenarios, dispatchers move by scooter, which may be influenced by multifarious factors such as road congestion, weather conditions, excessive traffic lights, unclear endpoint locations and other subjective and objective factors, leading to changes in movement time. Consequently, the actual movement time may deviate from the expected value, exhibiting a certain degree of volatility. We model the dispatchers’ movement time ${\pi }_{ij}$ as a budget-type uncertainty over the same rolling window ${\mathbb{N}}_t^{H^{\mathrm{win}}}$:

$$\begin{array}{c}\hfill {\Pi }_t^{ij}=\left\{{({\tilde{\pi }}_{t,r}^{ij})}_{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}\in {\mathbb{R}}^{|{\mathbb{N}}_t^{H^{\mathrm{win}}}|}|\begin{array}{c}{\displaystyle {\tilde{\pi }}_{t,r}^{ij}\in [l_{t,r}^{\pi },u_{t,r}^{\pi }],{\displaystyle \sum _{k=t}^r}\left|\frac{{\tilde{\pi }}_{t,k}^{ij}-{\overline{m}}_{t,k}^{\pi }}{{\widehat{m}}_{t,k}^{\pi }}\right|\le {\Gamma }_{t,r}^{\pi },\forall r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}\hfill \end{array}\right\}.\end{array}$$

(24)

(1) For each period $r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}$, $[l_{t,r}^{\pi },u_{t,r}^{\pi }]$ is the interval support of the uncertain movement time ${\tilde{\pi }}_{t,r}^{ij}$.

(2) ${\overline{m}}_{t,r}^{\pi }=(u_{t,r}^{\pi }+l_{t,r}^{\pi })/2$ and ${\widehat{m}}_{t,r}^{\pi }=(u_{t,r}^{\pi }-l_{t,r}^{\pi })/2$ are the mean and the half-interval length of ${\tilde{\pi }}_{t,r}^{ij}$, respectively.

(3) ${\Gamma }_{t,r}^{\pi }\in [0,r-t+1]$ is used to control the conservatism of the uncertainty set, representing the cumulative budget of uncertainty.

4.2.2. Robust Counterparts of Constraints

Proposition 1. Constraint (11) can be reformulated as constraint (25).

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{t}_a^s+{\overline{m}}_{t,n}^t+{\overline{m}}_{t,n}^{\pi }+\sum _{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}\left({\Gamma }_{t,r}^t{\mu }_{t,n,r}^t+{\gamma }_{t,n,r}^t+{\Gamma }_{t,r}^{\pi }{\mu }_{t,n,r}^{\pi }+{\gamma }_{t,n,r}^{\pi }\right)\le t_{a^{\prime }}^s+M(1-z_{a{a}^{\prime }}^{sh{h}^{\prime }}),\hfill \\ \forall a\in A,\forall a^{\prime }\in A_2^s,\forall s\in S,\forall h\in H,\forall n\in {\mathbb{N}}_t^{H^{\mathrm{win}}},\hfill \\ {\mu }_{t,n,r}^t+{\gamma }_{t,n,r}^t\ge {\widehat{m}}_{t,r}^t\forall r\in {\mathbb{N}}_t^{n-t+1},\forall n\in {\mathbb{N}}_t^{H^{\mathrm{win}}},\hfill \\ {\mu }_{t,n,r}^{\pi }+{\gamma }_{t,n,r}^{\pi }\ge {\widehat{m}}_{t,r}^{\pi }\forall r\in {\mathbb{N}}_t^{n-t+1},\forall n\in {\mathbb{N}}_t^{H^{\mathrm{win}}},\hfill \\ {\mu }_{t,n,r}^t,{\gamma }_{t,n,r}^t\ge 0\forall r\in {\mathbb{N}}_t^{n-t+1},\forall n\in {\mathbb{N}}_t^{H^{\mathrm{win}}},\hfill \\ {\mu }_{t,n,r}^{\pi },{\gamma }_{t,n,r}^{\pi }\ge 0\forall r\in {\mathbb{N}}_t^{n-t+1},\forall n\in {\mathbb{N}}_t^{H^{\mathrm{win}}}.\hfill \end{array}\right.\end{array}$$

(25)

Proof. Given the current period t, for each end period $n\in {\mathbb{N}}_t^{H^{\mathrm{win}}}$ we consider the sub-horizon ${\mathbb{N}}_t^{n-t+1}=\{t,\dots ,n\}$. Within the uncertainty set ${\mathcal{T}}_t^{ij}$ in Eq. (23), we define the scaled deviations for each $r\in {\mathbb{N}}_t^{n-t+1}$ as

$${\tilde{z}}_{t,r}^t=({\tilde{t}}_{t,r}^{ij}-{\overline{m}}_{t,r}^t)/{\widehat{m}}_{t,r}^t,{\tilde{z}}_{t,r}^t\in [-1,1].$$

Similarly, for the movement-time uncertainty set ${\Pi }_t^{ij}$ in Eq. (24), the scaled deviations are

$${\tilde{z}}_{t,r}^{\pi }=({\tilde{\pi }}_{t,r}^{ij}-{\overline{m}}_{t,r}^{\pi })/{\widehat{m}}_{t,r}^{\pi },{\tilde{z}}_{t,r}^{\pi }\in [-1,1].$$

Restricted to the sub-horizon ${\mathbb{N}}_t^{n-t+1}$, the sets ${\mathcal{T}}_t^{ij}$ and ${\Pi }_t^{ij}$ can be equivalently written as

$$\begin{array}{cc}\hfill & {\mathcal{T}}_{t,n}^{ij}=\left\{{\tilde{t}}_{t,n}^{ij}\in {\mathbb{R}}^{n-t+1}|\begin{array}{c}{\tilde{t}}_{t,r}^{ij}={\overline{m}}_{t,r}^t+{\widehat{m}}_{t,r}^t{\tilde{z}}_{t,r}^t,|{\tilde{z}}_{t,r}^t|\le 1,{\displaystyle \sum _{k=t}^r}\left|{\tilde{z}}_{t,k}^t\right|\le {\Gamma }_{t,r}^t,\forall r\in {\mathbb{N}}_t^{n-t+1}\hfill \end{array}\right\},\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & {\Pi }_{t,n}^{ij}=\left\{{\tilde{\pi }}_{t,n}^{ij}\in {\mathbb{R}}^{n-t+1}|\begin{array}{c}{\tilde{\pi }}_{t,r}^{ij}={\overline{m}}_{t,r}^{\pi }+{\widehat{m}}_{t,r}^{\pi }{\tilde{z}}_{t,r}^{\pi },|{\tilde{z}}_{t,r}^{\pi }|\le 1,{\displaystyle \sum _{k=t}^r}\left|{\tilde{z}}_{t,k}^{\pi }\right|\le {\Gamma }_{t,r}^{\pi },\forall r\in {\mathbb{N}}_t^{n-t+1}.\hfill \end{array}\right\}\hfill \end{array}$$

(27)

Therefore, constraint (11), $t_a^s+t_{ij}+{\pi }_{j{i}^{\prime }}\le t_{a^{\prime }}^s+M\left(1-z_{a{a}^{\prime }}^{sh{h}^{\prime }}\right),\forall s\in S,\forall a\in A,\forall a^{\prime }\in A_2^s,\forall h\in H,\forall h^{\prime }\in H$, is equivalent to the following inequality

$$\begin{array}{cc}\hfill & t_a^s+{\overline{m}}_{t,r}^t+{\widehat{m}}_{t,r}^t{\tilde{z}}_{t,r}^t+{\overline{m}}_{t,r}^{\pi }+{\widehat{m}}_{t,r}^{\pi }{\tilde{z}}_{t,r}^{\pi }\le t_{a^{\prime }}^s+M(1-z_{a{a}^{\prime }}^{sh{h}^{\prime }}).\hfill \end{array}$$

(28)

To obtain the robust counterpart, we maximize the deviation terms on the left-hand side. Introducing the auxiliary quantities

$$\begin{array}{cc}\hfill {\mathrm{AUX}}^t=& \underset{\{{\tilde{z}}_{t,r}^t\}}{\max }\left\{{\displaystyle \sum _{r=t}^n}{\widehat{m}}_{t,r}^t{\tilde{z}}_{t,r}^t|{\displaystyle \sum _{k=t}^r}\left|{\tilde{z}}_{t,k}^t\right|\le {\Gamma }_{t,r}^t,\left|{\tilde{z}}_{t,r}^t\right|\le 1,\forall r\in {\mathbb{N}}_t^{n-t+1}\right\},\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill {\mathrm{AUX}}^{\pi }=& \underset{\{{\tilde{z}}_{t,r}^{\pi }\}}{\max }\left\{{\displaystyle \sum _{r=t}^n}{\widehat{m}}_{t,r}^{\pi }{\tilde{z}}_{t,r}^{\pi }|{\displaystyle \sum _{k=t}^r}\left|{\tilde{z}}_{t,k}^{\pi }\right|\le {\Gamma }_{t,r}^{\pi },\left|{\tilde{z}}_{t,r}^{\pi }\right|\le 1,\forall r\in {\mathbb{N}}_t^{n-t+1}\right\}.\hfill \end{array}$$

(30)

Inequality Eq. (28) can be conservatively rewritten as

$$\begin{array}{c}\hfill t_a^s+\left({\overline{m}}_{t,n}^t+{\mathrm{AUX}}^t\right)+\left({\overline{m}}_{t,n}^{\pi }+{\mathrm{AUX}}^{\pi }\right)\le t_{a^{\prime }}^s+M\left(1-z_{a{a}^{\prime }}^{sh{h}^{\prime }}\right).\end{array}$$

(31)

Dual reformulation. Following the dualization approach of Bertsimas and Sim [41], we introduce dual variables ${\mu }_{t,n,r}^t,{\gamma }_{t,n,r}^t\ge 0$ for the relocation time, and ${\mu }_{t,n,r}^{\pi },{\gamma }_{t,n,r}^{\pi }\ge 0$ for the movement time. The dual forms of Eq. (29) and Eq. (30) are:

$$\begin{array}{cc}\hfill {\mathrm{AUX}}^t& =\underset{{\mu }_{t,n,r}^t,{\gamma }_{t,n,r}^t\ge 0}{\min }\left\{\sum _{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}{\Gamma }_{t,r}^t{\mu }_{t,n,r}^t+{\gamma }_{t,n,r}^t|{\mu }_{t,n,r}^t+{\gamma }_{t,n,r}^t\ge {\widehat{m}}_{t,r}^t,\forall r\in {\mathbb{N}}_t^{n-t+1}\right\},\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\mathrm{AUX}}^{\pi }& =\underset{{\mu }_{t,n,r}^{\pi },{\gamma }_{t,n,r}^{\pi }\ge 0}{\min }\left\{\sum _{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}{\Gamma }_{t,r}^{\pi }{\mu }_{t,n,r}^{\pi }+{\gamma }_{t,n,r}^{\pi }|{\mu }_{t,n,r}^{\pi }+{\gamma }_{t,n,r}^{\pi }\ge {\widehat{m}}_{t,r}^{\pi },\forall r\in {\mathbb{N}}_t^{n-t+1}\right\}.\hfill \end{array}$$

(33)

Substituting these dual forms into constraint (31) yields the linear robust counterpart of constraint (11):

$$\begin{array}{c}\hfill t_a^s+{\overline{m}}_{t,n}^t+{\overline{m}}_{t,n}^{\pi }+\sum _{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}\left({\Gamma }_{t,r}^t{\mu }_{t,n,r}^t+{\gamma }_{t,n,r}^t+{\Gamma }_{t,r}^{\pi }{\mu }_{t,n,r}^{\pi }+{\gamma }_{t,n,r}^{\pi }\right)\le t_{a^{\prime }}^s+M(1-z_{a{a}^{\prime }}^{sh{h}^{\prime }}).\end{array}$$

(34)

The associated feasibility constraints of the dual variables are:

$$\begin{array}{c}\hfill {\mu }_{t,n,r}^t+{\gamma }_{t,n,r}^t\ge {\widehat{m}}_{t,r}^t,{\mu }_{t,n,r}^{\pi }+{\gamma }_{t,n,r}^{\pi }\ge {\widehat{m}}_{t,r}^{\pi },{\mu }_{t,n,r}^t,{\gamma }_{t,n,r}^t,{\mu }_{t,n,r}^{\pi },{\gamma }_{t,n,r}^{\pi }\ge 0,\forall r\in {\mathbb{N}}_t^{n-t+1}.\end{array}$$

(35)

Thus, constraint (11) is replaced by its robust counterpart Eq. (34), ensuring that the temporal scheduling between consecutive relocation arcs remains feasible under any realization of uncertain relocation and movement times within their respective budget uncertainty sets.

Proposition 2. The robust dual linearized form of constraint (12) under the budget-based uncertainty set of dispatcher movement time ${\Pi }_t^{ij}$ is given as constraint (36).

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle t_a^s+{\overline{m}}_{t,n}^{\pi }+\sum _{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}\left({\Gamma }_{t,r}^{\pi }{\mu }_{t,n,r}^{\pi }+{\gamma }_{t,n,r}^{\pi }\right)\le t_{a^{″}}^s+M(1-z_{a{a}^{″}}^{sh{h}^{\prime }}),}\hfill \\ \forall a\in A_{start}^s,\forall a^{″}\in A_2^s,\forall s\in S,\forall h\in H,\forall n\in {\mathbb{N}}_t^{H^{\mathrm{win}}},\hfill \\ {\mu }_{t,n,r}^{\pi }+{\gamma }_{t,n,r}^{\pi }\ge {\widehat{m}}_{t,r}^{\pi },\forall r\in {\mathbb{N}}_t^{n-t+1},\forall n\in {\mathbb{N}}_t^{H^{\mathrm{win}}},\hfill \\ {\mu }_{t,n,r}^{\pi },{\gamma }_{t,n,r}^{\pi }\ge 0\forall r\in {\mathbb{N}}_t^{n-t+1},\forall n\in {\mathbb{N}}_t^{H^{\mathrm{win}}}.\hfill \end{array}\right.\end{array}$$

(36)

Proof. Given period t and end period $n\in {\mathbb{N}}_t^{H^{\mathrm{win}}}$, we again restrict ${\Pi }_t^{ij}$ in Eq. (24) to the sub-horizon ${\mathbb{N}}_t^{n-t+1}$ and define the scaled deviations as

$${\tilde{z}}_{t,r}^{\pi }=({\tilde{\pi }}_{t,r}^{ij}-{\overline{m}}_{t,r}^{\pi })/{\widehat{m}}_{t,r}^{\pi },{\tilde{z}}_{t,r}^{\pi }\in [-1,1],\forall r\in {\mathbb{N}}_t^{n-t+1}.$$

The set ${\Pi }_t^{ij}$ is represented equivalently by the following inequality.

$$\begin{array}{cc}\hfill & {\Pi }_t^{ij}=\left\{{\tilde{\pi }}_{t,n}^{ij}\in {\mathbb{R}}^{n-t+1}|\begin{array}{c}{\tilde{\pi }}_{t,r}^{ij}={\overline{m}}_{t,r}^{\pi }+{\widehat{m}}_{t,r}^{\pi }{\tilde{z}}_{t,r}^{\pi },|{\tilde{z}}_{t,r}^{\pi }|\le 1,{\displaystyle \sum _{k=t}^r}\left|{\tilde{z}}_{t,k}^{\pi }\right|\le {\Gamma }_{t,r}^{\pi },\forall r\in {\mathbb{N}}_t^{n-t+1}\hfill \end{array}\right\}.\hfill \end{array}$$

(37)

Therefore, constraint (12), $t_a^s+{\pi }_{ij}\le t_{a^{″}}^s+M\left(1-z_{a{a}^{″}}^{sh{h}^{\prime }}\right),\forall s\in S,\forall a\in A_{start}^s,\forall a^{″}\in A_2^s,\forall h\in H$ is equivalent to the following inequality

$$\begin{array}{c}\hfill t_a^s+{\overline{m}}_{t,r}^{\pi }+{\widehat{m}}_{t,r}^{\pi }{\tilde{z}}_{t,r}^{\pi }\le t_{a^{″}}^s+M(1-z_{a{a}^{″}}^{sh{h}^{\prime }}).\end{array}$$

(38)

By maximizing the deviation term on the left-hand side, we obtain the auxiliary maximization problem:

$$\begin{array}{c}\hfill {\mathrm{AUX}}_{start}^{\pi }=\underset{\{{\tilde{z}}_{t,r}^{\pi }\}}{\max }\left\{{\displaystyle \sum _{r=t}^n}{\widehat{m}}_{t,r}^{\pi }{\tilde{z}}_{t,r}^{\pi }|{\displaystyle \sum _{k=t}^r}|{\tilde{z}}_{t,k}^{\pi }|\le {\Gamma }_{t,r}^{\pi },\left|{\tilde{z}}_{t,r}^{\pi }\right|\le 1,\forall r\in {\mathbb{N}}_t^{n-t+1}\right\}.\end{array}$$

(39)

Dual reformulation. Following the dualization approach of Bertsimas and Sim [41], by introducing the dual variables ${\mu }_{t,n,r}^{\pi },{\gamma }_{t,n,r}^{\pi }\ge 0$, the dual form of the auxiliary problem is obtained as:

$$\begin{array}{c}\hfill {\mathrm{AUX}}_{start}^{\pi }=\underset{{\mu }_{t,n,r}^{\pi },{\gamma }_{t,n,r}^{\pi }\ge 0}{\min }\left\{\sum _{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}{\Gamma }_{t,r}^{\pi }{\mu }_{t,n,r}^{\pi }+{\gamma }_{t,n,r}^{\pi }|{\mu }_{t,n,r}^{\pi }+{\gamma }_{t,n,r}^{\pi }\ge {\widehat{m}}_{t,r}^{\pi },\forall r\in {\mathbb{N}}_t^{n-t+1}\right\}.\end{array}$$

(40)

Substituting the dual expression of ${\mathrm{AUX}}_{start}^{\pi }$ into constraint (38) yields the linearized robust dual constraint (41), that is:

$$\begin{array}{c}\hfill t_a^s+{\overline{m}}_{t,n}^{\pi }+\sum _{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}\left({\Gamma }_{t,r}^{\pi }{\mu }_{t,n,r}^{\pi }+{\gamma }_{t,n,r}^{\pi }\right)\le t_{a^{″}}^s+M(1-z_{a{a}^{″}}^{sh{h}^{\prime }}),\end{array}$$

(41)

together with the dual variable constraints

$$\begin{array}{c}\hfill {\mu }_{t,n,r}^{\pi }+{\gamma }_{t,n,r}^{\pi }\ge {\widehat{m}}_{t,r}^{\pi },{\mu }_{t,n,r}^{\pi },{\gamma }_{t,n,r}^{\pi }\ge 0,\forall r\in {\mathbb{N}}_t^{n-t+1}.\end{array}$$

(42)

This completes the robust dualization of constraint (12), ensuring temporal feasibility between successive virtual arcs under any realization within the uncertainty set ${\Pi }_t^{ij}$.

4.2.3. Robust Counterpart of the Objective Function

In addition to the temporal feasibility constraints, the relocation time $t_{ij}$ and dispatcher movement time ${\pi }_{ij}$ also appear in the generalized daily operating cost Eq. (1) as coefficients of e, y, and z. Hence, we further robustify the objective under the same rolling-window budget uncertainty sets ${\mathcal{T}}_t^{ij}$ and ${\Pi }_t^{ij}$.

Proposition 3.   Under the rolling-window budget uncertainty sets ${\mathcal{T}}_t^{ij}$ and ${\Pi }_t^{ij}$, the robust counterpart of the objective function Eq. (1) admits the following equivalent linear reformulation:

$$\begin{array}{cc}\hfill P=\min & \alpha \sum _{c\in C}\sum _{h\in H^c}(1-x_c^h)+\sum _{c\in C}\sum _{h\in H^c}{F}^h{x}_c^h+\sum _{i\in I}\sum _{j\in I}\sum _{t\in T}\sum _{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}\left({\overline{m}}_{t,r}^t{\Theta }_{t,r}^{t,ij}+{\overline{m}}_{t,r}^{\pi }{\Theta }_{t,r}^{\pi ,ij}\right)\hfill \\ \hfill & +\sum _{i\in I}\sum _{j\in I}\sum _{t\in T}\sum _{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}\left({\Gamma }_{t,r}^t{\mu }_{t,r}^{t,ij,\mathrm{obj}}+{\gamma }_{t,r}^{t,ij,\mathrm{obj}}+{\Gamma }_{t,r}^{\pi }{\mu }_{t,r}^{\pi ,ij,\mathrm{obj}}+{\gamma }_{t,r}^{\pi ,ij,\mathrm{obj}}\right),\hfill \end{array}$$

(43)

subject to the exposure-defining constraints

$$\begin{array}{cc}\hfill {\Theta }_{t,r}^{t,ij}& =\beta \sum _{h\in H}\sum _{a\in A_{t,r}^{ij}}{e}_a^h+\gamma \sum _{s\in S}\sum _{h\in H}\sum _{a\in A_{t,r}^{ij}}{y}_a^{sh},\forall i,j\in I,\forall t\in T,\forall r\in {\mathbb{N}}_t^{H^{\mathrm{win}}},\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill {\Theta }_{t,r}^{\pi ,ij}& =\delta \sum _{s\in S}\sum _{h\in H}\sum _{a\in (A_{start}^s\cup A_{end}^s)\cap A_{t,r}^{ij}}{y}_a^{sh}+\delta \sum _{s\in S}\sum _{h\in H}\sum _{h^{\prime }\in H}\sum _{a\in A_{1,t,r}^{s,i}}\sum _{a^{\prime }\in A_2^{s,j}}{z}_{a{a}^{\prime }}^{sh{h}^{\prime }},\forall i,j\in I,\forall t\in T,\forall r\in {\mathbb{N}}_t^{H^{\mathrm{win}}},\hfill \end{array}$$

(45)

where for an arc $a=(i(a),t(a),j(a),k(a))\in A$, we define $A_{t,r}^{ij}:=\{a\in A:i(a)=i,t(a)=t,j(a)=j,k(a)=r\}$, $A_{1,t,r}^{s,i}:=\{a\in A_1^s:t(a)=t,j(a)=i,k(a)=r\}$, and $A_2^{s,j}:=\{a^{\prime }\in A_2^s:i(a^{\prime })=j\}$.

The dual variables satisfy the following feasibility constraints:

$$\begin{array}{cc}\hfill {\mu }_{t,r}^{t,ij,\mathrm{obj}}+{\gamma }_{t,r}^{t,ij,\mathrm{obj}}& \ge {\widehat{m}}_{t,r}^t{\Theta }_{t,r}^{t,ij},\forall i,j\in I,\forall t\in T,\forall r\in {\mathbb{N}}_t^{H^{\mathrm{win}}},\hfill \\ \hfill {\mu }_{t,r}^{\pi ,ij,\mathrm{obj}}+{\gamma }_{t,r}^{\pi ,ij,\mathrm{obj}}& \ge {\widehat{m}}_{t,r}^{\pi }{\Theta }_{t,r}^{\pi ,ij},\forall i,j\in I,\forall t\in T,\forall r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}.\hfill \end{array}$$

(46)

All dual variables satisfy the following nonnegativity constraints:

$$\begin{array}{c}\hfill {\mu }_{t,r}^{t,ij,\mathrm{obj}},{\gamma }_{t,r}^{t,ij,\mathrm{obj}},{\mu }_{t,r}^{\pi ,ij,\mathrm{obj}},{\gamma }_{t,r}^{\pi ,ij,\mathrm{obj}}\ge 0,\forall i,j\in I,\forall t\in T,\forall r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}.\end{array}$$

(47)

Proof.   

$$\begin{array}{cc}\hfill P=\min \underset{\begin{array}{c}({\tilde{t}}_{t,r}^{ij})\in {\mathcal{T}}_t^{ij},\forall i,j,t\\ ({\tilde{\pi }}_{t,r}^{ij})\in {\Pi }_t^{ij},\forall i,j,t\end{array}}{\max }& \alpha \sum _{c\in C}\sum _{h\in H^c}(1-x_c^h)+\sum _{c\in C}\sum _{h\in H^c}{F}^h{x}_c^h+\sum _{i\in I}\sum _{j\in I}\sum _{t\in T}\sum _{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}\left({\tilde{t}}_{t,r}^{ij}{\Theta }_{t,r}^{t,ij}+{\tilde{\pi }}_{t,r}^{ij}{\Theta }_{t,r}^{\pi ,ij}\right).\hfill \end{array}$$

(48)

Following the scaled-deviation representation in Proposition 1–2, for each fixed $(i,j)$ and window starting at t, we write ${\tilde{t}}_{t,r}^{ij}={\overline{m}}_{t,r}^t+{\widehat{m}}_{t,r}^t{\tilde{z}}_{t,r}^{t,ij}$ and ${\tilde{\pi }}_{t,r}^{ij}={\overline{m}}_{t,r}^{\pi }+{\widehat{m}}_{t,r}^{\pi }{\tilde{z}}_{t,r}^{\pi ,ij}$, where $|{\tilde{z}}_{t,r}^{t,ij}|\le 1$, ${\sum }_{q=t}^r\left|{\tilde{z}}_{t,q}^{t,ij}\right|\le {\Gamma }_{t,r}^t$, and $|{\tilde{z}}_{t,r}^{\pi ,ij}|\le 1$, ${\sum }_{q=t}^r\left|{\tilde{z}}_{t,q}^{\pi ,ij}\right|\le {\Gamma }_{t,r}^{\pi }$ for all $r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}$.

Here, ${\Theta }_{t,r}^{t,ij}$ and ${\Theta }_{t,r}^{\pi ,ij}$ collect the decision-dependent exposure terms and are defined in constraints (44) and (45). Substituting the midpoint–deviation representations into Eq. (48) separates the objective into a nominal part and a deviation part. For each fixed $(i,j,t)$, the deviation part leads to the following two auxiliary maximization problems:

$$\begin{array}{cc}\hfill {\mathrm{AUX}}_{\mathrm{obj}}^{t,ij}=& \underset{\{{\tilde{z}}_{t,r}^{t,ij}\}}{\max }\left\{\sum _{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}{\widehat{m}}_{t,r}^t{\Theta }_{t,r}^{t,ij}{\tilde{z}}_{t,r}^{t,ij}|{\displaystyle \sum _{q=t}^r}\left|{\tilde{z}}_{t,q}^{t,ij}\right|\le {\Gamma }_{t,r}^t,\left|{\tilde{z}}_{t,r}^{t,ij}\right|\le 1,\forall r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}\right\},\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill {\mathrm{AUX}}_{\mathrm{obj}}^{\pi ,ij}=& \underset{\{{\tilde{z}}_{t,r}^{\pi ,ij}\}}{\max }\left\{\sum _{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}{\widehat{m}}_{t,r}^{\pi }{\Theta }_{t,r}^{\pi ,ij}{\tilde{z}}_{t,r}^{\pi ,ij}|{\displaystyle \sum _{q=t}^r}\left|{\tilde{z}}_{t,q}^{\pi ,ij}\right|\le {\Gamma }_{t,r}^{\pi },\left|{\tilde{z}}_{t,r}^{\pi ,ij}\right|\le 1,\forall r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}\right\}.\hfill \end{array}$$

(50)

Dual reformulation. Following the same dualization approach of Bertsimas and Sim [41] used in Proposition 1–2, introduce dual variables ${\mu }_{t,r}^{t,ij,\mathrm{obj}},{\gamma }_{t,r}^{t,ij,\mathrm{obj}}\ge 0$ for ${\mathrm{AUX}}_{\mathrm{obj}}^{t,ij}$, and ${\mu }_{t,r}^{\pi ,ij,\mathrm{obj}},{\gamma }_{t,r}^{\pi ,ij,\mathrm{obj}}\ge 0$ for ${\mathrm{AUX}}_{\mathrm{obj}}^{\pi ,ij}$. Their dual-feasibility conditions are exactly constraint (46), and the corresponding dual objective values yield the additional linear terms

$$\begin{array}{c}\hfill \sum _{r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}}\left({\Gamma }_{t,r}^t{\mu }_{t,r}^{t,ij,\mathrm{obj}}+{\gamma }_{t,r}^{t,ij,\mathrm{obj}}+{\Gamma }_{t,r}^{\pi }{\mu }_{t,r}^{\pi ,ij,\mathrm{obj}}+{\gamma }_{t,r}^{\pi ,ij,\mathrm{obj}}\right),\end{array}$$

(51)

in the outer minimization. Summing over all $(i,j,t)$ and combining with the nominal part completes the derivation, leading to the linear robust objective Eq. (43) with constraints (44)–(46).

Combining Proposition 1–3, the robust dual form of the car-sharing integrated scheduling problem is a mixed-integer linear programming model obtained by minimizing Eq. (43) subject to constraints (2)–(10), constraints (13)–(22), the robust scheduling constraints (25) and (36), and the robust-objective constraints (44)–(46).

5. Solution Method

The car-sharing integrated robust scheduling problem yields a large and highly constrained search space, and its difficulty grows rapidly with the fleet size. Moreover, execution-time uncertainty in vehicle relocation and dispatcher movement further complicates both feasibility and objective evaluation. To obtain high-quality solutions within acceptable computation times, we designed a hybrid algorithm for robust models, denoted as ALNS/GA. Unlike a standard ALNS or a purely deterministic GA, the proposed ALNS/GA incorporates additional procedures to account for the calibrated uncertainty sets: robustness is enforced during chromosome decoding and repair via robust time-feasibility checks, and solution quality is measured by the robust objective value. Algorithm 1 presents the resulting workflow, where GA provides global exploration and ALNS operators strengthen local improvement under the robust counterpart.

Algorithm 1 ALNS/GA

Require:  Original data; Population set P; Maximum iteration $T_{\max }$; Crossover rate $P_c$; Mutation rate $P_m$; Iteration counter $n_{ite}$; Removal limit n; Removal criterion $h_{\mathrm{rem}}$; Insertion criterion $h_{\mathrm{ins}}$ Ensure:  Final solution S 1: $\{{\mathcal{T}}_t^{ij},{\Pi }_t^{ij},{\Gamma }_t^t,{\Gamma }_t^{\pi }\}\leftarrow \mathrm{C}\mathrm{ALIBRATE}\mathrm{U}\mathrm{NCERTAINTY}\mathrm{S}\mathrm{ETS}(\mathrm{historical}\mathrm{data})$ 2: Generate initial population of set P 3: for each chromosome g in the population do 4:     $S\leftarrow \mathrm{R}\mathrm{OBUST}\mathrm{D}\mathrm{ECODE}\mathrm{R}\mathrm{EPAIR}(g,{\mathcal{T}}_t^{ij},{\Pi }_t^{ij})$ 5:     $F_{\mathrm{rob}}(S)\leftarrow \mathrm{R}\mathrm{OBUST}\mathrm{E}\mathrm{VALUATE}(S,{\mathcal{T}}_t^{ij},{\Pi }_t^{ij})$ 6:     $\mathrm{fit}(g)\leftarrow 1/(F_{\mathrm{rob}}(S)+\epsilon )$ 7: end for 8: while $n_{ite}<T_{\max }$do 9:     Select chromosomes using roulette wheel selection based on $\mathrm{fit}(·)$ 10:     if $P_c>$ a random number then 11:         Perform best cost route crossover and generate offspring chromosomes 12:     end if 13:     if $P_m>$ a random number then 14:         Perform single fragment mutation and generate offspring chromosomes 15:     end if 16:     for each newly generated chromosome g do 17:         $S\leftarrow \mathrm{R}\mathrm{OBUST}\mathrm{D}\mathrm{ECODE}\mathrm{R}\mathrm{EPAIR}(g,{\mathcal{T}}_t^{ij},{\Pi }_t^{ij})$ 18:         $F_{\mathrm{rob}}(S)\leftarrow \mathrm{R}\mathrm{OBUST}\mathrm{E}\mathrm{VALUATE}(S,{\mathcal{T}}_t^{ij},{\Pi }_t^{ij})$ 19:         $\mathrm{fit}(g)\leftarrow 1/(F_{\mathrm{rob}}(S)+\epsilon )$ 20:     end for 21:     Obtain a new population using elite strategy 22:     for each chromosome g in the best $20\%$ of the new population do 23:         $S\leftarrow \mathrm{R}\mathrm{OBUST}\mathrm{D}\mathrm{ECODE}\mathrm{R}\mathrm{EPAIR}(g,{\mathcal{T}}_t^{ij},{\Pi }_t^{ij})$ 24:         $(S,A)\leftarrow \mathrm{R}\mathrm{EMOVAL}(S,n,h_{\mathrm{rem}})$ 25:         $S\leftarrow \mathrm{I}\mathrm{NSERTION}(S,A,h_{\mathrm{ins}})$ 26:         $F_{\mathrm{rob}}(S)\leftarrow \mathrm{R}\mathrm{OBUST}\mathrm{E}\mathrm{VALUATE}(S,{\mathcal{T}}_t^{ij},{\Pi }_t^{ij})$ 27:         Re-encode S and update $\mathrm{fit}(g)\leftarrow 1/(F_{\mathrm{rob}}(S)+\epsilon )$ 28:     end for 29:     $n_{ite}\leftarrow n_{ite}+1$ 30: end while 31: Output the final solution S

5.1. Robustness Integration in ALNS/GA

We integrate robustness into our ALNS/GA algorithm at three levels. First, we calibrate budgeted interval uncertainty sets for vehicle relocation time and dispatcher movement time using a rolling window of historical observations, and treat the calibrated sets as fixed inputs when solving the robust scheduling instance at time index t. Second, robustness is embedded in solution construction: during chromosome decoding we enforce robust time-feasibility for both dispatcher movements and vehicle relocations under the calibrated uncertainty sets, and repair infeasible tasks by truncating relocation quantities and postponing start times to the earliest robust-feasible period. Third, the search is guided by the robust counterpart: each decoded schedule is evaluated by the worst-case cost $F_{\mathrm{rob}}(·)$ over the calibrated uncertainty sets, with a penalty applied to robust-infeasible schedules.

Motivated by the rolling data-driven calibration procedure of literature [40], we update the uncertainty-set parameters using the most recent historical observations. Specifically, the rolling window allows the interval supports and uncertainty budgets to adapt to time-varying traffic conditions, rather than remaining fixed over the entire planning horizon. At each time index t, we extract historical relocation times and dispatcher movement times within a window of length $H^{\mathrm{win}}$ and, for each OD pair $(i,j)$ and period r in the window, estimate the interval supports $[l_{t,r}^t,u_{t,r}^t]$ and $[l_{t,r}^{\pi },u_{t,r}^{\pi }]$ for ${\tilde{t}}_{t,r}^{ij}$ and ${\tilde{\pi }}_{t,r}^{ij}$, respectively. We then determine the budget parameters ${\Gamma }_t^t$ and ${\Gamma }_t^{\pi }$ via historical backtesting. Specifically, we compute cumulative normalized deviations over the window and choose the budgets to cover these deviations. The resulting sets $\mathcal{T}{t}^{ij}$ and $\Pi t^{ij}$ are fixed throughout the search and used in both robust decoding/repair and the worst-case evaluation $F_{\mathrm{rob}}(·)$. Algorithm 2 summarizes this rolling-window calibration procedure.

Algorithm 2 Data-driven calibration of budgeted interval uncertainty sets

Require:  Historical relocation times and movement times; window length $H^{\mathrm{win}}$; time index t Ensure:  Interval supports $\{l_{t,r}^t,u_{t,r}^t,l_{t,r}^{\pi },u_{t,r}^{\pi }\}$ and budgets ${\Gamma }_t^t,{\Gamma }_t^{\pi }$ 1: Extract observations in the rolling window $r\in {\mathbb{N}}_t^{H^{\mathrm{win}}}$ 2: Estimate interval supports $[l_{t,r}^t,u_{t,r}^t]$ and $[l_{t,r}^{\pi },u_{t,r}^{\pi }]$ for each $(i,j),r$ 3: Compute cumulative normalized deviations over the window implied by historical realizations 4: Set budgets ${\Gamma }_t^t$ and ${\Gamma }_t^{\pi }$ to cover the deviations (data-driven backtesting) 5: Construct ${\mathcal{T}}_t^{ij}$ and ${\Pi }_t^{ij}$ and return them as algorithm inputs

5.2. Decoding Method and Initial Population

According to the characteristics of the integrated scheduling problem, we encode each individual as a chromosome $g=(G,G_T,G_N,G_S,G_E)$, where G represents the distance between stations; $G_T$ denotes the relocation start times; $G_N$ denotes the number of vehicles to be scheduled; $G_S$ and $G_E$ denote the origin and destination stations of a relocation task, respectively.

Decoding and feasibility handling: A chromosome g induces K relocation tasks ${(G_S^k,G_E^k,G_T^k,G_N^k)}_{k=1}^K$. We sort the tasks in nondecreasing order of $G_T^k$ and decode them into a schedule S by stepping through the planning horizon, updating station inventories and enforcing station-capacity and robust dispatcher-feasibility constraints. If a task is infeasible, we repair it by truncating $G_N^k$ to the available vehicles at $G_S^k$ and, when needed, postponing $G_T^k$ to the earliest robust feasible period. Tasks that remain infeasible are removed; chromosomes that still yield an infeasible schedule are discarded and regenerated.

Initial population: Under the premise of satisfying the coding, we randomly generate n individuals to establish the initial population, which is denoted as $P=\{p_1,p_2,...,p_n\}$. Each chromosome is decoded using the above procedure; infeasible chromosomes are repaired when possible, otherwise discarded and regenerated, so that the initial population consists of feasible schedules.

Fitness function: The fitness function reflects the quality of individuals, where a larger value indicates a better individual. For each decoded schedule $S_i$, we evaluate its robust objective value $F_{\mathrm{rob}}(S_i)$. $F_{\mathrm{rob}}(S_i)$ is defined as the worst-case total cost over the uncertainty sets ${\mathcal{T}}_t^{ij}$ and ${\Pi }_t^{ij}$ under the robust counterpart. If $S_i$ is not robust feasible after decoding, a sufficiently large penalty is added to $F_{\mathrm{rob}}(S_i)$ to discourage infeasible individuals. Since the problem is formulated as a minimization model, the fitness is defined as the reciprocal of the robust objective value:

$$fit(i)={\displaystyle \frac{1}{F_{\mathrm{rob}}(S_i)+\epsilon }},$$

where $\epsilon >0$ is a small constant to avoid division by zero.

5.3. Operations of GA

The genetic algorithm seeks the optimal solution through operations such as selection, crossover, and mutation. The improved genetic operator utilized in the proposed ALNS/GA algorithm is outlined below.

Roulette wheel selection: This mechanism employs chromosome fitness values to assess chromosome quality. Higher fitness values indicate better chromosomes and a higher likelihood of inheritance in the next generation. However, this approach often leads to a high number of iterations with chromosomes of high fitness, limiting the algorithm’s search range. To address this limitation, the proposed ALNS/GA algorithm integrates roulette selection with an elite strategy to preserve outstanding individuals in the population and form a new one.

Initially, the fitness values of the current chromosomes $\mathrm{fit}(i)$ are obtained and summed as $\mathrm{sumfit}={\sum }_{i=1}^P\mathrm{fit}(i)$, which yields the selection probability $p(i)=\mathrm{fit}(i)/\mathrm{sumfit}$. The cumulative probability is then computed as $p_s(i)={\sum }_{j=1}^{i}p(j)$ for $i=1,\dots ,P$. A random real number u is generated within the range $[0,1]$, and the first chromosome is selected such that $p_s(i)\ge u$; otherwise, the chromosome i satisfying $p_s(i-1)<u\le p_s(i)$ is selected.

Best cost route crossover: The crossover operator helps to pass on good chromosomes to offspring and determines the global search capability of the algorithm [42]. The basic idea of the best cost route crossover operator is to seek a lower-cost solution based on the minimum-cost principle of insertion. The costing formula [43] for inserting customer c between nodes i and j is $cost(i,c,j)=d_{ic}+d_{cj}-d_{ij}$. After the crossover operation, in order to increase the search range of the algorithm and avoid becoming trapped in a local optimum, 80% of the population was subjected to mutation operations, and the remaining 20% population was subjected to worst removal and basic greedy insertion operations.

Single fragment mutation: Let the mutation probability of the population be $p_v$, and then generate a random number p in the interval (0,1), if $p>p_v$, the parent chromosome does not undergo mutations and does not produce new individuals. Otherwise, the selected chromosome fragment is mutated. A random integer $r(1\le r\le n)$ is generated, indicating that the fragment of the chromosome r is mutated. The fragment of the chromosome r is re-encoded to produce a new chromosome by replacing the original chromosome fragment.

Elite strategy: If exceptional individuals are lost during the iterative process, the algorithm may fail to achieve the global optimal solution. Therefore, following the operations of selection, crossover, and mutation, the elite strategy is employed to preserve outstanding individuals within the population. This strategy ensures the continual presence of the best solution throughout the iterative process by directly transferring the top-performing chromosomes to the next generation. Moreover, the quantity of elite individuals is determined by multiplying the population size by the elite rate. Upon incorporating elite individuals into a new population, an equivalent number of inferior individuals are removed to maintain a consistent population size.

5.4. Removal and Insertion Heuristics

According to the characteristics of the problem, adaptive large neighborhood search operators are used to enlarge the search space. The removal and insertion operators are inspired by the worst removal and basic greedy insertion of Ropke and Pisinger [44]. Different from a purely deterministic setting, the proposed operators are guided by the robust objective value $F_{\mathrm{rob}}(·)$ and the robust feasibility check. In particular, removal prioritizes tasks that contribute most to $F_{\mathrm{rob}}$, whereas insertion selects the task and insertion position that yields the smallest marginal increase in $F_{\mathrm{rob}}$ among robust feasible candidates.

Construct a set A containing all unserved tasks based on model constraints. During each iteration, the algorithm selects a task from A and inserts it into the current solution. An insertion is accepted only if the resulting schedule remains robust feasible under the uncertainty sets; otherwise, alternative positions are tested, and a new route is generated if needed. This process iterates until all tasks in A have been serviced.

Worst removal: This operator removes several genes from their current positions and stores them in a list. To embed robustness into the destroy phase, we evaluate the impact of removing a task on the robust objective. For a task k in the current solution S, define its robust marginal contribution as

$$\Delta F_{\mathrm{rob}}^{\mathrm{rem}}(k)=F_{\mathrm{rob}}(S)-F_{\mathrm{rob}}(S\setminus k).$$

Two removal heuristics are considered as follows.

Maximum Robust Cost Removal: Remove the task with the largest $\Delta F_{\mathrm{rob}}^{\mathrm{rem}}(k)$ from the current solution.

Maximum Robust Risk Removal: Remove the task that is most critical to robust feasibility, measured by the smallest robust time slack along the relocation sequence after decoding.

Basic greedy insertion: This operator reinserts removed tasks so that the search continues in an expanded neighborhood. Candidate insertions are evaluated using the robust objective and are restricted to positions that remain robust feasible after decoding. For inserting a task k at a candidate position p in the current solution S, we compute the robust marginal increase

$$\Delta F_{\mathrm{rob}}^{\mathrm{ins}}(k,p)=F_{\mathrm{rob}}(S\oplus (k,p))-F_{\mathrm{rob}}(S),$$

and only consider positions that remain robust feasible after decoding. Two insertion heuristics are used as follows.

Minimum Robust Cost Insertion: Insert the task and position that yield the minimum $\Delta F_{\mathrm{rob}}^{\mathrm{ins}}(k,p)$ among all robust feasible candidates.

Slack-oriented Insertion: Among robust feasible candidates, prefer the insertion that maximizes the minimum robust time slack after decoding.

6. Numerical Experiments and Case Study

This section evaluates the proposed integrated robust scheduling framework and the ALNS/GA solution method using a case study in Queens, New York. Specifically, we conduct two main sets of experiments to quantify the benefits of integration and robustness and to derive actionable insights for operational decision-making. First, we compare results under three scenarios. By comparing Scenario 1 and Scenario 2, we quantify the benefit of integrated scheduling in terms of total cost, order fulfillment rate, and generalized cost per met demand. By comparing Scenario 2 and Scenario 3, we assess the value of robustness, focusing on solution stability across repeated runs as well as cost and fulfillment performance. Second, we perform sensitivity analyses on the dispatcher scheduling window, fleet size, and heterogeneous vehicle-type availability. Detailed settings are given in the following subsections.

All numerical experiments were conducted on a LAPTOP-HAB1JSMT computer with a 1.80 GHz Intel(R) Core (TM) i5-8265U processor, 8 GB RAM, running Windows 10 64 bits. The proposed method is implemented using Python 3.3. For each instance, we ran the ALNS/GA 10 times.

6.1. Multi-Scenario Analysis

In the following case study, we use the real operation data from TLC in Queens, New York in June 2023 to test the proposed ALNS/GA. Example data can be found on https://www.nyc.gov/site/tlc/about/aggregated-reports.page. For ease of computation, 25 stations are marked, as shown in Figure 4. Here are the parameter configurations: (i) 25 stations; (ii) 5 dispatchers; (iii) customer demands of 100, 200, 300, 400 and 500; (iv) the number of vehicles for different demands are 75, 125, 175, 200 and 250; (v) the speed of vehicle is assumed as 30 km/h and the speed of scooter is assumed as 20 km/h; (vi) the unit penalty cost for unmet one demand is set to 4, the unit relocation cost for both vehicles and dispatchers is 0.3, while the unit movement cost is 0.1.

From June 1st to June 30th, 30 cases were selected for comparative analysis across three scenarios based on five types of customer demands. The comparison results of generalized daily operational costs are shown in Table 2. Bold entries indicate the best result for the same case and demand scale. The three scenarios are as follows:

Scenario 1: Conventional car-sharing system without integrated scheduling and robust optimization. The system processes requests sequentially in chronological order, prioritizing local assignment by matching an available vehicle at the request’s origin depot. If no vehicle is available at the origin, the system searches for available vehicles at neighboring depots and performs a relocation to serve the request; if no feasible vehicle can be found in the neighboring depots, the request is rejected and a penalty cost is incurred;

Scenario 2: Car-sharing integrated scheduling system without robust optimization;

Scenario 3: Car-sharing integrated robust scheduling system.

6.1.1. Comparison of Scenario 1 and Scenario 2

This section evaluates system performance under two settings: the conventional non-integrated scheduling system and the deterministic integrated scheduling system. As shown in Table 2, across all 30 daily cases from June 1 to June 30, Scenario 2 yields lower generalized daily operational costs than Scenario 1 under all five demand scales, with an average cost reduction of approximately 10.69%. A day-by-day comparison further indicates that, except for a few isolated cases with marginal reversals, Scenario 2 achieves cost savings in the vast majority of instances. These results suggest that integrating vehicle assignment with inter-station dispatch decisions can effectively mitigate redundant relocations and unmet-demand penalties induced by myopic matching in the non-integrated system, thereby reducing overall operating costs.

Figure 5 presents six representative days and compares Scenarios 1 and 2 in terms of generalized cost and order fulfillment rate under different demand levels. Overall, as demand increases from 100 to 500, the generalized cost rises in both scenarios. However, Scenario 2 yields lower costs than Scenario 1 at most demand levels. At the same time, the fulfillment rate in Scenario 2 remains consistently high (around 95% or above), whereas Scenario 1 shows a clearer decline and larger fluctuations as demand grows. In several cases under medium-to-high demand, the fulfillment rate in Scenario 1 drops to roughly 75%–85%, which widens the performance gap between the two scenarios.

Figure 5 further reports the generalized cost per met demand. Scenario 2 is generally lower and more stable across the selected days, while Scenario 1 exhibits more pronounced variability. It is also worth noting that, for June 1 with 500 demands, Scenario 1 shows a lower cost per met demand than Scenario 2. This pattern occurs together with a much lower fulfillment rate in Scenario 1, suggesting that this derived metric can decrease mechanically when fewer orders are served. In summary, Figure 5 provides visual evidence that deterministic integrated scheduling typically improves service performance while reducing generalized costs, and it delivers more stable operations under medium and high demand.

6.1.2. Comparison of Scenario 2 and Scenario 3

Based on the best results in Table 2 for 30 daily cases and five demand scales, Scenario 3 achieves a lower total cost than Scenario 2 in most case–scale combinations, and the advantage is more pronounced under medium and high demand. In contrast, under low demand, the cost difference between the two scenarios is usually small, and Scenario 2 is slightly better on some days. This suggests that when supply–demand pressure is limited, the benefit of adopting a robust strategy may be modest.

Figure 6 further reports the outcomes of ten independent runs for Scenarios 2 and 3 under five demand scales on six representative days. In Figure 6, the “blue dashed line” represents the results of ten optimizations using the robust optimization model, while the “purple solid line” represents the results of ten optimizations using the deterministic model. For the same day and demand level, the ten runs fluctuate only slightly around their respective cost levels. In most demand scales, the robust model tends to produce lower and smoother trajectories, indicating smaller dispersion and more stable outputs. Meanwhile, for a few demand levels or runs, the two curves are very close and may even cross. Overall, combining Table 2 and Figure 6 indicates that Scenario 3 is more likely to deliver a lower best cost than Scenario 2 and shows better stability across repeated runs, especially when demand is high.

Through multi-scenario experiments using real TLC operational data from Queens, both integrated scheduling and robust optimization are shown to improve system performance. Compared with the conventional non-integrated system in Scenario 1, the deterministic integrated scheduling system in Scenario 2 performs better across five demand scales and 30 daily cases: the average total cost decreases by 10.69% and the order fulfillment rate increases by 16.48%. This indicates that jointly optimizing vehicle assignment and inter-station dispatch decisions can effectively reduce unmet demand, thereby achieving both cost savings and service improvements. Building on this setting, the integrated robust scheduling system in Scenario 3 further reduces the average total cost by 1.96% relative to Scenario 2 and increases the fulfillment rate by an additional 0.99%. It also exhibits more stable fluctuations across repeated runs, suggesting that robust optimization enhances the reliability and adaptability of decisions when dispatcher movement times and vehicle relocation times are uncertain. Overall, integrated scheduling substantially improves efficiency and service quality, while robust integrated scheduling further strengthens operational stability under uncertainty, providing more reliable support for car-sharing operations.

6.2. Start and End Times of Scheduling

We examine the influence of varying scheduling start and end times on generalized daily operational costs and order fulfillment rate, using 400 customer demands and maintaining parameters consistent with those used in scenario analysis experiments. Figure 7 displays the customer demands time series over 30 days from June 1st to June 30th, utilizing data collected from all 25 stations. The peak periods are observed to be approximately at 6:00-7:00, 15:00-16:00, and 18:00. There are different patterns in both customer demands and vehicle availability across four distinct time stages:

Off-peak period (00:00-05:00): Typically around midnight, are characterized by minimal customer demands, with vehicles mostly remaining idle.

Morning rush period (05:00-10:00): Customer demand exhibits a steady growth trend, indicating high demand for vehicle utilization in this period. This observation aligns with customers’ daily activities, as most individuals typically drive to work in the morning.

Semi-peak period (10:00-16:00): Semi-peak hours in the middle of the day. During this period, customer demand displays a stable state pattern. During this phase, it can be considered another busy period characterized by constant relocation and movement of dispatchers.

Evening-peak period (16:00-24:00): There is a notable surge in customer demands around 18:00, followed by a gradual decline over time. This increase can be attributed to people’s activities after work, such as returning home, dining, and engaging in various social activities.

The aggregate count of vehicles picked up and dropped off at various stations during various time periods from June 1st to June 30th, totaling 30 days, is illustrated in Figure 8. It is evident that stations 9 and 11 have the highest customer demands. Therefore, we allocate dispatchers to stations 2 and 7, which are nearby and have lower demands.

Even though the division of the day into four stages doesn’t exactly match the traditional peak and off-peak traffic categorizations, there are still noticeable resemblances. This suggests that the patterns of pick-ups and drop-offs in car-sharing services closely mirror the typical travel and commuting behaviors seen with ordinary modes of transportation. Based on the analysis of customer demands throughout a day, the scheduling start times for the experiment are set at 8:00, 9:00, 10:00, and 11:00, while the scheduling end times are set at 20:00, 21:00, 22:00, and 23:00 in the evening.

To evaluate how the chosen dispatcher working window influences system performance, we vary the scheduling start time and end time and report the corresponding results in Figure 9. The figure indicates that both the generalized cost and the order fulfillment rate are sensitive to the scheduling window. Compared with the without scheduling setting, the with scheduling strategy generally maintains a higher fulfillment rate and achieves a lower generalized cost across all tested time settings. Moreover, under with scheduling, extending the end time tends to further improve overall performance. In contrast, the impact of changing the start time depends more on demand intensity and the temporal demand profile, and the resulting trends may therefore be non-monotonic or vary across instances.

Following these observations, we include additional results for two extra cases in the appendix, as shown in Figure A1 and Figure A2, to examine the robustness of the findings. Overall, these two cases exhibit patterns similar to those in Figure 9. Scheduling substantially improves the fulfillment rate and reduces the generalized cost relative to the no-scheduling baseline. In addition, under the large-demand setting with 500 requests, delaying the start time leads to a more pronounced deterioration in fulfillment performance, whereas extending the end time still tends to deliver better or more stable outcomes.

6.3. Sensitivity Analysis on Fleet Size

We conduct a fleet-size sensitivity analysis for the June 5 case with demand scale 500, where the fleet size is set to 210, 230, 250, 270, and 290 vehicles. Figure 10a illustrates the effects of fleet size on generalized total cost, satisfaction rate, and marginal satisfaction. The formula for marginal satisfaction concerning fleet size is given by Eq. (52). Here, $r(\xi )$ and $r(\xi +\Delta \xi )$ represent the satisfaction rates under fleet sizes $\xi$ and $\xi +\Delta \xi$, with $\Delta \xi$ set to 20 in this study. Figure 10b reports the impacts of fleet size on operating cost and the number of relocation tasks. Note that operating cost excludes the fixed cost of vehicles as well as the penalty cost incurred from unmet demand.

$$\begin{array}{cc}\hfill & MS(\xi )={\displaystyle \frac{r(\xi +\Delta \xi )-r(\xi )}{\Delta \xi }}\hfill \end{array}$$

(52)

Overall, as fleet size increases, the generalized total cost decreases and the satisfaction rate increases. Meanwhile, the marginal improvement in satisfaction is larger at smaller fleet sizes and then gradually tapers off, indicating diminishing returns from further expansion. A notable exception is that the satisfaction rate drops when the fleet size increases from 210 to 230. As shown in Figure 10b, the number of relocation tasks is much higher at a fleet size of 210, suggesting that under tight supply the system may rely on more aggressive relocation and dispatch actions to maintain a relatively high service level; this also leads to higher operating cost and higher generalized total cost. When the fleet size increases to 230, relocation intensity decreases, which may temporarily weaken the correction of station imbalances and result in a short-term decline in satisfaction.

Another important observation is that when the fleet size increases from 270 to 290, the generalized total cost continues to decrease, but the operating cost increases. This implies that the additional reduction in total cost is mainly driven by fewer unmet-demand penalties, while achieving the extra service improvement requires more relocation and dispatch effort. From an operational perspective, moderate fleet expansion is generally beneficial under high demand. However, when the fleet becomes large, expansion decisions should balance service targets against operating costs, so as to avoid introducing excessive operational burdens for limited marginal gains.

6.4. Sensitivity to the Available Vehicle Types

The data from June 1st to June 30th shows the distribution of customer numbers in orders as follows: 87.5% of orders are for 1-2 passengers, 10.75% are for groups of 3-5 people, and 1.75% are for groups of 6-7 people. Consequently, we establish three different types of vehicles to meet customer demands: two-seater economy cars, five-seater comfortable cars, and seven-seater family cars. It is worth noting that the fixed costs and relocation fees for different vehicles are both incremental.

As shown in Figure 11, restricting the available vehicle types can substantially increase the generalized cost. Under all demand scales, using only five-seater vehicles consistently results in the highest generalized cost. Introducing two-seater vehicles to form a mixed fleet leads to a clear cost reduction, and further adding seven-seater vehicles achieves the lowest generalized cost across all demand levels. Table 3 further clarifies the cost composition. Within the generalized cost, the fixed cost $C_f$ accounts for the largest share, followed by the relocation cost $C_{re}$, while the penalty cost $C_p$ is relatively small. Compared with VT1, VT2 and VT3 reduce $C_f$ noticeably by incorporating smaller-capacity vehicles. In addition, under medium and high demand, the richer vehicle mix provides better matching for heterogeneous requests, which helps reduce penalties and ineffective relocations caused by supply–demand mismatches. It is also worth noting that, under low demand, the cost difference between VT2 and VT3 is small. As demand increases, the cost advantage of VT3 becomes more pronounced, yielding a lower generalized cost than VT2 and VT1. Overall, vehicle-type configuration has a significant impact on both cost and service performance. However, more vehicle types are not always better. A more appropriate strategy is to select a vehicle mix that aligns with the demand scale and the demand structure, so as to achieve a better balance between cost and service.

7. Conclusions

7.1. Theoretical Implications

This study contributes to the car-sharing scheduling literature in three ways. Firstly, it formalizes an integrated scheduling problem for one-way station-based car-sharing by jointly considering order allocation, vehicle relocation, and dispatcher movement, while allowing a heterogeneous fleet. Secondly, it extends integrated scheduling to a data-driven robust setting by calibrating uncertainty sets from historical observations and accounting for time variability in both dispatcher movement and vehicle relocation, thereby reducing sensitivity to execution-time deviations. Thirdly, it develops a hybrid ALNS/GA solution method with robustness embedded in the search and evaluation, so the integrated robust model remains tractable on realistic instances.

7.2. Practical Implications

The empirical results provide several operational takeaways for station-based car-sharing operators. Integrated scheduling can materially improve both cost and service outcomes compared with non-integrated operations: across the tested cases, it reduces average total cost by about 10.69% and increases order fulfillment by about 16.48%. Building on this, the integrated robust system further reduces cost by 1.96% and improves fulfillment by an additional 0.99%, while producing more stable outcomes across repeated runs. Operators can therefore use integrated scheduling as the default approach, and introduce the robust variant when the system scales up or time-uncertainty risk increases to improve operational reliability. In addition, the sensitivity analyses suggest that fleet size and fleet composition should be adjusted with demand pressure: adding vehicles can reduce generalized cost up to a point, but excessive fleet size may cause station congestion and higher operating burdens; a mixed fleet can be beneficial under sufficiently diverse and large demand, but may increase cost under low demand.

7.3. Limitations and Future Research

This study has several limitations that open clear directions for future work. Firstly, the uncertainty treatment focuses on dispatcher movement time and vehicle relocation time; future work could also incorporate uncertainty in demand volume and cancellations within a unified robust framework. Secondly, the analysis is based on a specific operational setting and dataset; validating the approach across different cities, seasons, and station layouts would strengthen external validity and help calibrate uncertainty sets in practice. Thirdly, the current framework abstracts from several operational features that matter in many platforms, such as dynamic pricing, real-time re-optimization, and electric-vehicle charging constraints. Extending the model to include EV charging and energy feasibility is a natural next step.