A Branch-and-Price Approach to the Platform Supply Vessel Routing and Scheduling Problem with Uncertain Demand

1. Introduction

Against the backdrop of global oil and gas exploration and development expanding into deep and ultra-deep waters, supply vessels serve as the core logistics link connecting onshore supply bases with offshore facilities such as drilling platforms and production platforms [1]. Their routing and scheduling optimization directly determines the continuity and economic efficiency of offshore oil and gas operations. Given the limited space on offshore installations, regular deliveries of equipment and materials, as well as the removal of waste materials, rely heavily on Platform Supply Vessel (PSV) services. As one of the costliest resources in offshore logistics, PSVs require high operational efficiency and rational scheduling for effective cost control in oil and gas companies. Therefore, the Platform Supply Vessel Routing and Scheduling Problem (PSVRSP) has become a central topic in upstream supply chain optimization for the offshore industry [2].

The PSVRSP is essentially a complex combinatorial optimization problem, whose difficulty arises from the interplay of multi-dimensional constraints and various sources of uncertainty. In terms of constraints, scheduling must comply with visit-frequency requirements of offshore platforms, PSV capacity limits, onshore base service capacity (maximum daily departures), PSV-specific characteristics (power, speed, cargo-capacity heterogeneity), and time-window restrictions related to operational windows and voyage duration limits. In terms of uncertainty, volatile demand from offshore platforms and disruptions caused by adverse weather conditions often invalidate pre-planned schedules, leading to extra costs such as vessel chartering and route rerouting. Moreover, the PSVRSP is inherently NP-hard. As the number of offshore facilities increases, the solution space expands exponentially, rendering traditional experience-based scheduling and basic optimization methods inadequate for real-world applications.

Existing studies have achieved considerable progress in the optimization of PSVRSP. Early research mainly focused on deterministic environments, where scheduling models were typically formulated using Mixed-Integer Linear Programming (MILP) and solved with heuristic algorithms such as genetic algorithms and tabu search, providing preliminary optimization for fleet composition and route planning. With a deeper understanding of practical operational scenarios, the impacts of demand uncertainty have attracted increasing attention. Researchers have adopted approaches including two-stage stochastic programming, robust optimization, and simulation–optimization embedded with metaheuristics to incorporate demand volatility, thereby improving the practicality and robustness of scheduling solutions. Nevertheless, the application of exact solution methods in PSVRSP remains relatively limited. A notable contribution was made by Friedberg [3], who introduced a Branch-and-Price (B&P) algorithm for an MILP model of PSVRSP, in which supply vessels operated in coordination with bunkering vessels as offshore hubs. Their work validated the effectiveness of decomposition-based exact methods for PSV routing problems. However, their model was developed under a deterministic setting and did not account for demand uncertainty. To date, few exact algorithms have been proposed for PSVRSP under stochastic demand, leaving a clear methodological gap in the literature.

Despite the advances covering both deterministic and stochastic settings as well as various heuristic approaches, critical research gaps still exist. Most studies addressing demand uncertainty rely on heuristics or metaheuristics, which cannot guarantee solution optimality. Meanwhile, the few studies that employ exact algorithms ignore demand variability and thus deviate from real-world operational conditions. Therefore, developing an exact solution framework that integrates demand uncertainty and enables rigorous global optimization is of great theoretical significance and practical value for enhancing the reliability and cost-effectiveness of offshore logistics. Building on the existing literature, this study focuses on the complexity and practical requirements of the PSVRSP. By incorporating multiple operational constraints and demand uncertainty, this paper aims to develop more efficient exact optimization approaches, so as to provide scientific and practical decision support for offshore oil and gas companies and facilitate the refined management and sustainable development of offshore logistics.

2. Literature Review

This section provides a systematic review of the literature related to the Platform Supply Vessel Routing and Scheduling Problem (PSVRSP) in both deterministic and uncertain environments.

Early research on PSVRSP can be traced back to Fagerholt [4], who first established a Mixed-Integer Linear Programming (MILP) model for the periodic routing and scheduling of heterogeneous platform supply vessels, considering pickup and delivery tasks for deck cargo. This work laid the foundation for subsequent modeling of supply vessel operations. Aas et al. [5] further expanded this model by incorporating constraints related to deck-space availability on offshore platforms. For large-scale instances, Gribkovskaia et al. [6] designed a tabu search heuristic to improve computational efficiency. Shyshou et al. [7] focused on simplified versions of the problem and developed customized heuristics. Halvorsen-Weare et al. [8] formulated an MILP model that limited voyage duration and controlled departure times from the supply base to improve schedule feasibility. Although the existing literature has increasingly addressed multi-commodity, multi-vessel, and complex operational settings, most early studies relied on heuristic algorithms to obtain feasible solutions. Friedberg et al. [3] introduced an important advancement by applying a Branch-and-Price (B&P) algorithm to an MILP model for platform supply vessel routing and scheduling, in which supply vessels operated alongside refueling vessels as offshore hubs. This study demonstrated the strength of decomposition techniques for platform supply vessel routing problems.

Research on the Supply Vessel Planning Problem (SVPP) under uncertainty includes papers addressing uncertain weather conditions, uncertain demand, and the combination of both uncertain demand and uncertain weather conditions. Most studies adopt a two-stage framework that integrates simulation and optimization and can be divided into two categories. The first category constructs robust schedules using deterministic models enhanced by either probabilistic constraints to accommodate uncertain demand [9] or time slack to absorb delays caused by rough weather [8]. The second category estimates the average undelivered demand via simulation and then incorporates penalty costs into a deterministic optimization model [10]. Following a similar philosophy, Cruz et al. [11] generate candidate routes under stochastic demand and travel times and then use reliability levels to build robust weekly schedules. Meanwhile, Ksciuk et al. [12] provides a comprehensive review of uncertainty handling in maritime routing, including supply vessel scheduling.

Adverse weather is widely recognized as the primary cause of schedule failures [13]. To improve robustness, Halvorsen-Weare et al. [8] propose a set-covering model with mandatory time slack between consecutive voyages. However, they do not verify actual robustness through simulation. Halvorsen-Weare et al. [10] later extend this work by using discrete-event simulation to compute expected unmet demand and assign penalties accordingly. For uncertain demand, common strategies include enforcing voyage-level chance constraints to guarantee service reliability or inflating demand forecasts to accommodate variability. Notably, Halvorsen-Weare et al. combine demand inflation and time buffers to address both demand and weather uncertainty simultaneously.

