Mitigating the “black holes”: Periodic repair and maintenance problem of shared bikes

1. Introduction

The sharing economy demonstrates its value when prominent industry players embrace data science to enhance efficiency and reduce costs, particularly in the transportation sector (Zhong et al., 2024). A notable example is the free-floating bike-sharing industry, which has become integral to urban transportation. However, aggressive market expansions driven by venture capital have led to the so-called “Tragedy of the Commons,” resulting in significant metal consumption that undermines environmental goals (Kaspi et al., 2016). Furthermore, a substantial quantity of malfunctioning bikes not only amplifies risks for riders (Wang and Szeto, 2018) but also diminishes the available street space (Yin et al., 2019). Therefore, careful repair is crucial to ensure sustained operation in the bike-sharing industry, accounting for 22% of the overall operational costs (North American Bikeshare & Scootershare Association, 2022).

Through careful analysis of the available data, we discover a compelling computational phenomenon indicating that bikes are more likely to be broken in specific restricted areas. This issue is particularly prevalent in some urban areas with inadequate bike maintenance. For instance, malfunctioning bikes cannot be repaired promptly in some remote regions, or solar-powered locks cannot be recharged indoors. Consistent with our findings, certain neighborhoods exhibit disproportionately high failure rates, highlighting the need for geographically targeted maintenance. According to New York City Comptroller’s Office (2023), “Citi Bike users would encounter twice as many disabled or broken bikes at a typical station in the Bronx than in other boroughs [ ⋯ ] This points to a need for more geographically focused repair and maintenance operations.”

We therefore conceptualize a physical landscape analogous to black holes in cosmology, where certain regions on the bike-sharing map exhibit morbidly high failure rates. This conceptualization is inspired by a broader line of interdisciplinary research that borrows physical analogies to study complex spatial–temporal systems (e.g., Cai et al., 2025; Li et al., 2021) from an optimization perspective. From the user’s perspective, these regions pose significant service risks, as repeated failed pickup attempts may lead to frustration and eventual disengagement. From the operator’s standpoint, these regions are associated with severe asset underutilization. Recognizing this “black hole” phenomenon is thus critical for developing responsive and cost-effective operational strategies. Moreover, our data analysis reveals that failure rates are highly uncertain and shaped by exogenous covariates such as ride histories; for instance, short trips often signal higher failure likelihoods (Mobi, 2024). In light of these findings, we focus on optimizing repair decisions under the black hole phenomenon while accounting for such uncertainty in failure rates within our modeling framework.

Our study is motivated by a collaboration with a leading free-floating bike-sharing operator in China facing frequent bike malfunctions that severely degrade service quality, especially in high-failure-rate regions. Immediate repairs for every malfunction are infeasible due to practical constraints such as cost limitations and geographical distances. For instance, Mobi (2019) remarks that “[b]etween performing regular maintenance checks, we spend a lot of time commuting between stations.” Thus, determining the repair period—a key research direction highlighted by Shui and Szeto (2020)—is a critical decision that balances service quality and operational costs. Overly short repair periods drive up costs, whereas excessively long ones undermine service levels. Beyond balancing service levels and operational costs within individual regions, coordinating repair schedules across regions is equally important. Aligning repair periods for nearby regions or those with similar failure patterns enables route consolidation and reduces redundant travel, making spatiotemporal coordination a key lever for improving cost efficiency in geographically dispersed systems.

In this article, we study the adaptive periodic repair and maintenance problem, accounting for uncertainties in failure rates and related covariate information. The study focuses on a bike-sharing system maintained by a firm that operates across multiple regions, making strategic-level decisions on repair routes and service region assignments. Subsequently, operational-level decisions on adaptive repair periods and preventive maintenance are made to help mitigate the operational impact of the black hole phenomenon as captured through heterogeneous failure rates when stochastic covariate information becomes available. To address the lack of knowledge about the joint distribution of failure rates and covariates, we develop an operational data analytics (ODA) framework that boosts oracle solutions, such as those from distributionally robust optimization (DRO), to improve data-driven decision quality under limited distributional information.

The main results and contributions of this article are summarized as follows.

Conceptually, we conceptualize a physical landscape similar to black holes in cosmology to represent bike-sharing regions with exceptionally high failure rates. The phenomenon, referred to as a black hole in this study, plays a critical role in shared bike repair. In our framework, the black hole phenomenon is captured by regional failure-rate heterogeneity and the resulting penalty costs. To mitigate the impact of high-failure-rate regions and enhance operational efficiency, we propose two special strategies: the dedicated repair period and preventive maintenance intervention. The effectiveness of these two strategies is analytically illustrated within a two-region system.

Methodologically, departing from existing studies that primarily address static repair models, we draw upon the literature on adaptive operations management to introduce a periodic repair and maintenance model for shared bikes. We develop the ODA framework to exploit problem structure and improve decision quality under limited data for the periodic repair problem. To support this approach, nonlinearities in penalty costs are addressed using a linear decision rule (LDR), which is shown to yield optimal solutions in this setting.

Practically, our case study involving a real-world application in the free-floating bike-sharing system yields three interesting insights, utilizing an extensive industry-scale dataset. First, empirical validation demonstrates that the ODA approach generally improves solution quality, even when the oracle data-integration model underperforms. Second, we quantify regional disparities in penalty costs (or service levels) using the well-known Gini coefficient, and the results indicate that a dedicated repair period tends to mitigate the black hole phenomenon more effectively than preventive maintenance. Third, incorporating black hole information not only improves the service level in these high-failure-rate regions but also enhances the overall service level across all regions.

The article is organized as follows. Section 2 provides a comprehensive review of the related literature. In Section 3, we analyze historical operational data and identify two observations regarding the black hole phenomenon. Section 4 presents our stochastic periodic repair and maintenance model. In Section 5, we introduce an operational data analytics framework that improves data-integrated solutions by sequentially boosting a given oracle solution, such as one obtained from distributionally robust optimization. Section 6 discusses empirical results and practical insights derived from various numerical studies. Finally, we conclude the article in Section 7.

2. Literature review

This section provides a comprehensive review of the relevant literature on bike-sharing repair, joint replenishment, and small-sample data-driven optimization.

3. Black hole phenomenon: Data and observations

We obtained a comprehensive dataset from a leading industry partner in China, comprising operational records for 4,904 electric bikes in 2022. This dataset consists of order, battery swapping, user-reported failure, and repair data (see Appendix EC.1.1 for details). This section presents two key empirical observations from the repair data that motivate our modeling choices and reveal important operational challenges.

3.1. Failure rate differences in regions

To examine regional reliability, we define the regional failure rate s i as the expected proportion of bikes in region i transitioning to a failed state per unit time, estimated from operational data by dividing the number of new failures by the number of deployed bikes and the observation period length. Further details on interpretation and estimation of s i are provided in Appendix EC.1.5. Figure 1(a) shows the spatial distribution of user-reported failure records over a two-month period. However, absolute counts of reported failures can be misleading, as higher numbers may simply reflect larger bike populations rather than higher failure rates.

OPEN IN VIEWER

Figure 1.

Failure rate differences in regions.

To address this, we divide the entire operational area into 12 well-defined service regions and calculate failure rates for each. Figure 1(b) presents these region-specific failure rates. Notably, four adjacent regions exhibit considerably higher failure rates, suggesting localized operational issues where bikes are more prone to breakdowns. This pattern illustrates what we refer to as the “black hole” phenomenon, where certain areas consistently experience significantly high failure intensities. Notably, this phenomenon is not binary, but rather reflects a continuous spatial variation in failure rates across regions. To assess the generality of this finding, we analyze data from two additional bike-sharing firms in different cities, with Appendix EC.1.4 documenting similar region-specific failure patterns that support such heterogeneous failure structures. These observations underscore the need to identify and manage high-failure-rate regions in large-scale bike-sharing systems. Motivated by this, we integrate failure rates as a pivotal component in our subsequent modeling framework.

3.2. Impact of covariate information

In this section, we investigate how observable covariate information influences failure rate distributions in bike-sharing systems. We begin by constructing a set of eight covariates derived from raw operational data, including historical ride activity, repair logs, and bike age (detailed in Appendix EC.1.2). In addition to internal operational factors, we incorporate external factors such as temperature, precipitation, and wind, which may affect both bike usage and mechanical stress. We also include the “day of the week,” as usage patterns often vary significantly across weekdays and weekends. To identify the most relevant covariates, we apply Lasso regression with L 1 regularization and select the tuning parameter via cross-validation (Tibshirani, 1996). The Lasso regression yields an R 2 of 0.60, indicating a reasonably good fit to the data. Seven covariates are retained, with results and selected variables reported in Appendix EC.1.3.

