The VRP, originating from the seminal work of Dantzig and Ramser [5], seeks to determine optimal routes for a fleet of vehicles serving geographically dispersed customers from one or more depots. The VRP is classified as NP-hard [15], and its practical importance has spawned an extensive taxonomy of variants. Solomon [6] introduced time window constraints, yielding the VRP with Time Windows (VRPTW), and generated benchmark instances that remain standard references nearly four decades later. The Multi-Depot VRP (MDVRP), studied by Cordeau et al. [7] and later by Montoya-Torres et al. [16], extends the problem to multiple origin points, introducing depot assignment as an additional decision dimension. Heterogeneous fleet variants [8] capture realistic fleet compositions with varying capacities, speeds, and operating costs. Recent reviews by Braekers et al. [17] and Elshaer and Awad [18] provide comprehensive classifications of VRP variants and solution methodologies. The rich VRP framework proposed by Caceres-Cruz et al. [19] emphasizes that real-world routing problems typically involve multiple complicating features simultaneously time windows, heterogeneous fleets, multiple depots, and dynamic information yet most solution methods address these features in isolation. Our work contributes to this stream by incorporating disruption events as a first-class modeling construct within the multi-depot, time-windowed routing context.

The study of supply chain disruptions has expanded substantially since the early 2000s, driven by high-profile events that exposed fragilities in global logistics networks. Snyder et al. [4] provided a comprehensive review of OR/MS models for supply chain disruptions, categorizing nearly 150 scholarly works into strategic decisions, sourcing, contracts, inventory, and facility location. They identified a persistent gap in operational-level disruption management, particularly in routing and scheduling contexts. Tomlin [20] analyzed dual-sourcing strategies under disruption risk, demonstrating that tailored mitigation approaches outperform generic contingency planning. Sheffi and Rice [12] articulated the concept of supply chain resilience as encompassing redundancy, flexibility, and velocity of recovery, laying conceptual foundations that subsequent quantitative work has sought to operationalize.

Ivanov [2] introduced the concept of “viability” in supply chain management, arguing that disruptions induced by pandemics and ripple effects require models that go beyond traditional resilience to consider survival and adaptation. Hosseini et al. [3] proposed a Bayesian network approach for quantifying supply chain resilience, incorporating disruption probability, vulnerability, and recoverability as interconnected factors. While these contributions provide valuable strategic and analytical frameworks, the translation into routing-level algorithms that can execute in real-time operational settings remains limited. Our work addresses this gap by embedding disruption modeling directly within the routing optimization algorithm.

Metaheuristic methods dominate the practical solution landscape for VRP variants due to the computational intractability of exact approaches for realistic instance sizes. Gendreau et al. [21] provided an early taxonomy of metaheuristics for VRP, covering tabu search, simulated annealing, and Genetic Algorithms (GA). Tabu search, applied to VRP by Gendreau et al. [22] and further refined by Cordeau et al. [23] for multi-depot variants, has produced many best-known solutions on benchmark instances. GAs have been successfully applied to VRPTW by Berger and Barkaoui [24] and to MDVRP by Vidal et al. [25], whose Unified Hybrid Genetic Search (UHGS) represents a landmark contribution in the field.

ALNS, proposed by Ropke and Pisinger [9], has emerged as one of the most versatile and effective frameworks for VRP. ALNS maintains a portfolio of destroy and repair operators, adaptively selecting among them based on historical performance using a roulette-wheel mechanism. Pisinger and Ropke [26] demonstrated ALNS on five VRP variants simultaneously, establishing its generalizability. Subsequent extensions include the work of Demir et al. [27] on the pollution-routing problem, Ribeiro and Laporte [28] on the cumulative VRP, and Grangier et al. [29] on the electric VRPTW and nonlinear charging. The flexibility of the ALNS framework makes it particularly suitable for incorporating disruption-specific operators, a strategy we exploit in DRSO.