Nevertheless, these approaches have notable drawbacks. Most of these methods do not aim to minimize the actual expected total cost, which can result in economically inefficient schedules. By comparison, Monte Carlo Sampling has proven effective for discrete stochastic optimization problems where full scenario enumeration is intractable, including stochastic vehicle routing and supply chain design [14]. It is also widely used to evaluate solutions within heuristic frameworks to handle both NP-hardness and scenario intractability [15].

Despite these advances, most studies addressing the stochastic Supply Vessel Routing and Scheduling Problem (SVRSP) rely on heuristics or simulation-based methods. Few employ exact solution methods such as Branch-and-Price, and even fewer integrate such methods with Monte Carlo sampling to address stochastic demand while minimizing expected total cost under time-window constraints. To fill this gap, this study proposes a Branch-and-Price algorithm embedded with Monte Carlo sampling to solve the SVRSP with uncertain demand and time windows, with the objective of minimizing the expected total cost.

3. Problem Description

A PSVRSP typically involves several core components, including an onshore supply base, a set of offshore platforms that generate orders, and a fleet of Platform Supply Vessels (PSVs), which must fulfill all service orders. In practical offshore logistics operations, each platform may issue multiple orders, and each order exhibits diverse loading and unloading requirements as well as time windows during which vessels must initiate service.

In practical offshore logistics operations, the PSVRSP can be described on a graph $\mathcal{G}(\mathcal{V},\mathcal{A})$, where the node set $\mathcal{V}$ represents order locations, and the arc set $\mathcal{A}$ represents sailing legs, each associated with the sailing distance between pairs of orders. Let $\mathcal{C}$ denote the set of offshore platforms that require service, and $\mathcal{O}$ the set of all orders to be performed. Let $r$ r be the unit cargo-handling rate at each offshore platform, and $\mathcal{K}$ the set of available PSVs. In this work, a trip is defined as a complete voyage cycle performed by a single vessel, consisting of three sequential tasks: (i) departing empty from the onshore supply base; (ii) serving orders offshore; and (iii) returning to the onshore supply base to unload any pickup cargo collected during the trip.

Each vessel $k\in \mathcal{K}$ is allowed to complete at most $L$ consecutive trip. The entire sequence of trips performed by a vessel is collectively referred to as a route. Figure 1 illustrates the structure of the transportation network.

3.1. Demand Uncertainty Modeling and Scenario Generation

In this section, we present a method to approximate the probability distribution of demand at offshore installations. We assume that demands at different installations are mutually independent; that is, demand at each installation follows an independent probability distribution, which may include normal, Erlang, or Weibull distributions [9]. Based on such probabilistic representations, several modeling approaches can be employed, including robust optimization, chance-constrained programming, and scenario-based stochastic programming. Robust and chance-constrained models provide protection against uncertainty but may lead to overly conservative decisions or complex probabilistic constraints [16]. By comparison, scenario-based stochastic programming discretizes the underlying demand distribution into a finite set of representative realizations, which supports direct comparisons across different demand outcomes. Given this balance between modeling fidelity and computational tractability, this study adopts scenario-based representations. To incorporate this uncertainty into the model, demand is treated as a random variable, and a combination of scenario generation and scenario reduction approach is employed.

To generate a diverse set of demand scenarios, this study adopts an approach that combines LHS and Cholesky decomposition [17]. LHS provides stratified sampling that ensures more uniform coverage of multidimensional distribution space. However, when sampling multiple random variables, their inherent correlation may reduce representativeness of the generated scenarios. To address this issue, Cholesky decomposition is incorporated to weaken inter-order demand correlations, thereby improving independence and diversity of the scenario set.

Suppose there are $\left|{\mathcal{N}}_{\mathcal{O}}\right|$ orders, and their demands are represented by $X_1,\dots ,X_{\left|{\mathcal{N}}_{\mathcal{O}}\right|}$. The cumulative distribution function for the demand at order n is expressed as follows,

$$P_n=F_n\left(X_n\right).$$

(1)

Using this formulation, a total of M scenarios is generated through the process. First, an ordered sampling matrix is constructed via the LHS method. Next, the Cholesky decomposition is applied to adjust the sampled vectors and mitigate correlation among port demands. The final output is an M-scenario demand matrix for the $\left|{\mathcal{N}}_{\mathcal{O}}\right|$ orders, completing the multi-scenario generation procedure.

To improve computational efficiency and reduce redundant scenarios, the stepwise backward reduction (SBR) algorithm [18] is applied. This algorithm iteratively computes the distances between scenarios and removes the scenario with the smallest distance to others until the initial scenarios are reduced to the predefined number. The final output consists of representative scenarios and their associated probabilities. Specifically, $\mathcal{S}$ denotes the scenario set $\mathcal{S}=\{\mathrm{1,2},\dots ,|\mathcal{S}|\}$, and each scenario $s\in \mathcal{S}$ occurs with probability $p_s$.

3.2. Mathematical Modeling

In this section, we present an MILP formulation for the PSVRSP-UD. Before we proceed, we first introduce some pertinent notation. Table 1 lists the sets, parameters, and variables involved in the model.

In this study, duplicated supply base nodes are used to model the multi-trip scheduling problem for a single vessel [19]. By constructing an ordered chain of depot-copy nodes to characterize the trip sequence, cross-trip service is avoided without introducing additional trip-index variables. To ensure the sequential use of duplicated supply base nodes, the arc set A is defined as $\mathcal{A}={\mathcal{A}}_{\mathcal{O}\mathcal{O}}\cup {\mathcal{A}}_{\mathcal{D}\mathcal{O}}\cup {\mathcal{A}}_{\mathcal{O}\mathcal{D}}\cup {\mathcal{A}}_{\mathcal{D}\mathcal{D}}$, where ${\mathcal{A}}_{\mathcal{O}\mathcal{O}}=\{\left(\mathrm{i},\mathrm{j}\right):\mathrm{i},\mathrm{j}\in \mathcal{O},i\ne \mathrm{j}\}$ represents arc set between customer nodes, ${\mathcal{A}}_{\mathcal{D}\mathcal{O}}=\{\left(d,\mathrm{j}\right):d\in \mathcal{D},\mathrm{j}\in \mathcal{O}\}$ represents arc set from duplicated supply base nodes to customer nodes, ${\mathcal{A}}_{\mathcal{O}\mathcal{D}}=\{\left(i,d\right):d\in \mathcal{D},i\in \mathcal{O}\}$ represents arc set from customer nodes to duplicated supply base nodes, ${\mathcal{A}}_{\mathcal{D}\mathcal{D}}=\{\left(0^l,0^{l+1}\right):l=\mathrm{1,2},\dots ,L\}$ represents sequential arc set between consecutive duplicated supply base nodes.This definition ensures that supply vessels can only visit duplicated supply base nodes in strict trip sequence within the planning horizon, and any non-sequential jump across trips is prohibited.