Given the heterogeneity observed in failure patterns, we distinguish between deterministic covariates known in advance (e.g., seasonal patterns or long-term repair histories) and stochastic covariates realized over short operational periods (e.g., weather conditions or recent usage intensity). Deterministic covariates capture predictable operational cycles and influence high-level system planning, whereas stochastic covariates drive variability in failure rates and form the basis for scenario analysis. To analyze the impact of covariates, we partition the covariate space into distinct scenarios, each representing a specific operational regime. For instance, failure risks under extreme weather and high usage differ significantly from those under mild weather and low usage. Consequently, partitioning the covariate space into scenarios enables us to condition ambiguity sets on observable signals, thereby improving the structure of uncertainty modeling. Our approach aligns with recent advances in DRO that emphasize the value of covariate information in shaping uncertainty (e.g., Hao et al., 2020; Perakis et al., 2023).

Based on the seven selected covariates, we construct a multivariate regression tree to classify failure rate scenarios. Each leaf node of the tree represents a scenario, defined as a subset of covariate realizations. Figure 2 illustrates two representative scenarios identified from the tree. For each scenario, the figure shows the spatial distribution of failure rates across regions, suggesting that failure behavior varies meaningfully across scenarios. To rigorously confirm the significant difference in failure rates between the two scenarios, we conduct Hotelling’s T 2 test, a multivariate extension of the t-test (Hotelling, 1931). Given potential violations of normality, we perform a permutation test with 5,000 iterations based on Hotelling’s T 2 statistic. The observed T 2 statistic of 51.42 ( F = 3.04 , p = 0.0079 ) indicates statistical significance at the 1% level, confirming that failure rate distributions differ across scenarios. These findings provide strong support for a scenario-based modeling approach and underline the value of incorporating both deterministic and stochastic covariate information when analyzing failure rate heterogeneity.

OPEN IN VIEWER

Figure 2.

Failure rate pattern in various scenarios.

In this section, we formulate a data-driven repair and maintenance problem for a free-floating bike-sharing company operating over a defined service area. The problem is modeled on a complete directed graph G = ( I ∪ 0 , A ) , where I denotes regions to be served, 0 denotes the repair depot, and A = { ( i , 0 ) : i ∈ I } denotes the arcs between them. The distance associated with arc ( i , 0 ) ∈ A is denoted by d i , and the number of bikes in region i is denoted by N i . Opening a new repair route incurs a fixed cost denoted by ρ . Let J = { 1 , … , | J | } denote the index set of potential repair routes. This index set does not correspond to specific spatial paths in advance but rather serves as a placeholder for the set of routes to be selected.

Fully dynamic multi-period models can capture stochastic system evolution more accurately, but they are often computationally infeasible in practice. For example, Feng et al. (2015) formulate the joint replenishment problem as a Markov decision process and report that optimal policies are only computable for instances with few product types due to the curse of dimensionality. In light of these challenges, we adopt a periodic repair policy as a practical and computationally tractable alternative. Under this policy, each repair route j ∈ J is assigned a repair period T j , which is treated as a decision variable rather than a fixed input. All regions assigned to the same route are repaired at the same frequency, which enables coordinated routing and the consolidation of repair operations. The model determines how frequently each group of regions should be repaired by jointly optimizing these periods. Our approach builds on a well-established literature on periodic policies in inventory control, motivated by their analytical tractability and supported by strong theoretical and empirical performance; see Section 2.2 for a concise review. By incorporating repair periods as decision variables, our approach preserves key features of dynamic planning while ensuring that the resulting optimization problem remains computationally tractable, as detailed in Appendix EC.4. In the same appendix, we also present numerical evidence showing that periodic repair policies achieve small optimality gaps (below 2.5%) compared to fully dynamic multi-period optimal solutions under a stylized setting (i.e., a single-route, single-scenario case with a two-point failure-rate distribution). However, we explicitly acknowledge that this numerical study does not establish a theoretical relationship between the periodic cost rate and the fully dynamic optimal cost. Providing a rigorous theoretical guarantee for this approximation remains an open question for future research.

The objective is to determine the optimal number of repair routes (up to | J | ), repair assignments y i j , repair periods T j , and preventive maintenance interventions z i to minimize the total cost rate—that is, total cost per unit time—while covering all regions. For simplicity, all broken bikes in a region are assigned to the same route. Once the repair time T j arrives, a truck departs the depot to service broken bikes. In our industry partner’s operations, trucks perform on-site repairs for bikes with minor malfunctions, while those with more substantial issues are transported to depots for comprehensive repairs. Based on feedback from our industry partner, in-depot repairs rarely exceed truck capacity; thus, capacity constraints are omitted. Let c r denote the repair cost per bike and γ the truck routing cost per unit distance. The truck’s routing length consists of two parts: the shortest tour length l ( y j ) for each truck j , and the distance traveled within regions to locate broken bikes. The first part is an NP-hard traveling salesman problem (TSP), which we approximate using a linear model (Beardwood et al., 1959; Liu et al., 2021) to ensure compatibility with our optimization framework. This approach employs tractable predictors and approximation models for tour length, omitting the specific visiting sequence. Details are in Appendix EC.5.1. The second part scales linearly with the number of broken bikes, with coefficient ϕ . With slight abuse of notation, we represent failure rates by the random vector s = ( s i ) , where s i denotes the failure rate in region i , defined as the proportion of newly malfunctioning bikes per unit time. Based on the empirical findings in Section 3.2, failure rates may exhibit dependence on observable covariates. We represent these covariates by a V -dimensional vector v = ( v L , v S ) , where v L includes long-term covariates known in advance, and v S captures short-term covariates realized over operational cycles. We denote the joint failure rates and covariates as ( s , v ) ∈ R | I | × R V .

To mitigate the black hole phenomenon, we propose proactive preventive maintenance, a well-established practice in reliability literature aimed at preventing failures and extending the system lifespan of capital-intensive systems, thus lowering long-term costs (Li and Xu, 2004). The industry has also recognized its value and begun incorporating it into practical operations (Velco, 2022). The baseline failure rate of a repairable system can be decomposed into recoverable damage s i 1 and intrinsic fatigue s i 2 (Sun et al., 2018), with s i = s i 1 + s i 2 for each i ∈ I . Recoverable damage typically encompasses minor damage and can be rectified through preventive maintenance intervention; examples include brake failures and loose seats. Conversely, intrinsic fatigue represents the inevitable accumulation of damage and cannot be mitigated through preventive maintenance intervention; examples include tire damage or frame deformations. In practical operations, real-time preventive maintenance can be realized by rebalancing workers with a cost of c p per unit of time. We assume a uniform effectiveness ratio ω = s i 1 / s i across all regions, implying that preventive maintenance yields a proportional reduction of ω s i in the failure rate for each region i . This modeling choice is motivated by two considerations. First, preventive maintenance interventions—such as brake adjustments and seat tightening—are standardized and performed by trained personnel following uniform protocols, which supports a consistent effectiveness ratio. Correspondingly, in our model, preventive maintenance is treated as a system-wide intervention that proportionally reduces the hazard rate of operational bikes, following a proportional hazards model (Kumar and Klefsjö, 1994). Second, our empirical data lack sufficient granularity to reliably estimate region-specific effectiveness ratios ω i , given the limited availability of observed preventive maintenance outcomes. Accordingly, we model ω as a global parameter that strikes a balance between model tractability and empirical support.

To this end, we make four types of decisions: the number of repair routes, the assignment of repair routes to bike-sharing regions, the length of the corresponding repair period, and where to provide preventive maintenance. Specifically, we define the following decision variables:

x j : a binary variable, x j = 1 if the repair route j is selected, 0 otherwise;