Particle Swarm Optimization (PSO), introduced by Kennedy and Eberhart [30], has been adapted for discrete VRP variants through various encoding schemes [31]. The GWO by Mirjalili et al. [14] is a more recent swarm intelligence algorithm that has shown competitive performance on continuous optimization benchmarks. Extensions to discrete optimization, including GWO for TSP [32] and hybrid GWO for scheduling [33], suggest potential for VRP applications, though limited work has explored this direction. Our hybridization of GWO with ALNS is, to our knowledge, the first such combination proposed for any VRP variant.

Dynamic VRP (DVRP) settings involve information that is revealed or changes during route execution. Pillac et al. [13] provided a comprehensive classification of dynamic and stochastic VRP variants, distinguishing between problems with dynamic customer requests, dynamic travel times, and dynamic vehicle availability. Psaraftis et al. [34] surveyed dynamic vehicle routing with an emphasis on real-time information and the role of communication technologies. Ulmer et al. [35] addressed anticipatory routing in the DVRP, demonstrating that policies incorporating future information significantly outperform myopic re-optimization.

Stochastic VRP (SVRP) models uncertainty through probability distributions over problem parameters. Gendreau et al. [10] introduced the stochastic VRP with random demands, solved via an integer L-shaped method. Laporte et al. [36] provided exact and heuristic solutions for SVRP variants with stochastic demands and customers. More recently, scenario-based approaches have gained traction: Bent and Van Hentenryck [37] proposed a scenario-based online stochastic VRP, while Saint-Guillain et al. [38] developed a sample-average approximation method for the dynamic stochastic VRP.

Despite substantial progress in both dynamic and stochastic VRP research, few studies simultaneously address dynamic information revelation and disruption events within a unified framework. Notable exceptions include the work of Li et al. [39] on the VRP with disruptions and the contribution of Mu et al. [40] on disruption management for the VRPTW. However, these works do not incorporate multi-depot settings, heterogeneous disruption types, or resilience-oriented performance metrics. Our DMDVRP-D formulation and DRSO algorithm fill this gap by providing a comprehensive framework that addresses multiple disruption types, multiple depots, stochastic scenario modeling, and resilience evaluation within a single algorithmic architecture.

We formulate the DMDVRP-D as a MILP. Let G = (N, A) be a directed graph where N = D ∪ C is the set of nodes comprising depots D = {d1, d2, …, d|D|} and customers C = {c1, c2, …, c|C|}, and A = {(i, j): i, j ∈ N, i ≠ j} is the set of arcs. Let K = {k1, k2, …, k|K|} denote the set of vehicles, and T = {t0, t1, …, t|T|} denote discrete time periods within the planning horizon.

We define the following parameters:

I. qi: demand of customer i ∈ C.

II. Qk: capacity of vehicle k ∈ K.

III. cij: travel cost on arc (i, j) ∈ A.

IV. τij: travel time on arc (i, j) ∈ A.

V. [ei, li]: time window for node i ∈ N.

VI. si: service time at node i ∈ C.

VII. δRij,t ∈ {0, 1}: road closure indicator for arc (i, j) at time t.

VIII. δFd,t ∈ {0, 1}: facility failure indicator for depot d at time t.

IX. δDi,t ≥ 0: demand surge multiplier for customer i at time t.

X. δVk,t ∈ {0, 1}: vehicle breakdown indicator for vehicle k at time t.

Decision variables:

I. xijk ∈ {0, 1}: equals 1 if vehicle k traverses arc (i, j), 0 otherwise.

II. yik ∈ {0, 1}: equals 1 if customer i is assigned to vehicle k, 0 otherwise.

III. wi ≥ 0: arrival time of the vehicle at node i.

IV. ui ∈ {0, 1}: equals 1 if customer i is unserved, 0 otherwise.

V. zid ∈ {0, 1}: equals 1 if customer i is reassigned to backup depot d.

VI. fi ≥ 0: fraction of demand for customer i fulfilled through DS.

The objective function minimizes a weighted combination of routing cost, penalty for unserved customers, and disruption recovery cost:

Minimize Z = α ∑k∈K ∑(i,j)∈A cij xijk + β ∑i∈C pi ui + γ ∑i∈C ∑d∈D rid zid, (1)

where pi is the penalty for leaving customer i unserved, rid is the reassignment cost for redirecting customer i to backup depot d, and α, β, γ are weighting parameters reflecting the decision maker’s priorities.

Subject to the following constraints:

∑k∈K yik + ui = 1, ∀ i ∈ C, (2)

∑ q (1 + δD ) y ≤ Q , ∀ k ∈ K, (3)

i∈C i i,t ik k

∑j∈N xijk = yik, ∀ i ∈ C, ∀ k ∈ K, (4)

∑j∈N xijk = ∑j∈N xjik, ∀ i ∈ N, ∀ k ∈ K, (5)

x ≤ 1 − δR , ∀ (i,j) ∈ A, ∀ k ∈ K, ∀ t ∈ T, (6)

ijk ij,t ∑ x ≤ 1 − δF , ∀ d ∈ D, ∀ k ∈ K, ∀ t ∈ T, (7)

j∈N djk d,t

∑ x ≤ M (1 − δV ), ∀ k ∈ K, ∀ t ∈ T, (8)

(i,j)∈A ijk k,t

wi + si + τij − M(1 − xijk) ≤ wj, ∀ (i,j) ∈ A, ∀ k ∈ K, (9)

ei ≤ wi ≤ li, ∀ i ∈ N. (10)

Constraints (2) ensure each customer is either served by exactly one vehicle or marked as unserved. Constraints (3) enforce capacity limits adjusted for demand surges. Constraints (4) and (5) are assignment and flow conservation constraints. Constraint (6) prohibits traversal of arcs affected by road closures. Constraint (7) prevents route origination from failed depots. Constraint (8) deactivates broken-down vehicles. Constraints (9) and (10) enforce time window feasibility via big-M linearization.

The DMDVRP-D is NP-hard, as it generalizes the classical VRP. Moreover, the dynamic nature of disruption parameters and the stochastic uncertainty in their occurrence render exact solution intractable for realistic instance sizes, motivating the metaheuristic approach developed in Section 3.2.

The DRSO combines the neighborhood-based intensification capabilities of ALNS with the population-based exploration strengths of GWO. The algorithm operates in two interleaved phases: an ALNS phase that generates and refines candidate solutions through destroy-repair cycles, and a GWO phase that guides population-level search toward promising regions of the solution space.

To account for uncertainty in disruption occurrence, DRSO employs scenario-tree-based stochastic programming within its fitness evaluation. A scenario tree Τ is constructed with branching points at discrete time epochs t1, t2, …, tH within the planning horizon. Each branch represents a possible disruption realization characterized by its type (road closure, facility failure, demand surge, vehicle breakdown), location, severity, and duration.

Let Ω = {ω1, ω2, …, ω|Ω|} denote the set of scenarios generated from the scenario tree, with associated probabilities pω satisfying ∑ω∈Ω pω = 1. The stochastic fitness of a solution S is:

F(S) = (1 − η) Z(S, Δ0) + η ∑ω∈Ω pω Z(S, Δω), (14)

where Z(S, Δ0) is the objective value under the nominal (disruption-free) scenario, Z(S, Δω) is the objective value under disruption scenario ω, and η ∈ [0, 1] controls the trade-off between nominal optimality and disruption preparedness. Larger values of η produce solutions that are more resilient but potentially suboptimal under normal conditions.

Disruption probabilities are estimated from historical data and expert judgment. For the Nigerian case study, disruption probabilities are calibrated using five-year records of road closures (rainy season flooding, security incidents), facility outages (power failures, regulatory shutdowns), and demand surges (disease outbreaks, seasonal spikes). For the Chinese case study, probabilities are derived from logistics disruption databases maintained by the Shanghai Maritime University and publicly available incident reports from the YRD transportation authority.