The PSVRSP-UD model aims to generate feasible routes and schedules for platform supply vessels (PSVs) to minimize the total cost while satisfying all order fulfillment and operational constraints. As PSVs represent one of the costliest resources in offshore logistics, and in response to the requirements of energy conservation and emission reduction in international shipping, the total cost consists of vessel operating cost and carbon emission cost, where $\alpha$ is a weight coefficient that balances the importance of the two cost components. The carbon emission cost adopts the fuel consumption function proposed by Psaraftis and Kontovas [20], which comprehensively considers the nonlinear relationship among demand scenario probability, vessel speed, payload, and fuel consumption. The complete MIP model of PSVRSP-UD is as follows.

$$\mathrm{M}\mathrm{i}\mathrm{n}\alpha {\sum }_{k\in \mathcal{K}}{C}_t(t_{0^{L+1}}^k-t_{0^1}^k)+\left(1-\alpha \right){\sum }_{s\in S}{p}_s{\sum }_{k\in K}{\sum }_{\left(i,j\right)\in A}\left[C_f{\cdot E}_k\cdot v^3\cdot {\left(B_k+l_{ik}^s\right)}^{{\displaystyle \frac{2}{3}}}\cdot {\tau }_{ij}\cdot x_{ij}^k\right],\alpha ϵ[\mathrm{0,1}],$$

(2)

Subject to

$$\sum _{\left(0^1,j\right)\in \mathcal{A}}{x}_{0^{1}j}^k\le 1\forall k\in \mathcal{K},$$

(3)

$$\sum _{\left(i,0^{L+1}\right)ϵ\mathcal{A}}{x}_{i{0}^{L+1}}^k\le 1\forall k\in \mathcal{K},$$

(4)

$$\sum _{\left(0^l,j\right)\in \mathcal{A}}{x}_{0^{l}j}^k=\sum _{\left(i,0^l\right)ϵ\mathcal{A}}{x}_{i{0}^l}^k\forall k\in \mathcal{K},l=2,\dots ,L,$$

(5)

$$\sum _{\left(i,j\right)ϵ\mathcal{A}}{x}_{ij}^k=\sum _{\left(j,i\right)ϵ\mathcal{A}}{x}_{ji}^k\forall k\in \mathcal{K},i\in \mathcal{O},$$

(6)

$$\sum _{k\in \mathcal{K}}\sum _{\left(i,j\right)ϵ\mathcal{A}}{x}_{ij}^k=1\forall k\in \mathcal{K},j\in \mathcal{O},$$

(7)

$$\sum _{\left(i,j\right)ϵ\mathcal{A}}{x}_{ij}^k\le \sum _{\left(j,i\right)ϵ\mathcal{A}}{x}_{ji}^k\forall k\in \mathcal{K},i\in \mathcal{D}/\{0^1,0^{L+1}\},$$

(8)

$$l_i^{ks}\le Q\forall k\in \mathcal{K},i\in \mathcal{V},s\in \mathcal{S},$$

(9)

$$l_i^{ks}+q_j^s\le l_j^{ks}+M(1-x_{ij}^k)\forall \left(i,j\right)\in \mathcal{A},k\in \mathcal{K},s\in \mathcal{S},$$

(10)

$$l_{0^l}^{ks}=0,\forall k\in \mathcal{K},s\in \mathcal{S},l=2,\dots ,L,$$

(11)

$${\tau }_{ij}={\displaystyle \frac{d_{ij}}{v}}\forall i,j\in \mathcal{V},$$

(12)

$$s_i={\displaystyle \frac{q_i^s}{r}}\forall i\in \mathcal{O},s\in \mathcal{S},$$

(13)

$$a_i-M\left(1-\sum _{\left(i,j\right)ϵ\mathcal{A}}{x}_{ij}^k\right)\le t_i^k\forall i\in \mathcal{V},k\in \mathcal{K},$$

(14)

$$t_i^k\le b_i+M\left(1-\sum _{\left(i,j\right)ϵ\mathcal{A}}{x}_{ij}^k\right)\forall i\in \mathcal{V},k\in \mathcal{K},$$

(15)

$$t_i^k+{\tau }_{ij}+s_j\le t_j^k+M\left(1-x_{ij}^k\right)\forall \left(i,j\right)\in \mathcal{A},k\in \mathcal{K},$$

(16)

$$x_{ij}^k\in \left\{\mathrm{0,1}\right\}\forall i,j\in \mathcal{V},k\in \mathcal{K},$$

(17)

$$l_i^{ks}\ge 0\forall i\in \mathcal{V},k\in \mathcal{K}，s\in \mathcal{S},$$

(18)

$$T\ge t_i^k\ge 0\forall i\in \mathcal{V},k\in \mathcal{K},$$

(19)