We now describe the decision-making process. Figure 3 shows the sequence of periodic repair and maintenance operations. In the figure, the decisions are depicted using dashed lines, and information realizations are depicted with solid lines. Section 3.2 provides empirical evidence that both deterministic and stochastic covariate information can significantly influence failure rate behavior. Leveraging the information contained in the long-term covariates (e.g., seasonality indicators or long-term repair histories), the firm first determines the optimal number of selected repair routes ∑ j ∈ J x j and the corresponding repair assignment y i j between regions and repair routes. Subsequently, at the beginning of each operational cycle, short-term stochastic covariates (e.g., weather conditions, recent failure reports, and usage forecasts) are realized, where an operational cycle may correspond to a day, a week, a month, or even a quarter, depending on the company’s staffing structure and operational policies. After observing these short-term covariates, the firm determines the repair period T j for each repair route j and the real-time preventive maintenance intervention z i for each region i , where T j ≥ 0 for all j ∈ J and z i ∈ 0 , 1 for all i ∈ I . Finally, the uncertain regional failure rates are realized, which in turn determine the realized system outcomes under the implemented decisions.

OPEN IN VIEWER

Figure 3.

Sequence of periodic repair and maintenance operations.

This implies that the variables T and z are adaptive decision variables, specifically T ( v ) and z ( v ) , respectively. The decision variables T and z correspond to the two special strategies: dedicated repair period and preventive maintenance intervention, which are designed to mitigate the black hole phenomenon. The effectiveness of these strategies is first theoretically demonstrated within a two-region system in Section 4.1 and further numerically validated in the general system in Section 6. Note that although x and y could be viewed as adaptive to the long-term covariates v L , these covariates are fully known deterministic features, introducing no uncertainty. Therefore, we omit the explicit adaptive notation for x and y , embedding the deterministic covariates directly into the long-term decision-making stage.

At the beginning of each operational cycle—which may correspond to a day, a week, a month, or even a quarter—the operator determines a repair plan that includes both repair periods and preventive maintenance decisions, based on observed covariate information. The length of the operational cycle can be determined flexibly in practice. For firms characterized by high operational flexibility, urgent repair demands (as exemplified by our industry partner operating an e-bike sharing system), and sufficient workforce capacity, it may be feasible to update repair plans on a daily basis. In contrast, firms with more rigid scheduling, extended repair deadlines, or limited workforce capacity may prefer to revise repair decisions on a weekly or monthly basis. Our model accommodates this flexibility by treating the operational cycle as a generic time unit, enabling implementation across diverse operational environments.

Furthermore, the inter-cycle dependencies are relatively weak. Each cycle effectively resets the operational condition of the serviced bikes. In other words, the system “renews” itself through corrective repairs. As a result, the residual effect of one cycle on the next can be reasonably neglected. We intentionally adopt a cycle-based framework to reflect how operational decisions are commonly structured in practice. Many bike-sharing service providers plan repair activities on a rolling basis, leveraging the most up-to-date information to dynamically adjust decisions in response to evolving covariates—such as external factors (e.g., weather and seasonality) and observed failure patterns. Our model is designed to support such operational planning by balancing realism and computational tractability.

Our focus in this article is to provide a prescriptive analytics framework for the periodic repair and maintenance problem. Let P 0 ( R | I | × R V ) denote the set of all distributions of a random vector of dimension | I | + V . If the firm has perfect knowledge of the joint distribution of ( s , v ) , denoted by P ∈ P 0 ( R | I | × R V ) , we can formulate the periodic repair and maintenance problem as the following stochastic optimization problem.

min x , y ∑ j ∈ J ρ x j ⏟ (A) + E P [ Ψ ~ ( x , y , s , v ) ]

(1a)

s.t. ∑ j ∈ J y i j = 1 , ∀ i ∈ I ,

(1b)

y i j ≤ x j , ∀ i ∈ I , j ∈ J ,

(1c)

y i j ∈ { 0 , 1 } , ∀ i ∈ I , j ∈ J ,

(1d)

x j ∈ { 0 , 1 } , ∀ j ∈ J ,

(1e)

Ψ ~ ( x , y , s , v ) = min z ∈ { 0 , 1 } | I | T ≥ 0 ∑ i ∈ I c p N i z i ( v ) ⏟ (B) + ∑ i ∈ I δ N i ( ∑ j ∈ J s i ( 1 − ω z i ( v ) ) T j ( v ) y i j − C ) + ⏟ (C) + { ∑ j ∈ J γ l ( y j ) + ∑ i ∈ I ( γ ϕ N i s i ( 1 − ω z i ( v ) ) T j ( v ) + c r N i s i ( 1 − ω z i ( v ) ) T j ( v ) ) y i j T j ( v ) } ⏟ (D): the first two terms, and (E): the last term .

(2)

The objective function (1a) minimizes the expected total cost per unit of time; it is composed of five parts: (A) the fixed cost rate of repair route, (B) the cost rate of preventive maintenance, (C) the penalty cost rate of broken shared bikes, (D) the routing cost rate to collect the broken shared bikes, and (E) the repair cost rate. Note that when the broken bike proportion in a region exceeds a specified threshold, denoted by C , the user experience may be adversely affected due to reduced bike availability. To quantify this effect, we introduce a unit penalty coefficient δ to account for this dissatisfaction. Specifically, δ represents the penalty cost per unit of time for each bike when the broken bike proportion in a region exceeds the threshold C . Consequently, the penalty cost rate for region i , assigned to route j over the repair period T j , is given by δ N i ( s i ( 1 − ω z i ) T j − C ) + . This penalty structure can make high-failure-rate regions more costly when pooled with low-failure-rate regions under a shared (and typically longer) repair period, thereby implicitly capturing the operational mechanism associated with the “black hole” effect. Constraints (1b)–(1c) enforce that each region i is served by exactly one repair route, while constraints (1d)–(1e) specify the decision variable domains. This setting is typical in the uncapacitated facility location problem, which is a notorious NP-hard problem.

Cycle-based planning horizon: Recall that the full set of covariates is denoted by v = ( v L , v S ) , where v L represents long-term deterministic covariates known prior to the planning horizon. Accordingly, the long-term decisions x and y are made after observing v L and remain fixed over a planning horizon that spans multiple operational cycles (e.g., a week or a month). These deterministic covariates are therefore embedded directly into the long-term planning stage. In contrast, short-term covariates v S are stochastic and realized at the beginning of each operational cycle. After their realization, adaptive decisions T ( v ) and z ( v ) are made based on the full covariate information. This structure enables dynamic adjustment of short-term decisions in response to updated conditions, while preserving tractability in long-term planning.

Cluster-adapted policy: The decision variables T and z in the adaptive stochastic program (1) adapt to the covariate information v , while x and y , although conceptually adaptive to v L , are treated as non-adaptive since v L is deterministic covariate. To derive a practicable optimization model, we introduce the following cluster-adapted policy proposed by Perakis et al. (2023), which adapts to the realized covariate information cluster.

M ≜ { T ^ : Ω → R | J | , z ^ : Ω → { 0 , 1 } | I | | T ^ ( v ) = T l , z ^ ( v ) = z l if v ∈ Ω l , l ∈ L } .

The corresponding function Ψ ( x , y , s , v ) of the objective function is defined as follows:

Ψ ( x , y , s , v ) = min z l ∈ { 0 , 1 } | I | T l ≥ 0 ∑ i ∈ I c p N i z l i + ∑ i ∈ I δ N i ( ∑ j ∈ J s i ( 1 − ω z l i ) T l j y i j − C ) + + ∑ j ∈ J { γ l ( y j ) + ∑ i ∈ I ( γ ϕ N i s i ( 1 − ω z l i ) T l j + c r N i s i ( 1 − ω z l i ) T l j ) y i j T l j } .

The stochastic optimization model (1), however, poses the following two challenges. 1. Distributional ambiguity: The stochastic optimization model (1) assumes the precise knowledge of the joint distribution of failure rates and covariate information. Accurate estimation of the distribution is typically infeasible, and ignoring this uncertainty can lead to the “optimizer’s curse” (Smith and Winkler, 2006). 2. Computational challenges: Even with perfect knowledge of the true joint distribution, stochastic optimization still suffers from the “curse of dimensionality,” rendering it computationally intractable. To address these challenges, we formulate an operational data analytics framework for data-integrated repair decisions in Section 5.

4.1. Illustration of two special strategies

Here, we illustrate two specific strategies to mitigate the black hole phenomenon, namely the dedicated repair period and preventive maintenance intervention. To facilitate our analysis, we consider a two-region empirical model, as depicted in Figure 4. This setup features two distinct failure rates, s 1 and s 2 , where s 1 denotes a higher failure rate (i.e., black hole). The number of scenarios | L | is fixed at 1 to facilitate a clearer definition of the black hole region. With slight abuse of notation, we define S = s 2 s 1 ∈ ( 0 , 1 ] as an inverse indicator of the black hole effect, meaning that a smaller value of S indicates a stronger black hole effect. Each region has N shared bikes, and the distance between the two regions is denoted by d . To obtain structural results, we set the threshold for the proportion of broken bikes C and the route-opening cost ρ to 0. This eliminates nonlinear penalties and removes cost-related barriers to opening a dedicated route for the black hole region. A repair depot is located between the two regions, at a distance of m d from each region. The objective is to determine the optimal assignment, repair period, and preventive maintenance to minimize the overall cost rate of the two-region system.

