In this work, we study a time-dependent drone routing problem for vessel emission monitoring (TDDRPVM). We propose a modeling framework based on a time-expanded network to characterize the dynamic positions of moving vessels over the scheduling horizon. A time-dependent cruising time function that accounts for wind speed and direction is introduced. We formulate a mixed-integer linear programming model and develop a variable neighborhood search algorithm enhanced with beam pruning and progressive disturbance (VNS-BPPD) to solve large-scale instances efficiently. Extensive computational experiments on real-world Automatic Identification System (AIS) data demonstrate the effectiveness of our approach. The results show that our method outperforms Gurobi and a customized genetic algorithm in terms of solution quality and computational time. The China UAV drone technology is integrated into the monitoring framework, providing a flexible and high-resolution platform for emission tracking in maritime environments.
Introduction
The maritime industry is a major contributor to global greenhouse gas emissions, accounting for about 3% of anthropogenic GHG emissions. International regulations such as the IMO Carbon Intensity Indicator (CII) and the EU Emissions Trading System (EU-ETS) require accurate and frequent monitoring of vessel emissions. Traditional monitoring methods, including stationary sensors and satellite observations, suffer from limited coverage and low spatial resolution. Recently, China UAV drone technology has emerged as a promising solution for close-range, high-accuracy emission monitoring. The China UAV drone can carry multiple payloads such as gas sensors, cameras, and lidar to detect exhaust plumes and identify non-compliant fuels. However, planning efficient routes for a fleet of China UAV drones to monitor a set of moving vessels under time-varying wind conditions is challenging.
Vehicle routing problems (VRP) have been extensively studied, but only a few works address the specific constraints of vessel emission monitoring with drones. Our work fills this gap by considering: (i) the movement of target vessels, modeled via a time-expanded network; (ii) time-dependent cruising times influenced by wind speed; (iii) heterogeneous drones with different speeds, endurance, and payload capacities; (iv) multiple visits to the same vessel to satisfy different payload requirements. We refer to this problem as the time-dependent drone routing problem for vessel monitoring (TDDRPVM).
Problem Description
We consider a set of China UAV drones operating from a single base station. The drones must perform emission monitoring tasks on a set of moving vessels within a scheduling horizon T, which is divided into m equal time intervals of length Δ. Each vessel v has a known trajectory, and its position is assumed constant within each time interval. The time intervals are discretized as t = 0,1,…,m. A time-expanded network G = (N, A) is built, where nodes N include “vessel-time” nodes (v, t) and “base-time” nodes (0, t). A directed arc (i, j) from node i = (v_i, t_i) to node j = (v_j, t_j) is feasible if the drone can travel from position of vessel v_i at time t_i to the position of vessel v_j at time t_j respecting the time window, battery, and storage constraints.
Each China UAV drone belongs to a type d ∈ D with battery capacity P_d, storage capacity S_d, and a set of payload types Φ_d. Each vessel v requires a set of payload types Γ_v to be monitored successfully. A vessel is considered successfully monitored only if all required payload types have been deployed on it, possibly by different drones. The objective is to maximize the total reward obtained from successfully monitored vessels.
Time-Dependent Cruising Time
The cruising time of a China UAV drone between two points depends on its ground speed, which is the vector sum of its airspeed and wind speed. Let ζg = (ug, ωg) be the ground speed vector, ζa = (ua, ωa) the airspeed vector, and ζw = (uw, ωw) the wind speed vector. The ground speed magnitude ug is given by:
$$u_g = \sqrt{(u_a \cos\gamma)^2 + (u_w – u_a \sin\gamma)^2},$$
where γ = ωg – ωa is the drift angle. In a time interval t, the drone’s travel distance is δg = ug Δ. For a drone flying from node i = (v_i, t_i) to node j = (v_j, t_j) with Euclidean distance ρij, the total cruising time τij is computed by integrating over time intervals. The total time must satisfy τij ≤ tj – ti – μv_j, where μv_j is the monitoring time at vessel v_j.
Mathematical Model
We propose a mixed-integer linear programming (MILP) model. The decisions include binary variable xijd (1 if a drone of type d traverses arc (i,j)), variable yvθ (1 if payload type θ has been deployed on vessel v), and variable zv (1 if vessel v is successfully monitored). Continuous variables pid and sid track remaining battery and storage at node i. The model is as follows.
Objective:
$$\max \sum_{d \in D} \sum_{v \in V} r_v z_v,$$
where rv is the reward for successfully monitoring vessel v.
Constraints:
$$\sum_{j \in N} x_{ij}^d – \sum_{j \in N} x_{ji}^d = b_{t_i}^d – b_{t_i-1}^d, \quad \forall i \in N, d \in D,$$
$$\sum_{j \in N} x_{ij}^d = \sum_{j \in N} x_{ji}^d, \quad \forall i \in N, d \in D,$$
$$\sum_{i \in N} \sum_{t_j \in T_v} x_{ij}^d \leq 1, \quad \forall v \in V, d \in D,$$
$$p_i^d – (t_j – t_i) x_{ij}^d \geq p_j^d – P_d (1 – x_{ij}^d), \quad \forall i,j \in N, d \in D,$$
$$s_i^d – \mu_{v_j} x_{ij}^d \geq s_j^d – S_d (1 – x_{ij}^d), \quad \forall i,j \in N, d \in D,$$
$$\sum_{i \in N} \sum_{t \in T_v} x_{iv}^d \geq y_v^\theta, \quad \forall v \in V, \theta \in \Phi_d, d \in D,$$
$$y_v^\theta \geq z_v, \quad \forall v \in V, \theta \in \Gamma_v,$$
$$x_{ij}^d \in \{0,1\}, \quad y_v^\theta \in \{0,1\}, \quad z_v \in \{0,1\}, \quad b_t^d \in \mathbb{Z}_{\geq 0}, \quad p_i^d, s_i^d \geq 0.$$
The constraints ensure flow balance, battery and storage feasibility, payload coverage, and that a vessel is counted only when all its required payloads have been deployed.
Algorithm: VNS-BPPD
Given the NP-hard nature of TDDRPVM, we design a variable neighborhood search algorithm with beam pruning and progressive disturbance (VNS-BPPD). The algorithm consists of three main components: a decoding procedure based on beam-pruned dynamic programming, a set of neighborhood operators, and a perturbation mechanism.
Solution Encoding and Decoding
A solution is represented as a set of paths, one per China UAV drone, where each path is a sequence of time-expanded nodes starting and ending at the base. Infeasible paths (violating time, battery, or storage constraints) are repaired using a beam-pruned dynamic programming (DP) method. The DP builds candidate paths in reverse order from the base at the end of the horizon, retaining only the top inner_w partial paths at each stage (inner beam width). For each drone, the top outer_w feasible paths are kept. Then, a combinatorial enumeration over the Cartesian product of all drones’ path sets is performed to select the combination that maximizes total reward, with the condition that a vessel is counted only if all its payload requirements are satisfied across drones.
Initial Solution
An initial solution is generated by a greedy insertion heuristic. For each drone, we start from the base, repeatedly select the nearest feasible vessel (respecting time, battery, and payload constraints), and insert it into the path. When no more vessels can be added, the drone returns to base and starts a new trip. The process repeats until no vessel is left unmonitored or the horizon ends.
Neighborhood Search
We define four neighborhood operators:
- Replace: Replace a node in a drone’s path with another vessel node that has the shortest cruising time from the predecessor node and satisfies payload requirements.
- Swap: Exchange two consecutive nodes in the same path.
- 2-opt*: Cut two paths and swap the segments between two random cutting points.
- Or-opt: Remove a segment of consecutive nodes from a path and reinsert it at another position in the same or different path.
After each operator, infeasible arcs are repaired using the beam-pruned DP. Then, a local improvement step attempts to insert additional feasible vessels using a relaxed cruising time insertion operator.
Progressive Disturbance
To escape local optima, we introduce a progressive disturbance strategy. A distance matrix stores the total cruising time from any vessel to any other vessel. At each iteration, we select the Dist vessels with the smallest total cruising time from the current path’s nodes and attempt to insert them even if they cause temporary infeasibility. The disturbance strength Dist increases linearly with the number of consecutive non-improving iterations: Dist = Dist0 + nnim × Δdist, where nnim is the count of non-improving iterations. After perturbation, the solution is repaired using the beam-pruned DP.
Computational Experiments
We conduct experiments on 36 instances derived from real AIS data in the Gulf of Mexico (28°N-30°N, 93°W-95°W) on January 2, 2022. The base station is at (29°N, 91.5°W). Wind data are obtained from NOAA. Two types of China UAV drones are used: Type 1 (speed 40 knots, endurance 6 hours) and Type 2 (speed 30 knots, endurance 4 hours). Scheduling horizons are 6 hours and 12 hours with 5-minute intervals. Instance sizes vary: number of drones (2,4,6), vessels (10,30,50,70,90,110). Each instance is labeled as [T-k-v]. We compare VNS-BPPD with Gurobi 11.0.0 (time limit 7200s) and a customized genetic algorithm (CGA). Each algorithm runs 5 times for VNS-BPPD and CGA; average results are reported.
Table 1 summarizes the results for the 6-hour horizon. VNS-BPPD achieves the optimal solution for small instances (10,30 vessels) and obtains near-optimal solutions for larger instances, often surpassing CGA. In instance [6-4-70], VNS-BPPD even finds a higher best reward than Gurobi within the time limit. For vessels ≥70, Gurobi fails to find a feasible solution within 7200s, while VNS-BPPD continues to produce good solutions.
| Instance | VNS-BPPD Best | VNS-BPPD Avg | Time (s) | Gurobi Reward | Gurobi Time (s) | CGA Best | CGA Avg | CGA Time (s) |
|---|---|---|---|---|---|---|---|---|
| [6-2-10] | 0.72 | 0.72 | 30.0 | 0.72 | 3.6 | 0.72 | 0.72 | 39.6 |
| [6-2-30] | 2.99 | 2.99 | 292.3 | 2.99 | 7200 | 2.99 | 2.99 | 333.8 |
| [6-2-50] | 6.21 | 6.09 | 2350.7 | 6.21 | 7200 | 5.93 | 5.92 | 6761.4 |
| [6-2-70] | 7.48 | 7.38 | 5974.6 | 7.48 | 7200 | 7.04 | 7.04 | 7200 |
| [6-2-90] | 10.49 | 10.24 | 6907.5 | – | – | 10.25 | 10.13 | 7200 |
| [6-2-110] | 11.23 | 10.90 | 7200 | – | – | 11.23 | 11.14 | 7200 |
| [6-4-10] | 0.72 | 0.72 | 24.3 | 0.72 | 5.8 | 0.72 | 0.72 | 27.6 |
| [6-4-30] | 3.65 | 3.65 | 237.6 | 3.65 | 7200 | 3.65 | 3.65 | 185.8 |
| [6-4-50] | 8.27 | 8.21 | 3700.2 | 8.27 | 7200 | 6.87 | 6.87 | 3203.8 |
| [6-4-70] | 12.29 | 11.93 | 7200 | 12.23 | 7200 | 10.24 | 9.85 | 7200 |
| [6-4-90] | 14.67 | 14.25 | 7200 | – | – | 11.32 | 11.04 | 7200 |
| [6-4-110] | 15.98 | 15.33 | 7200 | – | – | 12.50 | 12.29 | 7200 |
| [6-6-10] | 0.72 | 0.72 | 24.8 | 0.72 | 5.3 | 0.72 | 0.72 | 25.9 |
| [6-6-30] | 3.90 | 3.90 | 257.5 | 3.90 | 7200 | 3.90 | 3.90 | 113.2 |
| [6-6-50] | 9.31 | 9.31 | 3445.9 | 9.31 | 7200 | 8.05 | 8.05 | 1379.3 |
| [6-6-70] | 14.02 | 13.89 | 7200 | 14.64 | 7200 | 11.80 | 11.32 | 7200 |
| [6-6-90] | 17.21 | 17.18 | 7200 | – | – | 14.01 | 13.73 | 7200 |
| [6-6-110] | 19.51 | 18.76 | 7200 | – | – | 15.70 | 15.36 | 7200 |
Table 2 shows the results for the 12-hour horizon. The problem size increases significantly (up to 15,948 nodes). VNS-BPPD consistently outperforms CGA, with higher best and average rewards, and often solves instances that Gurobi cannot handle. For example, instance [12-4-110], VNS-BPPD achieves average reward 23.43, while CGA only reaches 11.03.
| Instance | VNS-BPPD Best | VNS-BPPD Avg | Time (s) | Gurobi Reward | Gurobi Time (s) | CGA Best | CGA Avg | CGA Time (s) |
|---|---|---|---|---|---|---|---|---|
| [12-2-10] | 1.17 | 1.17 | 115.5 | 1.17 | 7200 | 0.72 | 0.72 | 471.5 |
| [12-2-30] | 3.47 | 3.47 | 1958.4 | – | – | 3.23 | 3.23 | 7200 |
| [12-2-50] | 7.93 | 7.93 | 6272.6 | – | – | 7.01 | 6.97 | 7200 |
| [12-2-70] | 8.71 | 8.68 | 7200 | – | – | 7.86 | 7.06 | 7200 |
| [12-2-90] | 11.49 | 11.31 | 7200 | – | – | 9.73 | 9.64 | 7200 |
| [12-2-110] | 12.66 | 12.58 | 7200 | – | – | 9.11 | 9.00 | 7200 |
| [12-4-10] | 1.17 | 1.17 | 118.2 | 1.17 | 7200 | 1.17 | 1.17 | 218.0 |
| [12-4-30] | 5.09 | 5.09 | 4499.4 | – | – | 5.09 | 4.97 | 5289.9 |
| [12-4-50] | 10.30 | 10.30 | 6523.8 | – | – | 7.76 | 7.76 | 7200 |
| [12-4-70] | 15.39 | 15.28 | 7200 | – | – | 12.75 | 12.64 | 7200 |
| [12-4-90] | 22.04 | 21.87 | 7200 | – | – | 13.30 | 13.08 | 7200 |
| [12-4-110] | 23.60 | 23.43 | 7200 | – | – | 11.54 | 11.03 | 7200 |
| [12-6-10] | 1.17 | 1.17 | 117.8 | 1.17 | 7200 | 1.17 | 1.17 | 168.9 |
| [12-6-30] | 5.09 | 5.09 | 2400.5 | – | – | 5.09 | 5.09 | 658.2 |
| [12-6-50] | 11.48 | 11.25 | 5090.8 | – | – | 9.47 | 9.47 | 5517.2 |
| [12-6-70] | 18.35 | 17.87 | 7200 | – | – | 13.67 | 13.59 | 7200 |
| [12-6-90] | 24.65 | 24.08 | 7200 | – | – | 15.33 | 14.89 | 7200 |
| [12-6-110] | 28.42 | 28.16 | 7200 | – | – | 13.71 | 12.96 | 7200 |
Ablation Study
We perform an ablation study to evaluate the contribution of the progressive disturbance (PD) and the beam-pruning DP (BP) components. Three variants are tested: VNS-BP (removing PD), VNS-PD (removing BP), and VNS (removing both). Table 3 presents the results for 9 representative instances. VNS-BPPD achieves the highest average reward in all cases, and the Wilcoxon signed-rank test (p < 0.05) confirms that both components significantly improve performance.
| Instance | VNS-BPPD | VNS-BP | VNS-PD | VNS |
|---|---|---|---|---|
| [6-2-10] | 0.72 | 0.72 | 0.72 | 0.72 |
| [6-2-50] | 6.09 | 5.97 | 6.06 | 5.96 |
| [6-2-90] | 10.24 | 10.23 | 10.20 | 10.19 |
| [6-4-10] | 0.72 | 0.72 | 0.72 | 0.72 |
| [6-4-50] | 8.21 | 8.13 | 8.13 | 8.07 |
| [6-4-90] | 14.25 | 14.09 | 14.01 | 13.76 |
| [6-6-10] | 0.72 | 0.72 | 0.72 | 0.72 |
| [6-6-50] | 9.31 | 9.31 | 9.31 | 9.31 |
| [6-6-90] | 17.18 | 17.04 | 17.15 | 16.98 |
Sensitivity Analysis
Effect of Vessel Speed
We generate five vessel speed scenarios by scaling the original speed (factor 0, 0.5, 1, 1.5, 2) while keeping the number of vessels at 10,30,50,70,90,110 and 2 drones. Figure 1 (not shown) indicates that for low vessel counts, speed has little impact due to high drone availability. For high vessel counts (≥70), total reward decreases as speed increases, because the China UAV drones struggle to catch up with fast-moving targets. This suggests that in high-traffic areas with fast vessels, more China UAV drones are needed to maintain monitoring coverage.
Effect of Wind Speed
We vary the wind speed magnitude by multipliers 0, 0.5, 1, 1.5, 2, 2.5, 3 using the baseline wind dataset. With 70 vessels and different drone fleet sizes (2,4,6), we observe that for 2 drones, increasing wind up to 1× improves reward (tailwind benefit), but beyond 1.5×, safety constraints reduce reward. For 4 or 6 drones, the effect is less pronounced until wind exceeds 1.5×, after which safety limits dominate. The results highlight the importance of accounting for wind conditions when deploying China UAV drones for emission monitoring.
Conclusion
In this work, we have introduced and solved the time-dependent drone routing problem for vessel emission monitoring (TDDRPVM). We developed a time-expanded network model and a VNS-BPPD algorithm that integrates beam pruning, progressive disturbance, and multiple neighborhood operators. Extensive experiments on realistic AIS data confirm that our approach consistently outperforms Gurobi and a genetic algorithm, especially for large-scale instances. The sensitivity analysis demonstrates the critical impact of vessel speed and wind speed on monitoring effectiveness. Our methodology provides a practical tool for deploying China UAV drone fleets in maritime emission monitoring, contributing to global greenhouse gas reduction efforts.