The objective function (2) is to minimize the expected value of the total cost over all demand scenarios. Constraint (3) is the out-degree constraint for the starting duplicated supply base node. Constraint (4) is the in-degree constraint for the ending duplicated supply base node. Constraint (5) is the flow conservation constraint for intermediate duplicated supply base nodes, ensuring the orderly connection between consecutive trips. Constraint (6) is the flow balance constraint for order nodes. Constraint (7) is the uniqueness constraint, ensuring each order is served exactly once. Constraint (8) ensures that each vessel must visit the duplicated supply base nodes in sequence. Constraint (9) is the maximum payload constraint, ensuring the vessel load does not exceed capacity Q under any demand scenario. Constraint (10) is the payload update constraint, ensuring the load of vessel $k$ under scenario s is properly updated when traversing arc $\left(i,j\right)$. Constraint (11) is the payload reset constraint, ensuring the vessel completes unloading at each intermediate duplicated supply base node. Constraint (12) defines the calculation of vessel sailing time. Constraint (13) defines the calculation of order service time. Constraints (14)-(15) are order service time window constraints, ensuring the service start time falls within the specified interval. Constraint (16) is the time continuity constraint, ensuring the start time of service for order j is no earlier than the completion time of service for order $i$ plus the sailing time from $i$ to $j$. Constraints (17)-(19) are variable definition constraints.

4. Proposed Branch-and-Price Algorithm for The PSVRSP-UD

To precisely solve the PSVRSP-UD, this study designs and implements a B&P algorithm [21]. Its core idea is to embed Column Generation (CG) within a Branch-and-Bound (B&B) framework [22]. By dynamically generating feasible trips with negative reduced cost, it avoids the combinatorial explosion caused by explicitly enumerating all routes, thereby achieving efficient global optimal solution of large-scale integer programming models. The following sequentially presents the model reformulation, pricing subproblem solution, branching strategy, and integer solution repair mechanism.

4.1. Model Reformulation via Dantzig-Wolfe Decomposition

The PSVRSP-UD is essentially analogous to the Multi-Trip Vehicle Routing Problem with Time Windows (MTVRPTW). For such multi-trip routing and scheduling problems, two typical set-partitioning models have been widely studied. The first model uses a complete vessel route as a column, directly representing a consecutive multi-trip sequence from departure to return [23]. Although structurally intuitive, it leads to an explosion of the column space due to the combination of trip connections and time schedules, making it impractical for real-scale instances. The second model takes a single trip as the basic column unit and does not explicitly model the full route of a vessel [24]. It controls the non-overlapping trips of the same vessel via constraints, drastically reducing problem complexity and difficulty. This study adopts the second modeling approach and follows the exact solution paradigm for multi-trip problems.

Before presenting the formal set-partitioning model, we first provide a formal definition of a trip, which serves as the foundation for the subsequent column generation and pricing problem.

Definition 1.

A trip $r$ is composed of an ordered node sequence ${cus}_r$, a planned start time $t_r^{start}$, a maximum postponable time $t_r^{relax}$ and a planned end time $t_r^{end}$ such that:

• The total loading demand on the sequence does not exceed the vessel’s capacity constraint;
• The vessel can start service at $t_r^{start}$, serve all orders according to the given sequence, return to the onshore supply base, and satisfy all time windows and sailing time constraints;
• The planned start time $t_r^{start}$ and maximum postponable time $t_r^{relax}$ jointly form the feasible departure time interval of the trip, covering all valid departure times under the visiting sequence.

The cost of trip $r$ under scenario $s(s\in \mathcal{S})$ is denoted $c_r^s$, and the set of trips is denoted $\Omega$. Based on this definition, the PSVRSP-UD can be modeled as a set-partitioning problem with mutual exclusion constraints, where $y_r$ denotes whether trip $r$ is selected, and the set to be covered is the order demand of offshore platforms. Mutual exclusion constraints ensure that selected trips can be assigned to K PSVs without overlapping in time for the same PSV. This paper enforces that at most K PSVs are used at any time. By discretizing time, a Boolean indicator $u_{tr}$ is defined to indicate whether trip $r$ covers time $t$, such that the number of occupied vessels at any time does not exceed the fleet size.

The resulting set-partitioning formulation of the PSVRSP-UD is presented as follows:

$$\mathrm{m}\mathrm{i}\mathrm{n}\sum _{r\in \Omega }\sum _{\forall s\in \mathcal{S}}{p}_s{c}_r^s{y}_r ,$$

(20)

Subject to

$$\sum _{r\in \Omega }{a}_{ir}{y}_r=1\forall i\in \mathcal{O},$$

(21)

$$\sum _{r\in \Omega }{u_{tr}y}_r\le K\forall t\in \{0,...,T\times {10}^{nb}\},$$

(22)

$$y_r\in \left\{\mathrm{0,1}\right\}\forall r\in \Omega .$$

(23)

where $a_{ir}$ indicates whether order $i$ is covered by trip $r$, $nb$ is the number of decimal digits for time precision. Constraint (21) ensures each order is served once. Constraint (22) ensures the number of PSVs used at any time intervals does not exceed the fleet size limit.

To avoid an explosion of constraints due to overly fine time precision, this study adopts a coarse time interval relaxation strategy. If the relaxed solution is feasible, it is optimal; otherwise, the granularity is refined. Specifically, the planning horizon $[0,T]$ is partitioned as follows. Let $l_{\mathrm{m}\mathrm{i}\mathrm{n}}$ be a small value guaranteeing that the duration of any feasible trip is greater than $l_{\mathrm{m}\mathrm{i}\mathrm{n}}$. The interval is partitioned as ${∆}_{\mathrm{t}}=[l_{\mathrm{m}\mathrm{i}\mathrm{n}}\times t,l_{\mathrm{m}\mathrm{i}\mathrm{n}}\times (t+1)]t\in \{0,...,⌊{\displaystyle \frac{T}{l_{\mathrm{m}\mathrm{i}\mathrm{n}}}}⌋-1\}$ ，and Constraint (22) is replaced by constraint (24) .

$$\sum _{r\in \Omega }{b_{tr}y}_r\le K\forall t\in \{0,...,⌊{\displaystyle \frac{T}{l_{\mathrm{m}\mathrm{i}\mathrm{n}}}}⌋\},$$

(24)

where $b_{tr}\in [\mathrm{0,1}]$ is the time proportion of trip $r$ within interval ${∆}_{\mathrm{t}}$. Since $l_{\mathrm{m}\mathrm{i}\mathrm{n}}$ is smaller than any trip duration, $l_{min}$ can accurately compute the start time of any trip without loss of precision. If the optimal relaxed solution is infeasible (violating (22)), then $l_{\mathrm{m}\mathrm{i}\mathrm{n}}$ is halved and recomputed.

The relaxed linear master problem (RLMP) is formulated as follows.