When a disruption event occurs during route execution, DRSO triggers a re-optimization procedure that exploits the structure of the pre-disruption solution to accelerate convergence. The re-optimization strategy operates as follows:

Disruption detection: the system monitors real-time data feeds for disruption signals. Upon detection, the current disruption state vector Δ(t) is updated.

Solution partitioning: the pre-disruption solution Spre is partitioned into an unaffected component SU (routes not impacted by the disruption) and an affected component SA (routes requiring modification).

Warm-start initialization: the GWO population is initialized with the unaffected component SU fixed, and the affected component SA is replaced with solutions generated by the resilience-specific destroy-repair operators. This warm-start reduces the search space substantially compared to re-solving from scratch. Focused ALNS: the ALNS component is restricted to operate only on the affected component, with resilience-specific operators receiving elevated initial weights proportional to the disruption type (e.g., FS receives higher weight for depot failures).

Convergence and deployment: the re-optimization runs for a reduced iteration budget Nreopt < Nmax, reflecting the real-time nature of the response. The best solution found is deployed, and its structure is cached for potential subsequent re-optimization.

The complete DRSO algorithm is presented in Algorithm 1.

Algorithm 1. DRSO.

Input: Graph G=(N, A), vehicles K, disruption scenarios Ω, parameters

Output: Best solution Sbest, Composite Resilience Score CRS

1: Initialize GWO population P = {S1, S2, ..., SNw} // nearest-neighbor + random

2: Initialize ALNS operator weights wω = 1 for all operators ω 3: Evaluate F(Si) for all Si ∈ P using Eq. (14)

4: Rank wolves: identify Sα, Sβ, Sδ

5: Sbest ← Sα

6: for iter = 1 to Nmax do

7: if DisruptionDetected(t) then

8: Update disruption state Δ(t)

9: Partition Sbest into SU and SA

10: WarmStartPopulation(P, SU, SA)

11: BoostResilienceOperators(Δ(t))

12: end if

13: Update convergence parameter: a = 2(1 - iter/Nmax)

14: for each wolf Si ∈ P do

15: // GWO position update

16: Si′ ← CrossoverSelect(Sα, Sβ, Sδ, a)

17: // ALNS refinement

18: for j = 1 to NALNS do

19: Select destroy operator d by roulette wheel on weights

20: Select repair operator d by roulette wheel on weights

21: Stemp ← d (d (Si′ ))

22: if Accept(Stemp, Si′ , TSA) then 23: Si′ ← Stemp

24: end if

25: UpdateOperatorScores(πω)

26: end for

27: if F(Si′ ) < F(Si) then Si ← Si′

28: end for

29: Rank wolves: update Sα, Sβ, Sδ

30: if F(Sα) < F(Sbest) then Sbest ← Sα

31: UpdateOperatorWeights(wω) using Eq. (11) every ψ iterations

32: Cool SA temperature: TSA ← φ · TSA

33: end for

34: Compute CRS(Sbest) using Eq. (15)–(18)

35: return Sbest, CRS

To enable comprehensive evaluation of solution quality under disruptions, we introduce the CRS, a multi-criteria metric that integrates three dimensions of supply chain resilience.

RT: the time elapsed between disruption occurrence and the point at which the re-optimized solution achieves a SL within ε of the pre-disruption level:

RT = min{t′ − tdis : SL(S(t′)) ≥ (1 − ε) · SL(Spre)}, (15)

where tdis is the disruption time, S(t′) is the solution at time t′, and SL(·) is the SL function.

SL: the proportion of total demand satisfied within time windows after re-optimization:

SL = ∑ q f (1 − u ) / ∑ q (1 + δD ). (16)

i∈C i i i i∈C i i,t

EC: the additional cost incurred due to disruption response, normalized by the pre-disruption cost:

EC = (Z(Spost, Δ(t)) − Z(Spre, Δ0)) / Z(Spre, Δ0). (17)

The CRS integrates these dimensions using a weighted geometric mean:

CRS = (RT )w1 · (SL)w2 · (1 − EC )w3, (18)