OPEN IN VIEWER

Figure 4.

Two assignment policies in a two-region system.

4.1.1. Illustration of dedicated repair period

This analysis aims to illustrate the effectiveness of dedicated repair period as a specific strategy to mitigate the black hole phenomenon. To focus on this strategy, we prohibit preventive maintenance by setting unit maintenance cost c p as infinity. For notational simplicity, we refer to the resulting optimization problem as Problem ( P 1 ). In the two-region system, there are two assignment policies: the pooling policy and the dedicated policy. Under the pooling policy, both regions are assigned to a repair route, resulting in a shared repair period, as shown in Figure 4(a). Conversely, the dedicated policy assigns a specific repair period to the black hole region, as shown in Figure 4(b). We highlight the key insights regarding the dedicated repair period strategy in the following proposition.

Proposition 1

When m = 1 2 , the dedicated policy is optimal. When m > 1 2 , if S ≤ S * , the dedicated policy is optimal, whereas if S > S * , the pooling policy is optimal, where S * satisfies 1 S * + S * = 4 m .

Proposition 1 provides key insights into the role of a dedicated repair period in addressing failure rate disparities between the two regions, represented by S . Specifically, when m > 1 2 , the decision to implement a dedicated repair policy is primarily determined by S . When the black hole phenomenon is pronounced, characterized by S ≤ S * , assigning a dedicated repair period to the black hole region is preferred. This preference stems from the fact that grouping regions with significant differences in failure rates together will significantly increase the penalty cost, as the repair period of black hole regions is prolonged, and broken bikes cannot be timely repaired. Conversely, if the black hole phenomenon is less severe, represented by S > S * , it becomes preferable to adopt the assignment scheme based on geographical location (pooling policy). In this case, the reduced penalty costs resulting from the dedicated policy fail to offset the increased routing costs. Thus, implementing a dedicated repair period for the black hole region may significantly contribute to minimizing the overall operational costs under certain conditions.

4.1.2. Illustration of preventive maintenance

Our analysis now focuses on preventive maintenance as a special strategy to mitigate the black hole phenomenon. In this context, we allow for preventive maintenance. For notational convenience, we refer to the corresponding optimization problem as Problem ( P 2 ). Regions with high failure rates could benefit from preventive maintenance, as it reduces these failure rates, thereby alleviating the impact of the black hole phenomenon. We highlight the key insights regarding preventive maintenance intervention in the following proposition.

Proposition 2

The optimal assignment policy for Problem ( P 2 ) satisfies the following properties: (a) When the pooling policy is optimal for Problem ( P 1 ), it remains optimal for Problem ( P 2 ); (b) When the dedicated policy is optimal for Problem ( P 1 ), the pooling policy becomes optimal for Problem ( P 2 ), if only the region with s 1 is maintained and the disparity in failure rates is small, characterized by S 1 − ω > S * ; otherwise, the dedicated policy remains optimal for Problem ( P 2 ).

Part (a) of Proposition 2 establishes that if the pooling policy is optimal for Problem ( P 1 ), the black hole region will not adopt a dedicated repair period, even when preventive maintenance intervention is introduced in Problem ( P 2 ). Part (b) illustrates that when the failure rate of the black hole region ( s 1 ) is not significantly higher than that of the non-black hole region ( s 2 ), represented by S 1 − ω > S * , and only black hole region undergoes preventive maintenance in Problem ( P 2 ), the region will no longer adopt a dedicated repair period due to the reduction in its failure rate. These findings underscore the pivotal role of preventive maintenance intervention in effectively mitigating the black hole phenomenon and guide managerial decisions on the adoption of dedicated repair periods for black hole regions.

In summary, we remark that the assumptions of a two-region system with deterministic failure rates are crucial for our characterization of the two special strategies in mitigating the black hole phenomenon. While this simplified system does not fully capture real-world complexities, analyzing it yields insights into the performance of strategies in more general systems, as discussed in subsequent sections.

5. Solution approach: An ODA-based boosting framework

In this section, we employ the operational data analytics (ODA) framework to structure the development of data-integrated repair decisions for bike-sharing systems. The motivation for using ODA is to directly link observed operational data to implementable decisions. From the ODA perspective, existing approaches such as predict-then-optimize (PTO) and distributionally robust optimization (DRO) are not competing methodologies; instead, they serve as data-integration modules that generate oracle solutions based on historical failure-rate and covariate information. Importantly, ODA does not rely on stronger or additional modeling assumptions than these approaches. Rather, it adopts a distinct integration-and-validation logic by treating decisions as adaptive functions of data—operational statistics—and selecting among them based on out-of-sample operational performance.

The ODA framework consists of two conceptual pillars. First, the data-integration model defines a class of admissible mappings from observed operational data, such as failure rates and covariates, to implementable repair decisions. This class incorporates partial structural knowledge of the system, such as homogeneity in a one-route setting or more flexible structures in general multi-route environments. In practice, this mapping class can be constructed by sequentially boosting a given oracle solution, thereby extending a single data-integrated rule into a structured family of operational statistics. Thus, ODA restricts the decision space to structurally meaningful data-to-decision mappings that reflect how operational signals should inform repair policies.

Second, the decision validation model selects among these candidate data-to-decision mappings by evaluating their resulting operational performance, such as the expected cost rate of periodic repair. Unlike approaches that optimize modified objectives (e.g., a worst-case criterion induced by an ambiguity set), ODA directly compares candidate data-integrated solutions based on validation performance and selects the one that performs best. This validation-driven perspective is particularly useful in our setting, as it can substantially improve solution performance when the available data are limited.

Together, these two components naturally lead to a sequential boosting-and-validation procedure. Starting from a baseline data-integrated solution generated by PTO, DRO, or other data-driven approaches, ODA constructs a structured subclass of operational statistics and employs validation to identify the solution with the best performance. We first analyze a homogeneous one-route case and then extend the analysis to develop an enhanced data-integrated solution for the general multi-route setting.

5.1. A homogeneous one-route case

In this section, we consider a special case of the periodic repair problem to illustrate its performance advantage of the ODA framework over conventional data-integrated benchmarks. To ensure analytical clarity and to align with the homogeneous property commonly assumed in this framework, we restrict attention to a simplified one-route setting. In addition, the threshold for the proportion of broken bikes is set at C = 0 , and the preventive maintenance decision z is held constant. Consequently, the objective function of the stochastic optimization problem (1) reduces to:

Γ ( T , s , v ) = ρ + E P [ ∑ i ∈ I c p N i z i ( v ) + ( γ ϕ + c r ) N i s i ( 1 − ω z i ( v ) ) + γ l T ( v ) + ∑ i ∈ I δ N i s i ( 1 − ω z i ( v ) ) T ( v ) ] = ρ + ∑ l ∈ L p l [ ∑ i ∈ I c p N i z l i + ( γ ϕ + c r ) N i s i ( 1 − ω z l i ) + γ l T l + ∑ i ∈ I δ N i s i ( 1 − ω z l i ) T l ] = ρ + ∑ l ∈ L p l [ ∑ i ∈ I c p N i z l i + ( γ ϕ + c r ) N i s i ( 1 − ω z l i ) ] ⏟ Γ 1 ( s , v ) + ∑ l ∈ L p l [ γ l T l + ∑ i ∈ I δ N i s i ( 1 − ω z l i ) T l ] ⏟ Γ 2 ( T , s , v ) .

The second equation holds under the adoption of the cluster-adapted policy. Before developing the solution approach, we first examine how the optimal decision should relate to the observed data. Notably, the function Γ 1 ( s , v ) is independent of the repair periods T , whereas the function Γ 2 ( T , s , v ) exhibits an appealing homogeneous property: Γ 2 ( c 0 − 1 2 T , c 0 s , v ) = c 0 1 2 Γ 2 ( T , s , v ) . This relation implies that if the observed failure rates s are scaled by a factor of c 0 , then the optimal repair periods T should adjust proportionally by a factor of c 0 − 1 2 . This insight motivates a desirable property for the data-integrated decision: it should preserve this scaling relationship between the repair periods and the observed failure rates.