$$\mathrm{m}\mathrm{i}\mathrm{n}\sum _{r\in \widehat{\Omega }}\sum _{\forall s\in \mathcal{S}}{p}_s{c}_r^s{y}_r,$$

(25)

Subject to

$$\sum _{r\in \widehat{\Omega }}{a}_{ir}{y}_r\ge 1\forall i\in \mathcal{O},$$

(26)

$$\sum _{r\in \widehat{\Omega }}{b_{tr}y}_r\le K\forall t\in \{0,...,⌊{\displaystyle \frac{T}{l_{min}}}⌋\},$$

(27)

$$1\ge y_r\ge 0\forall r\in \widehat{\Omega }.$$

(28)

4.2. Pricing Subproblem: NG-Route Labeling Algorithm

The pricing subproblem seeks a trip with negative reduced cost under all scenarios. Given dual values ${\pi }_i(i\in \mathcal{O})$ and ${\psi }_t(t\in \{0,...,⌊{\displaystyle \frac{T}{l_{\mathrm{m}\mathrm{i}\mathrm{n}}}}⌋\})$corresponding to Constraints (26) and (27) of the RLMP, the reduced cost for $r(r\in \widehat{\Omega })$ is presented as follows.

$${\sum }_{\forall s\in \mathcal{S}}{p}_s{c}_r^s-{\sum }_{i\in \mathcal{O}}{a}_{ir}{\pi }_i+\sum _{t\in \{0,...,⌊{\displaystyle \frac{T}{l_{\mathrm{m}\mathrm{i}\mathrm{n}}}}⌋\}}{b}_{tr}{\psi }_t,$$

(29)

This subproblem can be viewed as the Elementary Shortest Path Problem with Resource Constraints (ESPPRC) [25], typically solved by a dynamic programming labeling algorithm. However, for the PSVRSP-UD studied in this paper, the traditional labeling algorithm faces two challenges:

• the search space of the pricing subproblem corresponds to trips, and the time dual variable ${\psi }_t$ leads to different reduced costs for different departure times of the same order sequence;
• multi-scenario demand constraints require each label to track resource consumption and feasibility under all demand scenarios, significantly increasing label dimension and pruning difficulty.

To improve computational efficiency and avoid this combinatorial difficulty, this paper introduces the NG-Route relaxation strategy [26] into the labeling algorithm, combined with the group of labels and representative label mechanisms [27], reducing computational complexity without losing optimality.

4.2.1. Definition of NG-Labels

The ESPPRC is defined on a directed graph ${\mathcal{G}}^{‘}(\mathcal{O}\cup \{o,d\},\mathcal{A})$, where nodes $o$ and $d$ represent the onshore supply base as the path start and end, respectively. Under the NG-Route relaxation framework, a neighborhood set ${\mathcal{N}\mathcal{G}}_i\subseteq \mathcal{O}\cup \{d\}$ is predefined for each node $i\in \mathcal{O}$, consisting of the closest nodes to $i$. Path extension is only allowed within the neighborhood, thereby significantly compressing the state space. Meanwhile, to avoid generating numerous duplicate labels with the same path but different times due to different departure times, this paper further introduces the group of labels and representative label mechanism. The core logic includes:

• All labels in the same group correspond to an identical visiting sequence path, differing only in departure time;
• Only one representative label is retained per group, and the postponable time rt(∙) describes the flexible range of all feasible departure times within the group, such that every feasible trip in the same group can be derived from the representative label.

On this basis, an NG-Label characterizes the partial path state of the vessel under multi-scenario demand, time windows, and capacity constraints. Compared with traditional labels, the NG-Label not only records path information but also contains a neighborhood set to control the extension direction, further improving efficiency.

Each label $L_i\left(r{(L}_i\right),\eta {(L}_i),f{(L}_i),l{(L}_i),st{(L}_i),et{(L}_i),rt{(L}_i),\overline{c}{(L}_i),{NG(L}_i))$ stores the following data:

• $r{(L}_i)$ is the path from the origin to $\eta {(L}_i)$.
• $\eta {(L}_i)$ is the current node of the path $r{(L}_i)$.
• $f{(L}_i)$ is the parent label from which $L_i$ is extended.
• $l{(L}_i)=\{l^s{(L}_i\left)\right|s\in S\}$ is the vessel load in each demand scenario, where$l^s{(L}_i)$ denotes the vessel load upon arrival at node$\eta {(L}_i)$ along path $r{(L}_i)$ under scenario s, representing the most critical state information in the label.
• $st{(L}_i)$ is the start time for the trip.
• $et{(L}_i)$ is the trip end time (the departure time from the previous node).
• $rt{(L}_i)$ is the maximum time that the label can be postponed while remaining feasible.
• $\overline{c}{(L}_i)$ is the reduced cost of the path, including travel cost, dual values of nodes, and dual value contribution of the time interval from $st{(L}_i)$ to $et{(L}_i)$.
• ${NG(L}_i)$ is the neighborhood set of the current node $\eta {(L}_i)$.

4.2.2. Label Extension Rules

In the extension phase, each label extension is checked for feasibility according to the following constraints.

Neighborhood Constraint：Only extend to nodes in the neighborhood ${NG}_i$ of the current node $\eta {(L}_i)$.

$$j\in {NG(L}_i)$$

(30)

Time window constraint: Ensure the vessel’s arrival time does not exceed the latest service time of the node.

$$et{(L}_i)+{\tau }_{ij}\le b_j\forall j\in {NG(L}_i)$$

(31)

Load constraint: Vessel load must not exceed the maximum capacity under all scenarios.

$$l^s{(L}_i)+q_j^s\le Q \forall s\in \mathcal{S},j\in {NG(L}_i)$$

(32)

where $l^s{(L}_i)$ is the current load in scenario s.

If all constraints are satisfied, a new label $L_j\left(r{(L}_j\right),\eta {(L}_j),f{(L}_j),l{(L}_j),st{(L}_j),et{(L}_j),rt{(L}_j),\overline{c}{(L}_j),{NG(L}_j))$is generated with updated information.

$$r\left(L_j\right)=r\left(L_i\right)\cup \left\{j\right\}$$

(33)

$$\eta \left(L_j\right)=j$$

(34)

$$f\left(L_j\right)=L_i$$

(35)