norm norm

where RTnorm and ECnorm are min-max normalized values, and w1 + w2 + w3 = 1. Higher CRS values indicate greater resilience. Default weights are set as w1 = 0.3, w2 = 0.4, w3 = 0.3, reflecting the primacy of SL in supply chain performance.

We evaluate DRSO on two categories of test instances: modified Solomon VRPTW benchmarks and real-world case studies.

DRSO is compared against five baseline approaches:

ALNS-only: standard ALNS with the same operator portfolio but without GWO-guided search or resilience-specific operators.

GA: GA with order crossover, mutation, and elitism, calibrated following Vidal et al. [25].

PSO: discrete PSO with swap-sequence encoding, following Ai and Kachitvichyanukul [31].

GWO-only: standalone GWO adapted for VRP with the discrete position update of Section 3.2.3 but without ALNS refinement.

Gurobi (TL): the MILP formulation (Section 3.1) solved using Gurobi 10.0 with a time limit equal to the wall-clock time consumed by DRSO on the same instance.

DRSO parameters are calibrated through preliminary experiments on a subset of 10 Solomon-derived instances not included in the final test set. Table 4 summarizes the parameter settings.

All experiments are conducted on a workstation equipped with an Intel Xeon E5-2690 v4 processor (2.60 GHz, 14 cores) and 128 GB RAM, running Ubuntu 22.04 LTS. DRSO and all metaheuristic baselines are implemented in Python 3.11 with NumPy 1.26 and compiled Cython extensions for performance-critical operations. Gurobi 10.0.3 is used for MILP solving with default settings except for the imposed time limit. Each metaheuristic algorithm is executed 30 independent times per instance with different random seeds. Statistical significance of pairwise differences is assessed using the Wilcoxon signed-rank test at a significance level of p < 0.05. Effect sizes are reported using the Vargha-Delaney A-measure [42], where A > 0.71 indicates a large effect. We also report the Interquartile Range (IQR) of results to characterize solution consistency.

Table 5 presents the aggregated results across all 90 modified Solomon configurations, reporting mean values over 30 independent runs for each algorithm. DRSO consistently outperforms all baseline methods across the three primary metrics: disruption recovery cost, RT, and percentage of unserved customers.

DRSO achieves mean improvements of 20.9% over ALNS-only, 27.2% over GA, 30.0% over PSO, 25.4% over GWO-only, and 17.5% over Gurobi (TL) in recovery cost. RT improvements range from 23.5% (vs. Gurobi) to 41.0% (vs. PSO). The reduction in unserved customers is particularly notable, ranging from 26.9% to 46.5% across baselines. All pairwise differences are statistically significant at p < 0.05 via the Wilcoxon signed-rank test.

Table 6 presents results disaggregated by disruption severity level, revealing that DRSO’s advantage increases with disruption severity.

Under high disruption severity, DRSO achieves 28.1% cost reduction and 41.5% faster recovery compared to ALNS-only. This escalating advantage under severe conditions confirms the effectiveness of the scenario-tree stochastic evaluation and resilience-specific operators in anticipating and adapting to complex disruption patterns.

Table 7 presents results on the Nigerian pharmaceutical distribution network across three representative disruption profiles: rainy-season flooding (Profile A), combined flooding and security disruptions (Profile B), and compound disruptions including demand surge from disease outbreak (Profile C).

DRSO demonstrates substantial advantages across all profiles, with cost reductions of 19.2% (Profile A), 23.1% (Profile B), and 25.0% (Profile C) compared to ALNS-only. Under the most challenging compound disruption, Profile C, DRSO maintains a 96.4% average SL, compared to 89.7% for ALNS-only and 84.9% for PSO. The RT advantage (38.2 minutes vs. 58.6–68.1 minutes for baselines) is critical in the pharmaceutical context, where delayed deliveries of temperature-sensitive medications can result in spoilage and patient harm.