Proposition 3

A statistic T : R | I | × R V → R + | L | satisfies to:

T l * ( c 0 s , v ) = c 0 − 1 2 T l * ( s , v ) , ∀ c 0 ∈ R + , v ∈ Ω l , l ∈ L .

We define the following class of homogeneous operational statistics as

Q H = { T ∈ R + | L | : T l * ( c 0 s , v ) = c 0 − 1 2 T l * ( s , v ) , ∀ c 0 ∈ R + , v ∈ Ω l , l ∈ L } .

(3)

This scaling property guides the mapping of data into decisions, leading to the data-integration model. Importantly, no alternative data-integrated solution uniformly dominates this class of homogeneous operational statistics (Feng and Shanthikumar, 2023, Theorem 2). Accordingly, we restrict the data-integration model to this class when formulating the ODA framework. The decision validation model then selects the most effective data-integrated solution by evaluating its ultimate performance, measured by the expected cost associated with implementing the corresponding operational statistic.

Proposition 3 motivates the use of sequential boosting to further enhance the obtained solutions. In particular, the given oracle solutions generated by existing data-driven approaches can be viewed as baseline rules that may still leave room for systematic improvement. Sequential boosting builds directly on these existing approaches by offering a principled way to refine them iteratively. Specifically, starting from a candidate solution T c , we identify a subclass of homogeneous operational statistics within Q H that offers potential performance improvements over T c . Such a candidate solution can be obtained using any of the existing data-integration approaches in Section 5.2.1. We formulate the data-integration model as the boosted class of operational statistics based on the candidate solution:

Q B ( T c ) = { T ∈ R + | L | : T l ( s , v ) = β 0 T l c ( s , v ) , β 0 ∈ R + , v ∈ Ω l , l ∈ L } .

(4)

It is straightforward to verify that if the candidate solution T c ∈ Q H , then the corresponding boosted class is necessarily a subset of the class of homogeneous operational statistics; that is, Q B ( T c ) ⊆ Q H . The decision validation model then assesses each operational statistic based on its ultimate performance, measured by the expected implementation cost.

min T ∈ Q B ( T c ) E P [ Γ ( T , s , v ) ] .

(5)

Due to the lack of knowledge of P , direct evaluation of the expected objective may be infeasible. Instead, we employ a sample average approximation to assess decisions within the data-integration model. Operating in a nonparametric setting, we construct the model by sequentially boosting a baseline solution derived from existing approaches. We employ validation-based hyperparameter optimization to determine the optimal boosting parameter β 0 * and assess the performance of the resulting boosted solution. The following stylized example demonstrates the effectiveness of the proposed ODA framework in improving solutions for problems exhibiting the homogeneous property.

Example 1

Consider a case where s is one-dimensional and no covariate information is available. Suppose the true distribution of s is discrete with P [ s = 0.01 ] = 2 5 and P [ s = 0.10 ] = 3 5 , from which a training sample { 0.01 , 0.01 , 0.01 , 0.10 , 0.10 } is drawn. Using validation-based hyperparameter optimization, the optimal boosting parameter is determined to be β 0 * = 0.93 . Evaluated under the true distribution, the out-of-sample cost decreases from 21.86 to 21.69, showing that ODA improves solution quality for homogeneous problems.

5.2. The general multi-route case

The homogeneous one-route analysis in the previous section highlights the importance of selecting data-integrated solutions based on their realized operational performance, a central principle of the ODA framework. We now turn to the general multi-route setting, in which the objective function no longer exhibits homogeneity. Specifically, when the observed failure rates s are scaled by a factor c 0 , the corresponding optimal repair decisions do not scale proportionally. Instead, the resulting adjustments depend on the interaction among the scaling factor c 0 , the failure-rate realizations s , and the underlying operational structure. To accommodate this lack of homogeneity, we characterize operational policies through a general class of non-homogeneous operational statistics defined via admissible structure-preserving transformations. Formally, let π ( s , v ) = ( x , y , z , T ) denote an operational policy, and define the class of non-homogeneous operational statistics as

Q N H = { π : π ( c 0 s , v ) = χ c 0 , s , v ( π ( s , v ) ) , ∀ c 0 ∈ R + , v ∈ Ω l , l ∈ L } ,

where χ c 0 , s , v ( ⋅ ) denotes an admissible transformation, that preserves feasibility with respect to both continuous and binary decision variables.

5.2.1. The data-integration model

We begin with a baseline solution π c = ( x c , y c , z c , T c ) , derived from an existing data-driven approach such as PTO or DRO. Following the ODA construction, we generate data-integrated solutions by applying multiplicative scaling rules to the baseline solution. For continuous decisions, we adopt the standard multiplicative form:

T l j = β l j T T l j c , β l j T ∈ B l j T ,

(6)

where β T = ( β l j T ) is a nonnegative scaling matrix, and B T denotes a prespecified admissible family. This formulation preserves feasibility and enables refinements of candidate continuous decisions. For binary decisions, direct multiplicative scaling is generally incompatible with integrality requirements. To maintain consistency with the ODA framework while preserving integrality feasibility, we introduce the following structure-preserving transformation for binary decisions. Specifically, for each binary variable b ∈ { x j , y i j , z l i } , we define

b = β b b c + ( 1 − β b ) ( 1 − b c ) , β b ∈ B b ,

(7)

where β = ( β b ) denotes the binary scaling parameters, and B = ∏ b B b denotes the admissible family of binary scaling transformations. This transformation can be interpreted as a restricted multiplicative adjustment subject to integrality constraints, which maps a candidate binary decision to another feasible binary value and thereby extends the ODA scaling principle to binary settings. Collectively, these multiplicative transformations define the class of data-integrated solutions based on the baseline data-integrated solution π c :

Q B ( π c ) = { ( x , y , z , T ) : ( x , y , z , T ) as defined by the boosted rules (6) and (7) , ( β , β T ) ∈ B × B T , ( x , y ) satisfy (1b)-(1e) } .

(8)

This class provides a structured search space within which ODA selects improved solutions based on validation performance. While any data-integrated solution π c in Q N H does not satisfy the structured scaling relationship between the repair decisions and the observed failure rates, we can still identify a subclass of non-homogeneous operational statistics within Q N H . That is, the corresponding boosted class remains a subset of the general non-homogeneous operational statistics, that is, Q B ( π c ) ⊆ Q N H . While the boosted class Q B ( π c ) provides a flexible family of data-integrated solutions, it may be excessively large when used as a validation domain. A key principle of ODA is to refine the admissible class by incorporating structural insights derived from domain knowledge. In our setting, two such structural insights play a critical role.

Bounded adjustment of region-to-route assignments. In practice, baseline region-to-route assignments inferred from historical data are typically informative. Accordingly, within the data-integration model, we control the flexibility of region-to-route assignments by limiting the number of regions whose assignments are allowed to deviate from the baseline solution to at most κ . This constraint limits the flexibility of the corresponding adjustment parameters while preserving local adaptability. Formally, the admissible data-to-decision mappings are required to satisfy

∑ i ∈ I I { ∃ j ∈ J s.t. y i j ≠ y i j c } ≤ κ .

Failure-rate-ordered preventive maintenance decisions. To further impose operationally meaningful structure on the admissible class of preventive maintenance decisions, we exploit a key insight obtained from a deterministic case. Specifically, when a single repair route is selected and all regions are assigned to it, the resulting problem admits a simple and interpretable structure. As established in Proposition EC.2 in Appendix EC.3.2, the optimal preventive maintenance policy follows a monotone threshold structure with respect to failure rates: regions with sufficiently high failure rates are selected for preventive maintenance, while those below the threshold are not. Equivalently, if a region with a higher failure rate does not receive preventive maintenance, then all regions with lower failure rates also forgo maintenance. Accordingly, within the ODA data-integration framework, we restrict the admissible class of preventive maintenance decisions to those that preserve this failure-rate-ordered structure. Let μ l i denote the mean failure rate of region i ∈ I under scenario l . For regions assigned to route j , that is, { i ∈ I : y i j = 1 } , order them such that μ l i 1 ≤ μ l i 2 ≤ ⋯ ≤ μ l i | I j | , where I j := { i ∈ I : y i j = 1 } . Formally, we restrict the preventive maintenance decisions to satisfy a failure-rate-ordered structure:

z l i 1 ≤ z l i 2 , ∀ i 1 , i 2 ∈ I j such that μ l i 1 ≤ μ l i 2 .

Thus, incorporating the above restrictions, we define the refined validation domain Q ~ B ( π c ) as

Q ~ B ( π c ) := { π ∈ Q B ( π c ) : ∑ i ∈ I I { ∃ j ∈ J s.t. y i j ≠ y i j c } ≤ κ , z l i 1 ≤ z l i 2 , ∀ l ∈ L , j ∈ J , i 1 , i 2 ∈ I j = { i ∈ I : y i j = 1 } such that μ l i 1 ≤ μ l i 2 } .

(9)

By construction, we have Q ~ B ( π c ) ⊆ Q B ( π c ) , ensuring that the refined domain is a proper subset of the original validation domain while incorporating both prior knowledge and operationally motivated structural constraints. Before introducing the decision validation model, we first present existing solution approaches that induce a baseline data-integrated solution and thereby provide initial elements of the admissible class defined by the data-integration model.

Predict and Then Optimize (PTO). A conventional approach is to first estimate the joint distribution P of ( s , v ) from historical data and then optimize the resulting stochastic program to derive the optimal repair decisions. Specifically, the theoretically optimal decision given the knowledge of the joint distribution P can be obtained by solving the stochastic program as follows.

min x , y , z l , T l ∑ j ∈ J ρ x j + ∑ l ∈ L p l [ ∑ i ∈ I c p N i z l i + δ N i ( ∑ j ∈ J E P [ s i | v ∈ Ω l ] ( 1 − ω z l i ) T l j y i j − C ) + ∑ ∑ + ∑ j ∈ J { γ l ( y j ) + ∑ i ∈ I ( γ ϕ N i E P [ s i | v ∈ Ω l ] ( 1 − ω z l i ) T l j + c r N i E P [ s i | v ∈ Ω l ] ( 1 − ω z l i ) T l j ) y i j T l j } ] s.t. Constraints (1b) − (1e) .

(10)

Since the true distribution P is unknown, we rely on historical data to construct empirical estimators. Several classical approaches—such as the maximum likelihood estimator and the method of moments—are commonly employed for this purpose. We estimate the conditional expectation E P [ s ∣ v ∈ Ω l ] using the sample average μ l of the failure rate in each scenario l ∈ L .

Distributionally Robust Optimization (DRO). To mitigate overfitting under limited data, we adopt a widely used DRO approach as a baseline data-integrated solution, which models uncertainty in ( s , v ) by assuming that the true distribution P lies in an ambiguity set F ⊆ P 0 ( R | I | × R V ) . We formally define the adaptive data-integration DRO model as follows.

min x , y ∑ j ∈ J ρ x j + max P ∈ F E P [ Ψ ( x , y , s , v ) ]

(11a)

s.t. Constraints (1b) − (1e) .

(11b)

In Problem (12), the choice of ambiguity set F critically affects both model tractability and solution quality. Motivated by the empirical findings in Section 3.2, which show that failure rates depend on covariates such as historical riding patterns, repair activities, and weather conditions, we incorporate covariate information into the distribution of uncertain failure rates via a scenario-wise ambiguity set (Chen et al., 2020).

F := { P ∈ P 0 ( R | I | × R V ) | ( s , v ) ∼ P E P [ s ∣ v ∈ Ω l ] = μ l , ∀ l ∈ L E P [ ϱ l ( s ) ∣ v ∈ Ω l ] ≤ σ l , ∀ l ∈ L P [ v ∈ Ω l ] = p l , ∀ l ∈ L P ( s ∈ Z l ∣ v ∈ Ω l ) = 1 , ∀ l ∈ L } ,

where the support set is defined as Z l := { s ∈ R | I | : s _ l ≤ s ≤ s ¯ l } . The covariate space is partitioned into | L | non-overlapping regions Ω l , creating | L | scenarios with probabilities P ( v ∈ Ω l ) = p l and ∑ l ∈ L p l = 1 . For the ambiguity set F , we incorporate the mean, generalized moments, and support of s within each scenario l . Generalized moments, represented by a convex function ϱ l ( s ) , provide flexible statistical characterizations of the uncertainty in s , encompassing measures such as variance and mean absolute deviation (MAD) (Wiesemann et al., 2014). In this study, we adopt the operational data analytics framework developed by Feng and Shanthikumar (2023) to guide the selection of moment-based ambiguity sets; see the empirical evaluation in Section 6.2.

To solve the DRO model, we adopt a standard dual reformulation for tractability. The associated lemma and derivation are provided in Appendix EC.2.3, as they are not central to the ODA framework. The resulting problem (EC.6) is a semi-infinite program with infinitely many constraints and is not directly solvable. This nonlinear term ( ∑ j ∈ J s i ( 1 − ω z l i ) T l j y i j − C ) + in Problem (EC.6) prevents us from deriving a tractable robust counterpart. This penalty cost term, which adapts to the random variables s , can be represented by a function Υ l i ( s ) = ( ∑ j ∈ J s i ( 1 − ω z l i ) T l j y i j − C ) + for each l ∈ L , i ∈ I . To approximate the adaptive decision Υ l i ( s ) , we introduce an auxiliary random vector u and approximate Υ l i ( s , u ) by a linear decision rule, subject to the constraint Υ l i ( s , u ) ≥ ( ∑ j ∈ J s i ( 1 − ω z l i ) T l j y i j − C ) + , following Bertsimas et al. (2019) and Chen et al. (2020).

Y = { Υ ( ⋅ ) ∈ R | I | × R | I | → R | L | × | I | | Υ l i ( s , u ) = Υ l i 0 + ∑ k ∈ I Υ l i k 1 s k + ∑ k ∈ I Υ l i k 2 u k , ∀ l ∈ L , i ∈ I } ,

where all coefficients Υ l i 0 , Υ l i k 1 , Υ l i k 2 ∈ R . With the defined mapping functions for adaptive decisions, we now present the approximate adaptive DRO model as follows.

min ∑ j ∈ J ρ x j + ∑ l ∈ L ( α l + η l ′ μ l + θ l ′ σ l )

(12a)

s.t. Constraints (1b) − (1e) , (EC.6d) − (EC.6e) ,

(12b)

α l + η l ′ s + θ l ′ u ≥ p l [ ∑ i ∈ I c p N i z l i + ( γ ϕ + c r ) N i s i ( 1 − ω z l i ) + ∑ j ∈ J γ l ( y j ) T l j + ∑ i ∈ I δ N i Υ l i ( s , u ) ] , ∀ ( s , u ) ∈ W l , l ∈ L ,

(12c)

Υ l i ( s , u ) ≥ ∑ j ∈ J s i ( 1 − ω z l i ) T l j y i j − C , ∀ ( s , u ) ∈ W l , l ∈ L , i ∈ I ,

(12d)

Υ l i ( s , u ) ≥ 0 , ∀ ( s , u ) ∈ W l , l ∈ L , i ∈ I ,

(12e)

Υ ( ⋅ ) ∈ Y .

(12f)

Since all constraints in Problem (12) are linear in the random variables ( s , u ) , it can be reformulated as a tractable problem via standard duality techniques (He et al., 2020). Interestingly, even though we linearize the term Υ l i ( s , u ) using the LDR approximation, the following proposition states that the LDR approach still provides an optimal solution to Problem (11).

Proposition 4

Let Π and Π ′ represent the optimal values of Problem (11) and Problem (12), respectively. Thus, we have Π = Π ′ .

Since the seminal work of Ben-Tal et al. (2004), LDRs have been recognized as computationally efficient yet generally approximate, with exact optimality rarely guaranteed and typically dependent on specific structural properties. Subsequent studies identify conditions under which LDRs are exact: Bertsimas et al. (2010) and Iancu et al. (2013) for support-only ambiguity sets with structural assumptions; Bertsimas et al. (2019) and He et al. (2020) for two-stage problems with one- and multi-dimensional recourse; and Georghiou et al. (2026) for pointwise optimality, in contrast to the worst-case perspective of our study. Most results focus on single-scenario problems, with Hao et al. (2020) being a notable exception showing LDR optimality under scenario-wise ambiguity sets for multidimensional recourse. Our work contributes to this literature by identifying an appealing structure for our periodic repair problem with a scenario-wise ambiguity set under which the LDR approximation is provably optimal.