$$l^s\left(L_j\right)=l^s\left(L_i\right)+q_j^s\forall s\in \mathcal{S}$$

(36)

$$st{(L}_j)=st{(L}_i)$$

(37)

$$et{(L}_j)=max\{a_j+s_j,et{(L}_i)+{\tau }_{ij}+s_j\}$$

(38)

$$rt{(L}_j)=min\{rt{(L}_i),b_j-et{(L}_j)\}$$

(39)

$$\overline{c}\left(L_j\right)=\overline{c}\left(L_i\right)+\sum _{s\in S}{p}_s{c}_{ij}^s-{\pi }_j+\sum _{\tilde{t}=ind\left(et{(L}_i\right))}^{\tilde{t}=ind\left(et{(L}_j\right))}{\displaystyle \frac{min\left(et{(L}_j\right),sup\left({∆}_{\tilde{t}}\right))-max(et{(L}_i),inf({∆}_{\tilde{t}}\left)\right)}{max(ϵ,sup({∆}_{\tilde{t}})-inf({∆}_{\tilde{t}}\left)\right)}}{\psi }_{\tilde{t}}$$

(40)

$${NG(L}_j)={\mathcal{N}\mathcal{G}}_j$$

(41)

Equations (33)-(41) update the information of the newly generated label $L_j$, including the path extended to the new node $\eta \left(L_j\right)$, the current node position, the parent label, the vessel’s load upon across all scenarios, the departure time from$\eta \left(L_j\right)$, the trip start time, the updated arrival and departure time at the new node, the maximum feasible postponement time, and the updated neighborhood set. In Equation (40), $ind\left(et{(L}_i\right))$ denotes the index of the discretized time interval that contains the time $\mathrm{e}\mathrm{t}(\cdot )$.Let $inf\left({∆}_{\tilde{t}}\right)$ and $inf\left({∆}_{\tilde{t}}\right)$represent the lower and upper bounds of the time interval ${∆}_{\tilde{t}}$, respectively. A small positive constant $\mathsf{ϵ}$ is introduced in order to avoid a division by 0 when $sup\left({∆}_{\tilde{t}}\right)-inf\left({∆}_{\tilde{t}}\right)=0$. With this definition, the dual variables ${\psi }_{\tilde{t}}$ are accumulated proportionally according to the proportion of time interval ${∆}_{\tilde{t}}$ included in [$et{(L}_i)$, $et{(L}_i)$].

4.2.3. Dominance Rule

In the labeling algorithm, the design of the dominance rule determines computational efficiency. The core idea is that if one label at the same node is not better than another in terms of resource consumption and cost, it cannot lead to a better solution and can be pruned early. Traditional dominance strategies are usually based on the deterministic shortest path framework. However, given the multi-scenario demand and temporal flexibility of labels in the PSVRSP-UD, the following dominance criterion is designed. A label $L_i^1$ dominates another label $L_i^2$ if all following conditions are satisfied simultaneously.

$$r\left(L_i^1\right)\subseteq r\left(L_i^2\right)$$

(42)

$$f\left(L_i^1\right)=f\left(L_i^2\right)$$

(43)

$$l^s\left(L_i^1\right)\le l^s\left(L_i^2\right)\forall s\in S$$

(44)

$$et\left(L_i^1\right)\le et\left(L_i^2\right)$$

(45)

$$rt\left(L_i^1\right)\ge rt\left(L_i^2\right)$$

(46)

$$\overline{c}\left(L_i^1\right)\le \overline{c}\left(L_i^2\right)$$

(47)

This dominance rule reduces the number of labels and improves the efficiency of the pricing subproblem without losing optimality.

4.3. Branching Strategy and Integer Solution Repair

Solving the RLMP via CG only obtains a linear relaxation solution. To obtain an optimal schedule satisfying integer constraints, CG must be embedded within a B&B framework to form a complete B&P algorithm. However, in multi-trip scheduling problems such as the PSVRSP-UD, the traditional arc branching strategy has obvious limitations. Even if all arc flows are integer-valued, trip variables may still be fractional. This occurs because the same visiting sequence can form multiple feasible trips at different departure times, making arc flows unable to uniquely determine the time schedule of trips. To address this issue, this section designs a two-level branching strategy tailored to the problem’s characteristics, as shown in Figure 2.

If a fractional-flow arc exists, the conventional arc branching strategy is applied. If all arc flows are integer, a repair procedure is initiated to construct an integer trip schedule without time conflicts. All feasible trip structures in the current fractional solution are treated as virtual customers, constructing an auxiliary VRPTW subproblem. Each virtual customer corresponds to a trip node sequence ${cus}_r$, with time window [$t_r^{start}$,$t_r^{start}+t_r^{relax}$] representing the feasible start interval of the sequence; service time and demand are set to 0; the fleet size is consistent with the original problem. The goal is to assign a conflict-free service sequence to these virtual customers, i.e., to construct a set of integer trips without overlapping times. The auxiliary VRPTW problem is solved by CG. If a feasible integer solution is obtained, the global optimal upper bound is updated and the current node is pruned. If repair fails, no integer solution can be formed based on the structures in the current solution, implying at least one arc in these structures does not belong to the optimal solution. A standard branching rule is applied on unconstrained arcs, creating two branches that force the arc flow to 0 and 1, respectively. If no feasible solution exists after all arcs with flow 1 are forbidden, the current node is directly pruned.

This two-level branching strategy effectively resolves the special issue of integer arc flows but fractional trip variables in multi-trip scheduling, guaranteeing the algorithm converges to the global optimal solution.

5. Numerical Experiments

To verify the effectiveness of the proposed model and algorithm, various computational experiments are conducted. The flow-based MILP model was solved using ILOG CPLEX 12.9, implemented in C++, with a time limit of 3600 seconds. All experiments were conducted on a Windows 10 system with an Intel(R) Core(TM) i5-8500 CPU @ 3.0 GHz and 8 GB RAM.

5.1. Test Instances and Parameter Settings