Table 8 summarizes results on the YRD manufacturing network across four disruption profiles: typhoon-induced closures (Profile D), river flooding with modal disruption (Profile E), factory shutdowns with demand redistribution (Profile F), and compound multi-type disruptions (Profile G).

On the larger Chinese network (78 nodes), DRSO’s advantages are even more pronounced, with cost reductions of 19.1% to 23.2% compared to ALNS-only across disruption profiles. The MS operator proves particularly effective in the YRD context, where inland waterway alternatives along the Yangtze River provide viable backup routes when road or rail links are disrupted. DRSO achieves a 97.1% average SL, critical for just-in-time automotive manufacturing where line stoppages due to missing components incur costs of approximately ¥500,000–2,000,000 per hour.

Table 9 presents the Wilcoxon signed-rank test results and Vargha-Delaney A-measure effect sizes for pairwise comparisons between DRSO and each baseline, aggregated across all test instances.

All comparisons yield p < 0.05, confirming statistical significance. The A-measure values (0.74–0.93) consistently indicate large effect sizes, demonstrating that DRSO’s superiority is not merely statistically significant but practically meaningful. The smallest effect size is observed against Gurobi (TL), which benefits from the optimality guarantees of exact methods for the portion of the solution it can compute within the time limit.

Table 10 reports the average usage frequency and contribution of each ALNS operator across all DRSO runs, providing insight into the algorithm’s adaptive behavior.

The resilience-specific operators consistently achieve higher final weights than standard operators, indicating their effectiveness in the disruption context. FS achieves the highest contribution rate (34.8% of uses leading to improvement), followed by ER (32.1%). These results confirm that domain-specific destroy-repair operations provide substantially greater search guidance than generic operators when the problem involves disruption-induced structural changes to the solution.

Table 11 reports average wall-clock computation times across instance categories.

DRSO requires approximately 1.8–2.4 times the computation time of standalone methods for initial solution construction due to the overhead of GWO population management and scenario evaluation. However, the warm-start re-optimization mechanism delivers a critical advantage in dynamic settings: DRSO’s re-optimization time (18.4 seconds average) is 44–62% faster than baseline metaheuristics, which must re-solve the disrupted problem without structural guidance from the pre-disruption solution. Re-optimization speed is operationally significant, as re-optimization speed directly determines the delay between disruption detection and route adjustment deployment.

Table 12 presents an ablation study isolating the contribution of each DRSO component. Five variants are tested: DRSO without GWO (ALNS with resilience operators only), DRSO without resilience operators (ALNS+GWO with standard operators only), DRSO without scenario-tree evaluation (deterministic fitness), DRSO without warm-start re-optimization (cold restart), and the full DRSO.

The ablation study reveals that resilience-specific operators contribute the largest individual improvement (19.8% cost gap when removed), followed by scenario-tree evaluation (14.8%), GWO integration (12.3%), and warm-start re-optimization (8.3%). The warm-start mechanism, while having the smallest cost impact, has a disproportionately large effect on RT (18.2 TU without vs. 12.4 TU with), confirming its critical role in dynamic settings where rapid response is essential.

The computational results provide strong evidence for the effectiveness of the DRSO framework across diverse problem settings. Several factors contribute to DRSO’s superior performance. First, the hybridization of ALNS and GWO creates a synergistic search mechanism where GWO’s population-based exploration prevents premature convergence, a known weakness of standalone ALNS [26], while ALNS’s destroy-repair cycles provide the fine-grained neighborhood exploitation that pure population-based methods lack. The ablation study confirms this complementarity: removing either component degrades performance, but removing both would reduce DRSO to a basic heuristic.

Second, the resilience-specific operators encode domain knowledge about supply chain disruption response strategies into the search mechanism. Unlike generic VRP operators that treat solution modification as purely topological (rearranging visit sequences), the FS, ER, DS, and MS operators operate on the structural level of supply chain configuration. This semantic richness accelerates the discovery of high-quality solutions in disrupted scenarios, as evidenced by the operators’ consistently higher adaptive weights and improvement rates compared to standard operators (Table 10).