Note that the dual formulation of Problem (12) cannot be directly solved by commercial solvers (e.g., CPLEX and Gurobi) due to the nonlinear term γ l ( y j ) T l j in Constraint (12c). We reformulate Problem (12) as an equivalent quadratically constrained program (QCP) for fixed ( s , u ) ∈ W l (Appendix EC.3.1), and develop a tractable iterative constraint generation approach (Appendix EC.3.2) to address the QCP’s nonconvexity and resulting computational inefficiency.

5.2.2. The decision validation model

The decision validation model subsequently evaluates each operational statistic based on its ultimate performance. While the homogeneous property no longer holds, the key principle of ODA—selecting the most effective data-integrated solution through decision validation—remains applicable. In this sense, ODA serves as a unifying meta-framework: it systematically improves solution quality by iteratively refining candidate decisions and validating their empirical performance. We next demonstrate that the ODA framework continues to improve solution quality even for problems that inherently lack the homogeneous property, as illustrated in Section 6.2.

6. Numerical studies

This section discusses numerical studies conducted using real-world data provided by our industry partner, who manages a bike-sharing system in China. The entire operational area is divided into 12 well-defined service regions, with the repair depot marked in red, as shown in Figure 5. In the first experiment, we evaluate our proposed scenario-wise DRO model by comparing its performance to (a) the widely used PTO approach, and (b) the DRO model with a nominal ambiguity set, referred to as the nonadaptive DRO (NA-DRO). The second experiment assesses the performance of the ODA framework in enhancing oracle data-integrated solutions. Furthermore, we quantify the value of black hole information and preventive maintenance by comparing our model with two benchmark models, and we measure regional disparity in penalty costs using the Gini coefficient. Finally, we demonstrate the effectiveness of the dedicated repair period.

OPEN IN VIEWER

Figure 5.

Illustration of repair depot and bike-sharing regions.

We select failure rate data over a 61-day period (the length of the operational cycle aligns with our industry partner) in June and July 2022. To evaluate model performance, we randomly partition the data into training and test sets, using two-thirds for training and one-third for testing, respectively. The parameter settings are provided in Table EC.4 in Appendix EC.5.3. All numerical studies are conducted on a Dell desktop equipped with a 2.50 GHz Intel i5 − 14400F CPU, 16 GB of memory, and the Windows operating system. The algorithm is implemented in Python using Gurobi 10.0 as the solver. To construct the scenario-wise ambiguity set, we use a multivariate regression tree to partition the covariate samples in the training set into several scenarios. The probability, mean, variance, and mean absolute deviation of the failure rate for each scenario l are estimated as described in Appendix EC.5.2.

6.1. Impact of covariate information and robustness

In this section, we assess the advantages of incorporating covariate information and DRO by contrasting the out-of-sample performance of our scenario-wise adaptive DRO model against two benchmark models. One benchmark is the PTO model based on sample averages, while the other is the nominal DRO model with a single scenario (i.e., | L | = 1 ). For ease of exposition, we use DRO − 2-Var and DRO − 4-Var to denote the variance-based DRO models, and DRO − 2-MAD and DRO − 4-MAD to denote the MAD-based DRO models, using 2 and 4 scenarios, respectively. Similarly, NA-PTO, PTO − 2, and PTO − 4 denote the PTO models with 1, 2, and 4 scenarios, respectively. We enumerate δ from 1 to 5 with step 0.5.

Figures 6(a) and 6(b) depict the mean and standard deviation of out-of-sample costs across the models, respectively. These results are normalized by dividing them by those of the NA-PTO model for comparative analysis. First, incorporating covariate information leads to improved out-of-sample performance, as evidenced by DRO − 2-Var (DRO − 2-MAD) outperforming NA-DRO-Var (NA-DRO-MAD) in both mean cost and stability. Interestingly, DRO − 4-Var underperforms DRO − 2-Var, indicating that increasing the number of scenarios does not necessarily enhance robustness. These results highlight the importance of selecting an appropriate scenario size, which can be guided by cross-validation. Second, all robust models outperform all PTO models, demonstrating the value of accounting for distributional uncertainty in failure rates, despite the increased complexity compared to the PTO models. Third, the out-of-sample performance differs substantially across DRO models with varying ambiguity set structures. Specifically, DRO-Var models tend to achieve lower mean costs, while DRO-MAD models exhibit lower variability.

OPEN IN VIEWER

Figure 6.

Performance comparison (measured out-of-sample) of DRO, NA-DRO, PTO, and NA-PTO models.

6.2. Performance improvement with sequential boosting

The ODA framework enhances the quality of oracle solutions by selecting the most effective data-integrated solution based on the decision validation model. We implement validation-based hyperparameter optimization to assess the out-of-sample performance of each boosted solution and identify the optimal boosting parameter. To ensure the reliability of our conclusions, we consider (i) the maximum number of potential repair routes | J | ∈ { 1 , 2 , 3 } , and (ii) the number of scenarios | L | ∈ { 1 , 2 , 3 } , yielding a total of nine problem instances. Table 1 reports the performance of various data-integration models across these instances.

OPEN IN VIEWER

Table 1.

Improvement in out-of-sample mean cost under sequential boosting.

			ODA:PTO				ODA:DRO-MAD				ODA:DRO-Var
| J |	| L |	PTO	κ = 0	κ = 1	κ = 2	DRO-MAD	κ = 0	κ = 1	κ = 2	DRO-Var	κ = 0	κ = 1	κ = 2
1	1	68.92	68.73	68.73	68.73	68.48	68.75	68.75	68.75	68.92	68.73	68.73	68.73
	2	67.54	66.89	66.89	66.89	68.51	67.99	67.99	67.99	67.54	66.89	66.89	66.89
	3	69.50	69.10	69.10	69.10	69.93	68.40	68.40	68.40	69.50	69.10	69.10	69.10
2	1	65.51	61.54	61.54	61.54	60.35	60.34	60.34	60.34	59.97	59.85	59.54	59.54
	2	64.95	60.01	60.01	60.01	59.14	59.07	59.07	59.07	59.34	58.91	58.91	58.91
	3	62.25	60.94	60.94	60.94	61.49	61.58	61.58	61.58	60.89	60.22	60.22	60.22
3	1	64.94	58.81	58.77	58.63	57.92	57.96	57.52	57.52	56.87	56.91	56.91	56.91
	2	62.17	59.49	59.46	59.43	57.07	57.55	58.82	58.82	57.89	57.65	57.51	58.81
	3	59.59	60.64	59.65	59.65	58.98	57.98	57.63	57.63	58.67	58.55	59.11	59.11
avg.		65.04	62.91	62.79	62.77	62.43	62.18	62.23	62.23	62.18	61.87	61.88	62.02

Several key observations emerge from the results of sequential boosting. First, sequential boosting generally enhances solution quality: the boosted solutions consistently yield lower out-of-sample mean costs than their conventional counterparts. This indicates that sequential boosting helps identify more effective repair period decisions that better balance service levels and operational costs. Second, even when the oracle data-integration PTO model underperforms (e.g., with an out-of-sample mean cost of 65.04), the ODA framework enhances solution quality (e.g., reducing the cost to 62.77) by data-driven selection of the optimal boosting parameter, further demonstrating its practical effectiveness. Third, as the reassignment budget κ increases, the marginal performance improvement of the ODA solution diminishes and may even become negative. This pattern indicates that ODA achieves its strongest performance gains when the data-integration model is carefully structured using domain knowledge. Fourth, the ODA solutions exhibit performance improvements when applied to both DRO models based on mean absolute deviation and variance, without one consistently outperforming the other. This underscores the value of the ODA framework in guiding the appropriate choice of moment-based ambiguity set structures. In subsequent numerical studies, we use the ODA framework to select the best-performing data-integrated DRO solutions.

6.3. Impact of black hole information and preventive maintenance

As shown in Figure 1(b) in Section 3.1, regions with high failure rates incur disproportionately higher penalty costs in the absence of targeted measures such as preventive maintenance or dedicated repair periods, leading to pronounced disparities—that is, “inequality”—across regions. To quantify this disparity, we use the Gini coefficient, a standard inequality metric ranging from 0 (perfect equality) to 1 (perfect inequality), commonly applied to income or wealth distributions but equally suitable for comparing regional penalty cost differences. We compare our proposed model against two counterparts: ODA:DRO-WPM (without preventive maintenance) and ODA:DRO-WBH (without black hole information). In ODA:DRO-WPM, preventive maintenance is disabled by setting the unit maintenance cost c p to infinity. In ODA:DRO-WBH, regional failure rate heterogeneity is removed by assigning all regions the average failure rate from the training data, representing a setting where spatial variability is ignored.