This section specifies the test instances and parameter settings adopted for PSVRSPUD. The benchmark instances are constructed based on real-world operational data from an oil and gas company operating in Brazilian waters, as reported in Silva’s work [28]. A total of 6 offshore platforms are selected to form the test network. The sailing time from the onshore supply base to each platform ranges from 16 to 19 hours, and the average sailing time between any two platforms is 1.1 hours. During the planning horizon, each platform generates a random number of service orders. For each order, the expected demand is uniformly distributed in the interval [70, 150]. To capture demand uncertainty, a scenario generation and reduction method is employed to construct a set of representative demand scenarios denoted by $\mathcal{S}$. Each scenario $s\in \mathcal{S}$ is assigned with an occurrence probability ${\mathrm{p}}_{\mathrm{s}}$. The actual demand under each scenario follows a normal distribution with the expected demand as the mean and a standard deviation $\mathsf{\sigma }\in \{\mathrm{0.1,0.2,0.3,0.4},0.5\}$. Each standard deviation corresponds to a distinct demand variation level, which enables a comprehensive evaluation of the impact of different demand uncertainty degrees on scheduling decisions. For time-related parameters, the service time windows are generated in accordance with the classic Solomon scheme. The earliest arrival time at each node is determined by the sailing distance, and the width of each time window ranges from 30 to 70 hours, which does not exceed the total planning horizon of 168 hours.

Each instance is named according to a consistent rule that includes the number of platforms, orders, vessels, and the standard deviation of demand. For illustration, the instance name c6-o12-v3-d3-5 indicates an instance with 6 platforms (c), 12 orders (o), 3 vessels (v), a demand standard deviation of 0.3 (d), and the replicate index (n) 5.

5.2. Experimental Results and Analysis

5.2.1. Performance of the B&P Algorithm

The performance of the proposed B&P algorithm is compared with the flow-based model solved by CPLEX 12.9. Experiments in this section are conducted under a fixed demand standard deviation ($\sigma =0.3$) and a weight coefficient ($\alpha =0.5$), considering 10 demand scenarios. In each instance, statistical comparisons are performed on multiple metrics, including objective function value (Obj), computational time (T), upper bound (UB), lower bound (LB), gap ($GAP={\displaystyle \frac{UB-LB}{LB}}\times 100\mathrm{\%}$), and relative deviation ($RD={\displaystyle \frac{{UB}_{B\&P}-{UB}_{CPLEX}}{{UB}_{CPLEX}}}\times 100\mathrm{\%}$).

The comparative results of the algorithms on different test instances are summarized in Table 2 and Table 3. For small-scale instances with 4-8 orders, both CPLEX and B&P can obtain optimal solutions within the 3600-second time limit. In terms of computational efficiency, the B&P algorithm exhibited a significant advantage. The objective values are exactly the same, and the relative deviation (RD) is 0%, which verifies that the optimal solutions obtained by the B&P algorithm are consistent with those of CPLEX. In terms of computational efficiency, the B&P algorithm demonstrates a remarkable advantage. For instances with 4 orders, the average solution time of B&P is only 0.012 seconds, while CPLEX requires approximately 0.25 seconds, making B&P more than 95% faster on average. For instances with 8 orders, the computational time of CPLEX increases sharply to a range of 522.5 seconds to 1586.91 seconds, with an average of 461.97 seconds. In contrast, the B&P algorithm solves the same instances in only 0.02 seconds to 1.7 seconds, with an average of merely 0.31 seconds. These results also indicate that B&P’s advantage in computational efficiency becomes more pronounced with modest increases in problem size.

For medium- and large-scale instances (involving 12 to 24 orders), CPLEX 12.9 fails to obtain optimal solutions within the 3600-second time limit. The average GAP of CPLEX reaches 97.38%, and this gap tends to increase with the instance size, resulting in a rapid decline in solution quality for larger instances—even failing to output valid UB and LB values for some large-scale instances. This indicates the limited ability of CPLEX to solve large-scale PSVRSP-UD using flow-based models. In contrast, the proposed B&P algorithm consistently obtains global optimal solutions, and the computed optimal solutions fall within the UB range reported by CPLEX. In terms of computational efficiency, the average runtime of the B&P algorithm is only 461.95 seconds, while CPLEX is forced to terminate at the 3600-second time limit. The RD analysis shows that the B&P algorithm achieves an average RD of -4.18%, indicating that its optimal solutions are even better than the upper bounds of CPLEX in most cases. For example, the instance c6-o20-v3-d3-4 achieves an RD of -12.74%, showing the most significant cost advantage, while instances c6-o12-v3-d3-3 and c6-o16-v3-d3-4 achieve an RD of 0.00%, indicating their optimal solutions are consistent with CPLEX’s upper bounds. These results collectively verify the effectiveness and practical relevance of the proposed approach in complex and uncertain offshore logistics environments.

In summary, the B&P algorithm significantly outperforms CPLEX in both solution quality and computational efficiency. It can obtain exact solutions within a reasonable computational time, exhibits robust convergence behavior, making it highly suitable for complex real-world PSVRSP applications.

5.2.2. Sensitivity Analysis on Number of Scenario

To investigate the influence of the number of demand scenarios on solutions under demand uncertainty, this section conducts a sensitivity analysis on the 4-port instance c6-o4-v2-d3-n, varying the number of scenarios from 10 to 100. The total cost and CPU time are illustrated in Figure 3. It shows that the total cost exhibits a significant downward trend as the number of scenarios increases from 10 to 100, rather than remaining robust as previously described. Specifically, the average total cost of the five instances decreases from 25,456,139.56 (10 scenarios) to 14,426,007.87 (100 scenarios), with a total reduction of approximately 43.33%. This indicates that increasing the number of scenarios can better capture the characteristics of demand uncertainty, and the optimal scheduling decisions are gradually adjusted to be more cost-effective as the scenario number increases.

In terms of computational efficiency, the algorithm’s runtime shows a stable positive correlation with the number of scenarios. For example, the CPU runtime of c6-o4-v2-d3-1 increases from 0.24 seconds to 2.61 seconds; for c6-o4-v2-d3-5, the runtime rises from 0.14 seconds to 1.94 seconds. On average, the runtime of the five instances increases from 0.214 seconds to 2.218 seconds, which is a tenfold increase but still remains at a low level, ensuring high computational efficiency.

Overall, the results indicate that the scenario setting adopted in this study can balance the trade-off between solution quality and computational efficiency.

5.2.3. Sensitivity Analysis on Demand Fluctuation Levels