Third, the scenario-tree stochastic fitness evaluation enables DRSO to construct solutions that inherently hedge against potential disruptions, even before they occur. This proactive resilience contrasts with purely reactive approaches that optimize for the nominal scenario and rely on re-optimization after disruptions. The value of anticipatory optimization is particularly evident under high-severity disruption conditions (Table 6), where DRSO’s advantage over deterministic approaches is most pronounced.

The two real-world case studies demonstrate DRSO’s applicability to distinct supply chain contexts with fundamentally different disruption profiles. The Nigerian pharmaceutical case study highlights the algorithm’s value in infrastructure-constrained environments where disruptions are frequent, diverse, and often compound (e.g., simultaneous flooding and security incidents). The ability to maintain 96.4% SLs under compound disruptions is particularly significant for pharmaceutical logistics, where delivery failures can have direct health consequences. The 38.2-minute average re-optimization time is well within the operational decision window for pharmaceutical delivery managers who typically have 1–2 hours between disruption detection and route commitment.

The Chinese manufacturing case study illustrates DRSO’s scalability to larger, more complex networks with multi-modal transportation options. The MS operator’s effectiveness in the YRD context, where Yangtze River waterways provide natural backup for disrupted road and rail links, demonstrates the importance of embedding mode-specific knowledge into the optimization framework. The 97.1% SL and 22.7-minute average re-optimization time meet the stringent requirements of JIT manufacturing supply chains where component delays translate directly to production line stoppages.

From a managerial perspective, DRSO’s CRS metric provides decision-makers with a single, interpretable measure of solution resilience that balances the competing objectives of rapid recovery, service continuity, and cost control. The configurable weights (w1, w2, w3) allow organizations to tailor the metric to their specific priorities; for example, pharmaceutical distributors may weight SL more heavily, while cost-sensitive manufacturers may emphasize EC minimization.

DRSO’s consistent superiority over Gurobi (TL) merits discussion. Gurobi, as a state-of-the-art commercial solver, provides optimality guarantees for the portion of the MILP it can solve within the time limit. However, the dynamic, stochastic, and multi-objective nature of the DMDVRP-D challenges exact methods in several ways: 1) the scenario-tree evaluation increases the effective problem size by a factor of |Ω|, 2) the dynamic re-optimization requirement demands rapid solution generation that exceeds Gurobi’s warm-start capabilities for modified MILPs, and 3) the non-linear resilience metric (CRS) cannot be directly optimized within the MILP framework. These observations align with the broader finding in the metaheuristic literature that well-designed heuristics outperform exact methods under time constraints for large-scale, structurally complex optimization problems [43].

Several limitations of this study should be acknowledged. First, the disruption scenarios, while calibrated from real data, are generated synthetically and may not capture the full complexity of real-world disruption dynamics, including cascading effects and interdependencies between disruption types. Second, the scenario-tree approach with |Ω| = 20 scenarios represents a trade-off between computational tractability and disruption coverage; larger scenario sets may improve hedging quality at the cost of increased computation time. Third, the real-world case studies, while based on authentic network topologies and disruption patterns, use partially simulated demand data due to confidentiality constraints. Fourth, DRSO’s computational overhead compared to standalone methods (Table 11) may limit its applicability in extremely time-critical settings where sub-second re-optimization is required. Fifth, the current formulation assumes complete and accurate disruption detection; in practice, disruption information may be partial, delayed, or noisy, requiring robust optimization extensions.

Internal validity is supported by the use of 30 independent runs with controlled random seeds and non-parametric statistical testing. However, parameter sensitivity analysis, while conducted informally during calibration, is not reported exhaustively. External validity is enhanced by testing on both standardized benchmarks and real-world instances from two continents and two industry sectors, though generalizability to other supply chain types (e.g., cold chain, humanitarian logistics) remains to be demonstrated. Construct validity of the CRS metric depends on the appropriateness of the weight settings; while default weights are justified by supply chain management literature, alternative weighting schemes may be more suitable for specific applications.