The Gini coefficient can be illustrated using the Lorenz curve. Figure 7(a) plots the Lorenz curve of penalty cost distributions for the case where | J | = 3 and | L | = 2 . The x-axis shows regions ranked by penalty costs, and the y-axis reports the corresponding cumulative share of total penalty costs. The ODA:DRO curve lies closer to the line of perfect equality, with a Gini coefficient of 0.43, compared to 0.51 and 0.70 for the benchmark models. This indicates that integrating both preventive maintenance and black hole information improves the equity of the penalty cost distribution. Preventive maintenance alone yields only a modest reduction in inequality, likely because it is deployed selectively in regions with extremely high failure rates (see Figure 8(a) in Section 6.4), often concentrated in a single region. By contrast, incorporating black hole information substantially narrows regional disparities, as it enables the operator to explicitly target persistently high-risk regions and mitigate excessive accumulation of failed bikes in localized areas. This result highlights that neglecting black holes causes a few regions to suffer most of the breakdowns, while accounting for them helps balance the load across regions and leads to a more equitable and spatially balanced system performance.

OPEN IN VIEWER

Figure 7.

Gini coefficient (measured out-of-sample) under different ODA:DRO models and dedicated policies.

OPEN IN VIEWER

Figure 8.

Repair plan and expected penalty cost (measured out-of-sample) at each region.

We further compare the out-of-sample statistics for component costs in Table 2, where bold values indicate superior performance. First, the ODA:DRO model achieves total cost rate reductions of 2.50% and 10.14%, respectively, relative to the two benchmark models, underscoring the value of incorporating both preventive maintenance and black hole information in bike-sharing repair operations. Second, while ODA:DRO diminishes penalty cost inequality across regions (see Figure 7(a)), it also notably reduces the overall penalty cost compared to ODA:DRO-WPM and ODA:DRO-WBH. This finding suggests that ODA:DRO not only enhances the service level in high-failure-rate regions but also delivers system-wide efficiency gains through more balanced resource allocation across regions. Third, interestingly, as shown in the column labeled Routing, the DRO model has a higher routing cost rate than the DRO-WBH model. This discrepancy appears to stem from DRO-WBH’s tendency to determine repair assignments only based on proximity due to the lack of black hole information, which results in a lower routing cost.

OPEN IN VIEWER

Table 2.

Out-of-sample statistics of component costs under different ODA:DRO models.

Model	Statistic	Opening	Routing	Maintenance	Repair	Penalty	Sum
ODA:DRO	Mean	3.00	38.85	2.00	4.99	7.62	56.46
	Var@90%	3.00	44.80	2.19	6.55	18.46	74.29
	Var@95%	3.00	45.81	2.19	6.89	19.50	76.35
	Standard deviation	0.00	5.15	0.27	1.28	7.22	12.84
ODA:DRO-WPM	Mean	3.00	40.89	0.00	5.08	8.94	57.91
	Var@90%	3.00	47.18	0.00	6.66	18.03	75.60
	Var@95%	3.00	48.12	0.00	6.97	22.82	84.53
	Standard deviation	0.00	5.62	0.00	1.32	9.14	14.72
ODA:DRO-WBH	Mean	2.00	35.00	10.87	4.67	10.28	62.83
	Var@90%	2.00	39.13	10.87	5.89	23.61	83.09
	Var@95%	2.00	40.31	10.87	6.29	26.24	83.92
	Standard deviation	0.00	3.92	0.00	1.18	12.76	17.34

6.4. Impact of dedicated repair period

Note that two key strategies are employed to address challenges in high-failure-rate regions: preventive maintenance and dedicated repair periods. Section 6.3 examines the impact of preventive maintenance on both the overall penalty cost (Table 2) and regional penalty cost disparities (Figure 7(a)). We now turn to repair periods and the resulting repair patterns. When | J | = 1 , all regions share a common repair period, whereas | J | > 1 allows high-failure-rate regions to adopt dedicated repair periods. When | J | = 2 , preventive maintenance is disabled by setting the unit maintenance cost c p to infinity, so that the effect of dedicated repair periods can be examined separately. For clarity, we focus on cases with | J | = 1 and | J | = 2 , corresponding to a single and a dedicated repair period, respectively. Results are reported for | L | = 1 ; those for | L | > 1 are omitted for brevity.

Figure 8 illustrates the relationship between the repair plan and assignment selection. For each region, the length of the adjacent green rectangle represents its repair period, while the radius of the red circle below indicates its expected penalty cost. Regions are marked with a red icon if preventive maintenance is applied and a blue icon otherwise. The figure shows that when | J | = 1 , the firm implements preventive maintenance for the region with the highest failure rate. In contrast, for | J | = 2 , the firm opts for a shorter repair period by opening a dedicated route for four high-failure-rate regions, none of which undergo preventive maintenance. Furthermore, compared to the single repair period, the dedicated repair period proves more effective in reducing the inequality of penalty costs across service regions, as evidenced by the uniformity of the radii of the red circles in Figure 8(b) compared to those in Figure 8(a). The corresponding Gini coefficient, depicted in Figure 7(b), highlights this disparity. Notably, the dedicated repair period exhibits a Gini coefficient of 0.43, indicating greater equality than the single repair period with a coefficient of 0.69. This insight underscores the efficacy of a dedicated repair period over preventive maintenance in addressing the challenges posed by the black hole phenomenon.

7. Conclusion

This article addresses a fundamental operational challenge in the bike-sharing industry: reducing maintenance costs through timely and cost-effective repair planning. We propose a data-driven periodic repair framework that integrates uncertain failure rates and covariate information into the repair operations of shared bikes. The key contribution of our work lies in the conceptualization of a physical landscape analogous to the black hole phenomenon in cosmology, where certain locations exhibit morbidly high failure rates. To address the undesirable black hole phenomenon, we introduce two special strategies: preventive maintenance and a dedicated repair period. Analytical results derived from a stylized two-region system demonstrate the theoretical value of these strategies. By leveraging black hole information, we formulate a novel periodic repair problem for shared bikes and address it within the operational data analytics framework, which refines data-to-decision mappings, including those derived from scenario-wise distributionally robust optimization models (DRO) that admit exact linearization under linear decision rules. Through extensive numerical studies, we validate the effectiveness of the periodic repair model and demonstrate the valuable insights gained from considering the black hole phenomenon.

While our approach leverages regional failure rates to mitigate the adverse effects of the black hole phenomenon, it does so in a continuous and stochastic manner without explicitly delineating these regions. This modeling choice provides greater flexibility in capturing spatial heterogeneity without resorting to rigidly binary region classifications. Nonetheless, explicitly identifying high-risk regions may offer additional insights and practical value. A key reason for not explicitly modeling black hole regions lies in several nontrivial challenges. First, the identification of such regions is inherently data-driven and depends on spatial aggregation choices, which may lead to instability in the resulting optimization decisions. Second, explicitly representing black hole regions typically requires binary decision structures, which would significantly increase computational complexity under the DRO framework. Third, black hole effects are intrinsically dynamic and may evolve over time, making static region-based representations potentially inadequate.

To partially address the first challenge, we introduce the concept of malfunction gathering regions (MGRs) in Appendix EC.6 as a complementary diagnostic tool. MGRs are defined as spatial clusters characterized by high densities of historical malfunctions, constructed using a combination of Kernel Density Estimation and a convex hull approach. This provides an interpretable partition of the operational area into malfunction-prone and normal regions, which can help localize persistently high-risk areas beyond the continuous modeling framework. Building on this, future research can further address these challenges in two directions. First, to address the challenges of identification and tractability, one may develop robust methods for identifying black hole regions under uncertainty and incorporate them into the DRO framework via region-dependent ambiguity sets. Second, to address the dynamic nature of black hole effects, one may extend the model to capture their temporal evolution, such as through time-dependent clustering.

Supplemental Material

sj-pdf-1-pao-10.1177_10591478261451546 - Supplemental material for Mitigating the “black holes”: Periodic repair and maintenance problem of shared bikes

Supplemental material, sj-pdf-1-pao-10.1177_10591478261451546 for Mitigating the “black holes”: Periodic repair and maintenance problem of shared bikes by Chengcheng Yu, Lan Lu, Lindong Liu and Qiao-Chu He in Production and Operations Management