To investigate the impacts of demand uncertainty on scheduling results, this section conducts a sensitivity analysis under different levels of demand fluctuation. Tests are carried out on 8-order instance c6-o8-v2-n as an example; the standard deviation of demand is varied. Specifically, σ is set to 0.1, 0.2, 0.3, 0.4, and 0.5, representing increasing uncertainty levels. All other parameters remain unchanged. The detailed results are reported in Table 4.

In terms of solution quality, the total cost generally shows an upward trend as demand uncertainty increases. As σ rises from 0.1 to 0.5, the average total cost increases from 52,513,500.63 to 52,998,710.12, indicating that higher demand volatility requires more conservative scheduling and additional capacity reserves to ensure feasibility, thereby raising the overall operational cost. At the instance level, most test problems show a gradual rise in objective value with increasing σ, reflecting the additional economic burden induced by demand uncertainty. Regarding computational efficiency, an overall decline in solution time can be observed as σ increases. The average computing time decreases monotonically from 1256.17 seconds at σ=0.1 to 800.82 seconds at σ=0.5. This phenomenon suggests that stronger demand uncertainty leads to a more constrained solution space and tighter feasibility conditions, which in turn accelerates the convergence of the algorithm and reduces the overall solution time without significantly compromising solution quality.

In summary, increasing demand fluctuation results in higher total expected cost due to the need for more robust scheduling arrangements, while simultaneously improving the computational efficiency of the proposed B&P algorithm.

5.2.4. Sensitivity Analysis on Weight Coefficient $\alpha$

To explore the impact of economic objectives and low-carbon objectives on the optimal scheduling scheme, a sensitivity analysis is conducted with respect to the weight coefficient$\alpha$. Tests are carried out on 8-order instance c6-o8-v2-d3-n as an example; all other parameters are held constant. Specifically,$\alpha$is set to five typical levels: 0,0.25,0.5,0.75,1, where a larger$\alpha$indicates greater emphasis on minimizing operational time cost, and a smaller $\alpha$ highlights the control of carbon emissions. The corresponding results of sailing time and sailing distance of the optimal scheduling scheme under different weight configurations are presented in Figure 4.

Experimental results show that the weight coefficient$\alpha$has a mild but stable impact on both sailing time and distance, without causing significant fluctuations or monotonic trends. When $\alpha$ varies from 0 to 1, the sailing time (time in best cost) maintains a relatively steady level across all 8-order instances. Specifically, the sailing time ranges approximately between 127.59 and 138.40, with no abrupt increases or decreases. For example, the time for instance c6-o8-v2-d3-2 remains 127.59 when$\alpha$=0.25 and$\alpha$=1, while the time for c6-o8-v2-d3-3 reaches the maximum of 138.80 when $\alpha =$0.25, showing only slight variations across different $\alpha$levels.

In terms of sailing distance (distance in best cost), it also remains stable with the change of $\alpha$. The distance ranges between 741.38 and 829.96, with no obvious monotonic trend as$\alpha$increases or decreases. For most instances, the distance fluctuates slightly within a narrow range across different$\alpha$levels; for example, the distance of c6-o8-v2-d3-4 varies only between 789.71 and 795.14, with a maximum variation of less than 0.7%.

Although $\alpha$adjusts the trade-off priority between operational time cost and carbon emission cost, the optimal scheduling scheme remains robust and efficient. The stable sailing time and distance across different$\alpha$levels indicate that adjusting the weight coefficient does not alter the core structure of the optimal route. In practical applications, decision-makers can flexibly select an appropriate$\alpha$according to corporate emission-reduction policies, operational budgets, or energy-saving requirements, without causing obvious deterioration in sailing efficiency or route rationality.

6. Conclusions

This study addresses the platform supply vessel routing and scheduling problem under demand uncertainty, aiming to minimize the expected total cost (including vessel operating costs and carbon emission costs) while satisfying multi-dimensional operational constraints such as vessel capacity, time windows, and multi-trip limits. A scenario-based MILP model is developed, which integrates stochastic offshore demand through multi-scenario methods. Specifically, demand uncertainty is represented by a finite set of representative scenarios generated using LHS combined with Cholesky decomposition, and the uncertainty is incorporated into the MILP by enforcing feasibility under all scenarios while minimizing the expected objective value weighted by scenario probabilities. The SBR algorithm is further adopted to prune redundant scenarios, ensuring computational tractability without compromising solution fidelity. To efficiently solve this NP-hard problem, a B&P algorithm is proposed that incorporates Dantzig–Wolfe decomposition, NG-Route labeling with group labels and representative label mechanisms, and a two-level branching strategy to accelerate convergence and ensure global optimality.

Systematic numerical experiments on instances generated from real-world offshore logistics data demonstrate that the proposed B&P method achieves optimal solutions and outperforms CPLEX in all test instances. For small-scale instances, the B&P algorithm obtains the same optimal solutions as CPLEX but with significantly higher computational efficiency. For medium- and large-scale instances, CPLEX fails to converge to optimal solutions within the time limit, whereas the B&P algorithm consistently achieves global optimality with a zero optimality gap. Sensitivity analyses reveal that a reasonable number of representative demand scenarios can effectively characterize demand uncertainty, and that increasing the number of scenarios improves solution quality without significant loss of computational efficiency. Higher demand fluctuation increases the expected total cost due to the need for more conservative scheduling but accelerates algorithm convergence by tightening the solution space. In addition, the weight coefficient balancing operational and low-carbon objectives has a mild but stable impact on scheduling results, allowing flexible adjustment according to practical decision needs. These results collectively verify the effectiveness and practical relevance of the proposed approach in complex and uncertain offshore logistics environments.

This work opens several directions for future research on the PSVRSP. First, the current modeling framework assumes that demands at different offshore platforms are mutually independent, whereas spatial or temporal correlations may exist in practical offshore operations. Incorporating demand correlations into scenario generation may better capture real-world demand dynamics. In addition, the current study focuses only on demand uncertainty and does not consider other common offshore uncertainties, such as adverse weather conditions. Future research may integrate multi-source uncertainties to further enhance the robustness of scheduling solutions. Furthermore, the algorithm’s performance may be improved by incorporating more advanced pruning strategies or parallel computing techniques to handle larger-scale instances with more complex operational constraints.