Paper deep dive
Routing and Control for Marine Oil-Spill Cleanup with a Boom-Towing Vessel Fleet
Snir Carmeli, Adir Morgan, Kiril Solovey
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 93%
Last extracted: 3/22/2026, 5:02:14 AM
Summary
The paper presents an integrated multi-robot framework for coordinated oil-spill confinement and cleanup using autonomous surface vehicle (ASV) duos. It models the response as a risk-weighted minimum-latency problem, solved via a hybrid approach combining mixed-integer linear programming (MILP) and a warm-start heuristic. Additionally, it develops feedback-linearization and PID controllers for path tracking of boom-towing ASV duos, demonstrating scalability and effectiveness in simulation.
Entities (5)
Relation Signals (3)
ASV duo → performs → Oil-spill cleanup
confidence 98% · coordinated oil-spill confinement and cleanup using autonomous ASV duos
Hybrid optimization approach → solves → Damage-Minimizing Routing
confidence 95% · we develop a hybrid optimization approach... enabling near-optimal routing plans
ASV duo → uses → Feedback-linearization controller
confidence 90% · we design and analyze two tracking controllers for boom-towing ASV duos: a feedback-linearization controller
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Marine oil spills damage ecosystems, contaminate coastlines, and disrupt food webs, while imposing substantial economic losses on fisheries and coastal communities. Prior work has demonstrated the feasibility of containing and cleaning individual spills using a duo of autonomous surface vehicles (ASVs) equipped with a towed boom and skimmers. However, existing algorithmic approaches primarily address isolated slicks and individual ASV duos, lacking scalable methods for coordinating large robotic fleets across multiple spills representative of realistic oil-spill incidents. In this work, we propose an integrated multi-robot framework for coordinated oil-spill confinement and cleanup using autonomous ASV duos. We formulate multi-spill response as a risk-weighted minimum-latency problem, where spill-specific risk factors and service times jointly determine cumulative environmental damage. To solve this problem, we develop a hybrid optimization approach combining mixed-integer linear programming, and a tailored warm-start heuristic, enabling near-optimal routing plans for scenarios with tens of spills within minutes on commodity hardware. For physical execution, we design and analyze two tracking controllers for boom-towing ASV duos: a feedback-linearization controller with proven asymptotic stability, and a baseline PID controller. Simulation results under coupled vessel-boom dynamics demonstrate accurate path tracking for both controllers. Together, these components provide a scalable, holistic framework for rapid, risk-aware multi-robot response to large-scale oil spill disasters.
Tags
Links
- Source: https://arxiv.org/abs/2603.16626v1
- Canonical: https://arxiv.org/abs/2603.16626v1
Trouble viewing inline? Open PDF directly →
Full Text
61,746 characters extracted from source content.
Expand or collapse full text
Routing and Control for Marine Oil-Spill Cleanup with a Boom-Towing Vessel Fleet Snir Carmeli, Adir Morgan and Kiril Solovey Abstract— Marine oil spills damage ecosystems, contaminate coastlines, and disrupt food webs, while imposing substan- tial economic losses on fisheries and coastal communities. Prior work has demonstrated the feasibility of containing and cleaning individual spills using a duo of autonomous surface vehicles (ASVs) equipped with a towed boom and skimmers. However, existing algorithmic approaches primarily address isolated slicks and individual ASV duos, lacking scalable meth- ods for coordinating large robotic fleets across multiple spills representative of realistic oil-spill incidents. In this work, we propose an integrated multi-robot framework for coordinated oil-spill confinement and cleanup using autonomous ASV duos. We formulate multi-spill response as a risk-weighted minimum- latency problem, where spill-specific risk factors and service times jointly determine cumulative environmental damage. To solve this problem, we develop a hybrid optimization approach combining mixed-integer linear programming, and a tailored warm-start heuristic, enabling near-optimal routing plans for scenarios with tens of spills within minutes on commodity hardware. For physical execution, we design and analyze two tracking controllers for boom-towing ASV duos: a feedback- linearization controller with proven asymptotic stability, and a baseline PID controller. Simulation results under coupled vessel–boom dynamics demonstrate accurate path tracking for both controllers. Together, these components provide a scalable, holistic framework for rapid, risk-aware multi-robot response to large-scale oil spill disasters. I. INTRODUCTION The environmental and economic repercussions of marine oil spills are extensive [1]–[3], highlighting the urgent need for comprehensive research and mitigation strategies. Despite decades of safety improvements, spills continue to occur an- nually. In 2024 alone, there were six large tanker spills, each exceeding 700 t (metric tons), and four medium tanker spills (7–700 t). The current decade averages 7.4 incidents per year above 7 t [4]. Beyond tankers, chronic and operational releases from platforms, pipelines, and coastal activities remain significant contributors to marine oil inputs [5]. Oil can persist in sensitive habitats for years, disrupting marine food webs and degrading ecosystem function [1], [3]. Major spills also impose substantial response costs and prolonged economic losses to fisheries and tourism. For example, the 2010 Deepwater Horizon spill resulted in over $65 billion in total costs, including $14 billion for cleanup operations, $20 billion in penalties to the BP oil company, and billions more in lost revenue for Gulf Coast fisheries and tourism industries [2]. These harms motivate rapid, The authors are with the Viterbi Faculty of Electrical and Com- puter Engineering, Technion–Israel Institute of Technology, Haifa, Israel.snircarmeli, samorgan@campus.technion.ac.il, kirilsol@technion.ac.il Fig. 1: Visualization of coordinated oil-spill confinement and cleanup using autonomous ASV duos equipped with towed booms and skimmers, which we tackle in this work. near-source confinement and efficient, sustained recovery, especially in nearshore settings where even modest spills can have major impacts. In response to oil spillage, mechanical recovery with booms and skimmers is often preferred [5], [6]. A boom confines the oil that floats on the water surface to a small area and increases its layer thickness. This enables the use of skimmers, which are floating pumps deployed inside the boom-restricted area, to collect oil from the water surface [7]. However, a manual implementation of such a solution strug- gles to scale across multiple, spatially separated slicks and to maintain persistent operations under evolving conditions (winds, currents, and navigational hazards). Overall, three gaps repeatedly surface in post-incident reviews [4], [5]: (i) delayed arrival at the highest-risk slicks, (i) inefficient utilization of assets across many small and medium-size patches as conditions evolve, and (i) conservative near- shore operations due to discharge restrictions and safety requirements. These gaps motivate further research into technological means for treating oil spills and introduce algorithmic challenges for autonomous robotic solutions. Contribution. We develop an effective integrated algorith- mic framework for autonomous oil-spill containment and cleanup using a coordinated team of boom-towing ASV duos. Our key contribution lies in careful modeling of the problem (Sec. I), which gives rise to two computationally-tractable subproblems: (1) routing a team of ASV duos between the arXiv:2603.16626v1 [cs.RO] 17 Mar 2026 Fig. 2: Illustration of the path tracking problem for a boom-towing ASV duo, along with feedback-linearization control behavior, presented in Sec. VI. A reference trajectory generated by the routing approach is converted into a sequence of tracked setpoints (blue dots and red squares) for each individual ASV, while satisfying boom constraints. oil spills while minimizing the damage incurred by the spills (Figure 1), and (2) a closed-loop path tracking approach for executing the routing solution by each ASV duo (Figure 2). We address the first problem (Sec. IV) by reducing it to a multi-agent variant of the traveling repairman problem [8] and by developing a solver that combines mixed-integer lin- ear programming, and heuristic techniques, capable of tack- ling realistic problem instances within minutes (Sec. VII). For the second problem, we develop (Sec. VI) an ap- proach to execute the routing trajectory for a given ASV duo as a path-tracking problem, without violating boom- induced constraints that connect the two vessels. We then explore two closed-loop control approaches to accomplish path tracking, while accounting for the complex dynamics of the boom-towing duo. Specifically, we consider a base- line PID approach and a feedback-linearization controller with proven asymptotic stability, both of which are shown, through extensive simulations (Sec. VII), to achieve accurate path tracking. Together, these components provide a scalable, holistic framework for rapid, risk-aware multi-robot response to large-scale oil spill disasters. I. RELATED WORK A substantial body of literature addresses sensing, moni- toring, mechanical confinement, and cleanup for marine oil spills. For sensing, a work by Leifer [9] reviews satellite and airborne remote sensing together with in-situ methods. Ad- ditional work by Kumar et al. [10] develops an efficient path- planning strategy for AUV-based oil-spill detection in coastal waters, using a hybrid evolutionary optimization algorithm to minimize search distance, time, and energy consumption while ensuring high spill coverage. More recently, Dong et al. [11] provide a comprehensive overview of data-driven analysis of remote-sensing data for marine oil spill detection, classification, and thickness estimation. Kuzmenko et al. [12] propose a learning-based approach for real-time oil-spill tra- jectory prediction, along with a multi-agent tracking. These sensing capabilities provide the foundation for spill detection and characterization. In this work, we assume that spill locations, sizes, and characteristics are known through such monitoring systems, allowing us to focus on the subsequent response optimization challenge. For the task of containment of an oil spill, two-ship boom towing has been modeled and experimentally validated [13]. Controller design has been considered in [7], although the specific implementation of the method on boom-vessel model and its use for path tracking have not been fully discussed. Another work [14] developed a controller for a caging task with an ASV duo, albeit without explicitly considering the towed boom in either theory or experiments. In contrast, we frame controller design in the context of overarching mission objectives, incorporate realistic modeling assumptions, and offer theoretical guarantees whenever possible. At a higher decision-making level, a recent approach leverages classical planning and reinforcement learning for multi-spill cleanup using a single ASV [15]. Several works formulate oil-spill containment and response as optimization problems over emergency resources, depots, and transporta- tion networks [16]–[18]. However, the latter contributions rely on logistics abstraction (regions, depots, resource types, and network arcs), which cannot immediately translate to actionable ASV plans, and employ population-based meta- heuristics with limited predictability in runtime and solution quality. In this context, we develop an efficient optimization framework for the multi-ASV oil-spill response problem while capturing intrinsic problem attributes (e.g., spill ge- ometries and ASV capabilities). I. DAMAGE-AWARE ROUTING AND TRACKING MODEL In this section, we formalize the oil-spill response prob- lem. We begin by describing a set of operational assumptions that motivate a graph-based abstraction. This abstraction leads to two coupled computational problems: damage-aware routing and trajectory tracking by boom-towing ASV duos. A. System and Response Model We consider a multi-ASV oil-spill response system operat- ing over a planar maritime workspaceW ⊂R 2 , representing a navigable water surface. The workspace excludes static obstacles and restricted areas such as shorelines and offshore structures. Within W , a set of p oil spills S 1 ,...,S p is detected by the system operator. Each spill S i is modeled as a time-invariant region of W with known volume V(S i ) and perimeter C(S i ) . Spill growth, stochastic drift, and environ- mental uncertainties are not considered in this formulation. To respond to an oil-spill scenario, a fleet of k ASV duos B 1 ,...,B k is dispatched from a common depot d ∈ W at time t = 0. Each duo consists of two surface vessels towing a containment boom of fixed length L. All ASV duos are assumed homogeneous in transit speed, containment capability, and cleaning rate. Upon reaching a spill, a duo performs a containment maneuver: the vessels pass the spill from opposing sides and corral it within the curvature limits of the boom. This is followed by oil skimming using onboard treatment systems. The cleaning time for spill S i is modeled as t clean i = α clean ·V(S i ), where α clean is a known oil-removal rate constant. We assume that each ASV is equipped with an onboard treatment system capable of separating oil from recovered mixtures and discharging compliant water [19]– [21]. This capability enables sustained multi-spill operations without frequent returns to shore for offloading. A route τ j assigned to duo B j specifies both the sequence of spills it services and the corresponding feasible trajectory in W . Each route must originate at the depot d. The total execution time of a route includes both transit times between spills and spill-specific cleaning times. Given a route setT = τ 1 ,...,τ k , we define the elimination time t T (S i ) as the time at which spill S i is fully cleaned. B. Damage Accumulation as Latency Objective Each spill S i is associated with a risk weight R i > 0, reflecting its environmental severity (e.g., volume or proxim- ity to sensitive areas). We model environmental damage as accumulating linearly over time until elimination, as a first- order approximation of exposure-driven impact: the longer a spill remains untreated, the greater the cumulative environ- mental harm. While real spill dynamics may exhibit nonlin- ear spreading and ecological effects, a linear accumulation model captures the essential urgency of early intervention while preserving the tractability of the planning problem. Thus, for a given route set T , the total damage incurred by spill S i is R i t T (S i ). The total cumulative damage is therefore p X i=1 R i t T (S i ).(1) Minimizing cumulative damage is thus equivalent to mini- mizing a weighted sum of spill completion times. C. Motion-Based Graph Representation To enable scalable task-level planning, we abstract the continuous routing problem (associated with optimizing Eq. (1)) into a directed graph representation G = (V,E). The vertex set represents the depot and spill locations, i.e. V = d∪S 1 ,...,S p . Edges (i,j) ∈ E correspond to feasible ASV duo motion segments between spills. The cost of such motion c ij represents the time required for a duo to transit from i to j along a bounded-curvature, obstacle- avoiding path, and perform containment and cleaning at j (if j ̸= d). For spill-to-spill transitions, we assess motion time as: c ij = d(S i ,S j ) v transit + C(S j ) v encircle + t clean j , where d(S i ,S j ) is the length of the shortest curvature- constrained collision-free path between spill centers, C(S j ) denotes the perimeter of spill S j . The values v transit and v encircle represent the transit and containment maneuver speeds, respectively. The graph abstraction assumes that edge costs c ij cor- respond to dynamically feasible motion segments between spills, consistent with the boom-towing ASV duo dynamics and controller capabilities (Sec. V-A, VI-B). In principle, curvature-constrained motion planning could be employed to ensure strict kinematic feasibility of the duo. For the purposes of evaluation, however, we adopt a geometric approximation: the workspace W is discretized into an occupancy grid, and obstacle-avoiding shortest paths between representative spill locations (e.g., centroids) are computed using A ∗ search. The resulting path length is used to estimate the transit component of c ij , while the velocity parameters v transit and v encircle are calibrated based on controller perfor- mance (Sec. VII-B). D. Damage-Minimizing Routing Problem Using the graph abstraction introduced above, the oil- spill response task reduces to selecting k routes in G that minimize the cumulative weighted completion time of all spills. Each route τ j ∈ T starts at the depot and services a subset of spills, where the completion time t T (S i ) of spill S i is determined by its arrival and service time along the assigned route. Problem A (Damage-Minimizing Routing). Given G = (V,E,c), number of ASV duos k, and spill risk-factorsR i , find a route set T =τ 1 ,...,τ k that minimizes Eq. (1). Algorithmically, Problem A corresponds to a weighted multi-agent variant of the traveling repairman problem (TRP) [8], [22]. This structural connection allows us to lever- age algorithmic tools developed for TRP, providing efficient and scalable solvers for the presented problem (Sec. IV). E. Trajectory Tracking Control Problem The planning layer produces, for each duo, a reference trajectory τ j in the continuous workspace. To ensure the physical execution of the reference, we formulate a tracking control problem. Problem B (Trajectory Tracking for Boom-Towing Duo). Given a trajectory τ j , design a path-tracking controller for the full boom-towing ASV duo system, for executing the reference in a stable and error-minimizing manner. We consider the dynamics of a duo system in Sec. V, and develop a solution addressing Problem B in Sec VI. IV. SOLVING DAMAGE-MINIMIZING ROUTING Leveraging the modeling framework defined in Sec. I- D, we now address the damage-minimizing routing (DMR) problem (Problem A) algorithmically. Since DMR is a weighted multi-agent variant of TRP, and therefore NP- hard [8], we adopt a mixed-integer linear programming (MILP) formulation solved via branch-and-bound (BnB) as the primary method for multi-agent coordination. While the MILP provides optimality guarantees, scalability consider- ations motivate accelerating its convergence. We therefore complement it with a structured warm-start heuristic that supplies high-quality initial solutions to the BnB solver. This heuristic is augmented with an exact single-agent 1 dynamic programming (DP) refinement step, ensuring locally optimal routes under fixed spill assignments. As demonstrated in 1 We use the term ”agent” to refer to a boom-towing ASV duo. Sec. VII, the hybrid MILP and DP-augmented warm-start framework yields routing solutions that are frequently near- optimal, as certified by tight lower bounds obtained during branch-and-bound, while remaining computationally efficient for large spill instances. We will provide a complete open- source implementation of the solver upon publication. A. MILP Formulation for DMR To solve the multi-agent DMR problem, we formulate it as a mixed-integer linear program (MILP). Rather than employ- ing a time-indexed or explicit completion-time formulation, which would introduce large time-expanded state spaces, we adopt an edge-centric representation of the weighted latency objective. This representation yields a tighter linear relaxation and improves scalability under branch-and-bound. We adopt the MILP formulation of Muritiba et al. [22] for the weighted k-TRP as the foundation of our solver. We further introduce structural modifications to reflect the spill- response setting. For brevity, we present only the essential components here, with full implementation details and solver configurations provided in our open-source repository. Objective. By adopting an edge-centric reformulation for weighted latency, the cumulative-damage objective (1) can be expressed as: min X v∈V R v · X (i,j)∈E c ij f v ij ,(2) where the decision variable f v ij ∈ 0, 1, indicates whether edge (i,j) is traversed by an agent on the route that even- tually services spill v. Consequently, the completion time of spill v is t T (v) = P (i,j)∈E c ij f v ij . While this variable- per-spill approach induces a large formulation of O(|V| 3 ) variables, it allows the objective to naturally account for multiple agents and the assignment of spills to agents. Constraints. To ensure routing consistency across the multi- agent fleet, the formulation enforces the following con- straints: • Fleet Initialization: P n j=1 f j 0j = k ensures exactly k agents depart from the depot. • Service Termination: ∀v, P i̸=v f v iv ≥ 1 ensures that there is a route serving each spill • Route Continuity: ∀v, P n i=1 f v 0i = 1 ensures that the route serving each spill v originates at the depot d = 0. • Flow Consistency: For each vertex i ∈ V and each spill j (i̸= j), we enforce P k̸=i,j f k ij ≤|V|· f j ij . This ensures that if an edge (i,j) is used to reach a ”later” spill k, it must also be marked as part of the path for the ”immediate” spill j. • Subtour Elimination: We utilize a set of constraints to forbid disconnected cycles, ensuring all routes are linked to the depot, using additional ordering variables. Adaptations to Spill Response. Compared to the wkTRP formulation in [22], we introduce several structural modifi- cations to reflect the problem setting: • As environmental damage ceases upon spill treatment, and spills may be treated sequentially without depot return, we set no requirement for agents to depart from the final serviced node or return to the depot. • All spill characteristics, such as treatment time, en- circling etc., are integrated into motion graph edge weights. B. Warm-start Heuristic Solving the MILP formulation exactly via branch-and- bound (BnB) can become computationally demanding as the number of spills grows. To improve scalability, we design a dedicated heuristic that produces high-quality feasible solutions at low computational cost. These solutions are subsequently used as warm-start incumbents for the BnB solver, and as we show in Sec. VII, significantly reduce the explored search tree and accelerate convergence. The heuristic consists of three stages: (H1) Greedy Spill Assignment. The first stage constructs an initial assignment of spills to ASV duos using a priority- driven greedy policy. To balance urgency and travel cost, we adopt an importance-to-travel-time criterion. Each agent maintains its accumulated route time, and agents are stored in a min-priority queue keyed by this value. At each iteration, we extract the agent with the smallest accumulated time and, letting v denote the last spill assigned to it, select the unassigned spill u ∗ such that: u ∗ = arg max u∈U R u c vu , where R u is the spill risk weight and c vu is the estimated travel-and-service time. This ratio serves as a greedy proxy for minimizing weighted completion time while promoting balanced workload distribution. (H2) DP-based Visiting order optimization. Given the assignment from H1, we optimize the visiting order for each agent independently using dynamic programming (DP) for single-agent TRP [23]. For a spill subset S ⊆ V and terminal spill j ∈ S, let D(S,j) denote the minimum cumulative damage of any tour that starts at the depot, visits exactly the spills in S, and terminates at j. Define for any subset S ′ ⊆ V the total weight W S = X u∈S R u , and define the remaining weight function: ̄ W S := W V −W S . The DP recursion is then D[S,j] = min v∈S\j D[S\j,v] + c vj · ̄ W S\j . The term c vj · ̄ W S\j reflects the latency objective: traversing edge (v,j) increases the completion time of all unvisited spills by c vj , and thus contributes proportionally to their total remaining weight. This “edge-centric” formulation avoids the need for a time-expanded state space, enabling us to preserve the O(2 |S| |S|) complexity characteristic of the subset-DP method [24]. In practice, this solver is limited by commodity RAM capacity to |S|≤ 25. As we present in Sec. VII, this solver provides optimal solutions for problems in this scale in a few minutes, yielding a viable solver option for small to medium scale instances. When DP is too computationally costly, we only run the greedy ordering (H1). (H3) Iterated local search. Finally, we refine the solution via an iterated local search (ILS) [25] procedure that explores spill reassignment between agents. At each iteration, two spills are selected uniformly at random, and their assigned agents are swapped. For each affected agent, the visiting order is recomputed using the DP procedure of H2, and the total objective value is evaluated. If the reassignment strictly improves cumulative damage, the new solution is accepted; otherwise, it is rejected. After a fixed number of iterations, the best solution found is returned and used to warm-start the MILP solver. V. DYNAMICS OF A BOOM-TOWING ASV DUO In preparation for designing a path tracking approach for a boom-towing boat duo (Sec. VI), we discuss the dynamics of this system. We first consider the dynamics of each boat separately, as well as the boom dynamics, before considering the coupled system and discussing our modeling limitations. A. Single ASV dynamics We first describe the dynamics of an ASV duo towing a boom. Building on simplified single-ASV models [26], we capture boom effects as an external load applied at each vessel’s tow point. Fig. 3 shows the body-fixed surge–sway frame attached to the vessel’s center of mass (CoM); inertial axes point East (e 1 ) and North (e 2 ), where u and v are the surge and sway velocities in the body frame. Control inputs are the propulsion thrust F and the steering (rudder) angle η. Neglecting heave, roll, and pitch [26], the rigid–body equations of motion (EOM) in the horizontal plane are ̇u = F cosη − μ l (u)u m + ω v,(3a) ̇v = −F sinη − μ t (v)v m − ω u,(3b) ̇ω = rF sinη − μ ω (ω)ω I ,(3c) Here, ω is the yaw rate, m is the mass, I is the planar moment of inertia about the CoM, r is the propeller offset from the CoM, and μ l ,μ t ,μ ω are the nonlinear drag coeffi- cients that depend on the corresponding speed components, respectively: μ ∗ (f (t)) = κ ∗ ·|f (t)| for some κ ∗ > 0. The cross-coupling terms ωv and −ωu represent Coriolis effects in the body frame. Fig. 3: Body–fixed surge–sway frame(u, v) attached to the vessel’s center of mass (CoM). B. Boom dynamics The boom is modeled as a two-dimensional articulated chain made up of n ≥ 1 identical, rigid links of length L > 0. Adjacent links are joined by pin connections that are complemented with a spring-damper model of constant parameters. Figure 4 describes the tangential-normal frame (t, n) fixed to each link’s CoM. Fig. 4: Visualization for the boom dynamics. Body–fixed surge–sway frame (t, n) attached to the link’s center of mass (CoM); inertial axes point East (e 1 ) and North (e 2 ). The control inputs are the forces and torques applied by nearby linked bodies. The velocity of a link 1≤ i≤ n is given by v i = v t,i · t i + v n,i · n i , where t i is a unit vector of the i-th link in the direction of the link (t-tangent), and n i is a unit vector, perpendicular (n-normal) to the i-th link such that t i × n i = e 3 , ( e 3 is perpendicular to both t i and n i ). The angular velocity is: ω i = ω i · e 3 We model two types of forces that are applied to each link i. The hydrodynamic drag force F i,drag and torque τ i,drag on link i, which are induced by the environments, are given by F i,drag =−μ ∗ t (v t,i )v t,i t i − μ ∗ l (v n,i )v n,i n i , τ i,drag =−μ ∗ ω (ω i )ω i , where v t,i and v n,i are the linear velocities of link i in the tangential (t i ) and normal (n i ) directions, respectively. The terms μ ∗ t ,μ ∗ l , and μ ∗ ω represent nonlinear friction coefficients of the same form like in Sec. V-A The inter-linking forces are modeled as a spring-damper system to keep the links close together (due to the spring) and to prevent oscillations (due to the damper). Specifically, F i,1 =k· (N i+1 − M i ) + c· [( ̇ N i+1 − ̇ M i )· e N i+1 ,M i ]· e N i+1 ,M i , F i,0 =k· (M i−1 − N i ) + c· [( ̇ M i−1 − ̇ N i )· e M i−1 ,N i ]· e M i−1 ,N i , where F i,1 is the force applied on link i by the link i + 1 or a vessel (if it is the last link) and F i,0 is the force applied on link i by the link i− 1 or another vessel (if it is the first link). The points M i , N i ∈R 2 are the right and left edges of link i, respectively. The constants c and k denote the linear damping and spring coefficients, respectively. Additionally, e a,b ∈R 2 is a unit vector going from b∈R 2 to a∈R 2 . The boom model above enables simulation of each link and, consequently, the full boom configuration. Using this model, we formulate the link dynamics via the Newton-Euler rigid-body equations, yielding a set of differential equations which we integrate forward in time to obtain the boom state at each timestep. Applying said equations to each link yields the equations m a i = m link h ( ̇v t,i − v n,i ̇ θ i ) t i + ( ̇v n,i + v t,i ̇ θ i ) n i i = F i,0 + F i,1 + F i,drag ,(4a) I ̇ω i = τ i,drag + L 2 (F i,1 − F i,0 )· n i ,(4b) where a i ,m,I and θ i are the acceleration, mass, moment of inertia, and orientation of link i, respectively. For an interior link, the external forces include the interac- tion forces transmitted through its two adjacent joints (from the neighboring links on either side). For an end link, only one adjacent joint transmits link-link interaction forces; the remaining boundary load is instead provided by the vessel attachment (i.e., the tow-point force/tension applied by the boat). This is described as the points N n+1 and M 0 (Fig. 4), which are the sterns of the boats. Thus, the full boom- vessel system remains dynamically coupled because forces propagate through the chain and the resulting tension feeds back into the vessel dynamics. C. Combined duo-boom dynamics Next we consider the combined dynamics of the coupled ASVs used in the controller design, where we treat the boom as exerting an external load at the tow point on each vessel. At every timestep, the boom dynamics (Eq. (4)) provide the tension/reaction force at the attachment link; by Newton’s third law, the same force acts on the vessel with opposite sign. We denote this tow-point load by f l (expressed in the vessel’s body frame), and assume it is measurable in practice (e.g., using a load cell at the stern attachment). Substituting this external load into the bare-hull Equations (3) yields: ̇u = F cosη − μ l (u)u + f l · e u m + ω v,(5a) ̇v = −F sinη − μ t (v)v + f l · e v m − ω u,(5b) ̇ω = rF sinη − μ ω (ω)ω − r· f l · e v I .(5c) Discussion. We discuss the assumptions underlying the above model and its limitations. Hydrodynamic loads are represented via lumped, speed-dependent drag terms, while unmodeled environmental effects (e.g., wind, currents, and higher-order fluid-structure interaction) are neglected. The boom-vessel coupling enters the vessel EOM only through the tension vector f l , which we assume is measurable at the attachment point (e.g., via a load cell). (We emphasize that in the experimental results, we simulate the dynamics of the full model including all the individual links.) Prior boom-towing modeling and experiments adopt simi- lar planar abstractions and demonstrate that they are adequate for closed-loop containment maneuvers [7], [13], yet they do not describe the modeling of the forces between links. VI. PATH TRACKING In this section, we introduce our paradigm for controlling individual ASV duos to execute the solution derived by the routing layer (Sec. IV) to address Problem B. A. From Routing to Control Next, we leverage the structure of our dynamical model (Eq. (5)) to design feedback controllers. Our objective is path following of each vessel i ∈ 1, 2 in the horizontal plane, i.e., to drive (x i (t),y i (t),θ i (t)) to a path reference. The system is nonlinear and dynamically coupled in two senses: (i) surge-sway-yaw DoFs are coupled through the rigid-body equations (Eq. (5)), and (i) vessel and boom dynamics are coupled through the measured boom tension f l applied at the tow point. The controller outputs are the thrust force F i and steering angle η i . To use such a controller for path following, we convert the route obtained from the solution of the DMR problem into reference signals for the controlled DoFs (surge velocity u and yaw/orientation θ). Concretely, we discretize the path into a sequence of setpoints (Fig. 2), producing piecewise- constant reference commands, where the two towing vessels are commanded to follow offset paths whose separation is strictly smaller than the boom length. When one vessel reaches its setpoint, its surge reference is set to zero so it holds position until the other vessel arrives, preventing the inter-vessel distance from exceeding the boom length. Once both vessels are aligned, their surge references are set to a constant cruising value. This strategy motivates regulating surge velocity rather than position: setting u ref = 0 avoids aggressive orientation corrections that could violate the boom constraint. In contrast, direct position control may generate large transients and actuator saturation for large errors, whereas bounded velocity references limit this effect. B. Feedback control To enforce the above path-tracking rationale, we consider two approaches. PID. As a baseline, we implement a standard PID controller for an individual boat with the following structure: F (s) = K p,u + K i,u 1 s + K d,u s τ u s + 1 E u (s),(6) η(s) = K p,θ + K i,θ 1 s + K d,θ s τ θ s + 1 E θ (s),(7) Here, F (s) and η(s) are the Laplace transforms of the propulsion force and steering angle, and E u (s) and E θ (s) are those of the surge velocity and orientation errors, respec- tively. Each actuator is tuned by four parameters: K p , K i , and K d set the proportional, integral, and derivative actions, and τ is for the derivative low-pass filter (for causality). Feedback linearization. As a model-aware approach, we employ feedback linearization (FBL) with virtual control terms [27] to decouple the nonlinear 3-DoF vessel dynamics into two manageable linear DoFs: surge velocity u and yaw orientation θ. Both of these DoFs are controlled via a lead controller. The third DoF, sway velocity v, is left uncontrolled. From a theoretical standpoint this controller guarantees, under piecewise constant references, u(t) − u ref → 0,θ(t)− θ ref → 0, and v(t) → 0, as t → ∞ (proof in Sec. VI-C). Moreover, as we will see below, it requires only 4 tuning parameters. Starting from the surge u and yaw ω equations of dynam- ics in Eq. (5), the inputs F (t) and η(t), act through the body- frame components F (t) cosη(t) (surge) and F (t) sinη(t) (yaw). We split each actuation component into a virtual control, α u (t) and α ω (t), and a lumped term, d u (t) and d ω (t), capturing known/estimated disturbances in Eq. (5): F (t) cos(η(t)) = α u (t) + d u (t),(8a) F (t) sin(η(t)) = α ω (t) + d ω (t).(8b) To make those equations explicit, we define the lumped terms d u (t) and d ω (t) by collecting all non-actuation con- tributions from the surge and yaw equations (Eq. (5)), i.e., the projections of the boom tension and the hydrody- namic/coriolis terms. This yields d u (t) =− f l (t)·e u (t) + μ l u(t) u(t)− mω(t)v(t), (9a) d ω (t) = f l (t)·e v (t) + μ ω ω(t) ω(t) r .(9b) Substituting Eq. (8) and (9) into Eq. (5), we obtain two linear and one nonlinear EOMs, respectively: ̇u(t) = α u (t) m , ̇ω(t) = ̈ θ(t) = r I · α ω (t),(10a) ̇v(t) = −α ω (t)− μ t (v(t))v(t)− μ ω (ω(t))ω(t) r m − ω(t)u(t). (10b) Given the parameters α u ,α ω ,d u , and d ω we reconstruct the required propeller force and steering angle from Eq. (8a)- (8b), F =± p (α u + d u ) 2 + (α ω + d ω ) 2 , η = atan2 α ω + d ω F , α u + d u F , where the sign of F is selected to be positive, so that the steering angle stays within − π 2 ≤ η ≤ π 2 , ensuring the propeller continues to push from behind the hull. Now we can address the two linear systems (Eq. (10a)) with the virtual control terms α u and α ω as input, and we can design those terms to achieve a desired behavior. For the angular velocity ̇ θ, we choose to control the orientation, so the plant is of the form of a double integrator: ̈ θ(t) = γ ω α ω (t), γ ω = r I > 0. For the surge u, we choose to control the speed, so the plant is also of the form of a single integrator: ̇u(t) = γ u α u (t), γ u = 1 m > 0. In terms of control in the Laplace domain, denote by Θ(s) and U (s) the Laplace transforms of θ(t) and u(t), respectively. In addition, α ω (s) and α u (s) are the Laplace transforms of α ω (t) and α u (t), respectively. Hence, the above equations can be written as Θ(s) α ω (s) = γ ω s 2 (11) U (s) α u (s) = γ u s (12) Both plants in Equations (11), (12) can be controlled in a closed loop with a lead controller C(s) = K p · √ β· s + Ω c s + √ β· Ω c (13) Here, Ω c denotes the desired crossover frequency, set- ting the transient speed. We choose K p = Ω 2 c γ to set the desired crossover frequency, and set β > 1 to provide the required phase lead (phase margin). In addition, we implement this controller via a normalized feedback form, which preserves the stability and improves the transient response (See App. A). C. Theoretical guarantees of feedback-linearization con- troller We prove that the tracking error of the FBL controller tends to zero. Theorem 1 (Asymptotic stability). Consider the feedback- linearized closed-loop dynamics in Equations (10), with the virtual inputs α u ,α ω produced by the (normalized) lead controllers C u (s) = K p u √ β u s + Ω c u s + √ β u Ω c u , K p u = Ω c u γ u , γ u = 1 m , (14) C ω (s) = K p ω √ β ω s + Ω c ω s + √ β ω Ω c ω , K p ω = Ω 2 c ω γ ω , γ ω = r I . (15) Additionally, assume the following: (A1) the values u ref and θ ref are constants; (A2) the conditions Ω c u > 0, β u > 0, Ω c ω > 0, and β ω > 1 are met; (A3) the sway-drag satisfies μ t (v)≥ μ > 0 for all v ̸= 0; (A4) the boom-tension f l and drag terms μ l (u)· u, μ t (v)· v, μ ω (ω)· ω are bounded. Then, for any initial condition, the tracking errors satisfy u(t)−u ref → 0,θ(t)−θ ref → 0, and v(t)→ 0,as t→ ∞. Assumption A2 is solely based on the formation of a lead controller, and Assumption A3 is based on the non-negativity of the drag function μ t (v) as defined in Sec. V-A. We obtain the proof through the following claims. Claim 1 (Surge and yaw tracking are exponentially stable). Under Assumptions (A1)-(A2), the surge loop (u,α u ) and the yaw loop (θ,α ω ) are exponentially stable about the equilibrium induced by the constant references u ref and θ ref . Proof. We derive the characteristic polynomials for the surge and yaw loops and prove that they are Hurwitz under (A1) and (A2). This implies that the corresponding linear closed- loop dynamics are exponentially stable and track constant references with zero steady-state error. For the surge loop, with open-loop transfer function L u (s) = C u (s)P u (s), the closed-loop characteristic equation 1 + L u (s) = 0 yields χ u (s) = s 2 + 2 p β u Ω c u s + Ω 2 c u . For a quadratic p(s) = a 2 s 2 + a 1 s + a 0 , the Routh-Hurwitz criterion [28] reduces to a 2 > 0, a 1 > 0, a 0 > 0. Here a 2 = 1, a 1 = 2 √ β u Ω c u , and a 0 = Ω 2 c u , which are strictly positive for Ω c u > 0 and β u > 0; hence the surge closed- loop polynomial is Hurwitz [29]. For the yaw loop, with L ω (s) = C ω (s)P ω (s), the closed- loop equation 1 + L ω (s) = 0 gives the characteristic polynomial χ ω (s) = s 3 + p β ω Ω c ω s 2 + p β ω Ω 2 c ω s + Ω 3 c ω . For a cubic s 3 + a 2 s 2 + a 1 s + a 0 , the Routh–Hurwitz conditions are a 2 > 0, a 1 > 0, a 0 > 0, a 2 a 1 > a 0 . Here a 2 = p β ω Ω c ω , a 1 = p β ω Ω 2 c ω , a 0 = Ω 3 c ω , which are positive for Ω c ω > 0 and β ω > 0, and a 2 a 1 = β ω Ω 3 c ω > Ω 3 c ω = a 0 ⇐⇒ β ω > 1. Given that the conditions in Claim 1 are satisfied, we can now address the third DOF of the sway velocity v(t). Claim 2 (Sway stability under exponentially decaying dis- turbance). Consider the sway dynamics under feedback lin- earization: ̇v =− 1 m μ t (v)v + φ(t),(16) where φ(t) = −α ω (t) − μ ω (ω(t))ω(t)/r m − ω(t)u(t). (17) Assume: (S1) There exists μ> 0 such that μ t (v)≥ μfor all v ̸= 0; (S2) The function φ(t) decays exponentially: there exist ̄ φ > 0 and λ > 0 such that |φ(t)|≤ ̄ φe −λt , ∀t≥ 0.(18) Then v(t) → 0 as t → ∞. Moreover, v(t) converges exponentially. Proof. Choose the quadratic Lyapunov function V (v) = 1 2 mv 2 , which is positive definite and radially unbounded. Differen- tiating along Equation (16) yields ̇ V = mv ̇v =−μ t (v)v 2 + mvφ(t)≤−μ v 2 + m|v|·|φ(t)|. Using v 2 = 2 m V and |v| = q 2 m √ V yields the ISS-type bound ̇ V ≤− 2μ m V + √ 2m √ V |φ(t)|.(19) By standard ISS comparison arguments (e.g., Proposi- tion 4.19 in Khalil [27]), there exist class-KL and class-K functions β and γ such that |v(t)|≤ β |v(0)|,t + γ sup 0≤τ≤t |φ(τ )| , ∀t≥ 0. (20) Using Eq. (18), we have sup 0≤τ≤t |φ(τ )|≤ ̄ φ for all t, and sup 0≤τ≤t |φ(τ )| → 0 as t → ∞. Therefore, (20) implies |v(t)|→ 0. Furthermore, combining Equations (18) and (20), yields an explicit exponential estimate of the form |v(t)|≤ k 1 |v(0)|e −ct + k 2 ̄ φe −λt , k 1 ,k 2 ,c > 0, so v(t) converges exponentially to zero. Claims 1 and 2 imply u(t) → u ref , θ(t) → θ ref , and v(t)→ 0, completing the proof of Theorem 1. VII. EXPERIMENTAL EVALUATION In this section, we evaluate our routing (Sec. IV) and path- tracking (Sec. VI-B) approaches. 12510 0 2500 5000 7500 10000 12500 15000 17500 20000 Total Damage (lower is better) 25 spills 12510 0 20000 40000 60000 80000 100000 120000 Total Damage (lower is better) 50 spills 12510 Number of agents 0 25000 50000 75000 100000 125000 150000 175000 200000 Total Damage (lower is better) 100 spills Lower Bound Greedy Greedy + Optimal Subpaths + ILS MILP - no warm start MILP - with warm start Fig. 5: Final objective values for 25, 50, and 100 spills instances across varying numbers of agents. Shaded gray regions indicate lower bounds. A. Damage Minimizing Routing Solution Evaluation To assess the routing framework, we generate synthetic multi-spill scenarios over a bounded planar workspace W . Static polygonal obstacles are randomly generated to impose navigation constraints. The depot location is fixed, and spill locations are sampled uniformly at random within the free space of W . Each spill is assigned a risk weight R i drawn independently from a predefined range to represent heteroge- neous environmental severity. Transit distances d(S i ,S j ) are computed using occupancy-grid discretization of W followed by obstacle-avoiding A ∗ search between spill representatives (Sec. I-C), yielding a DMR instance. We evaluate two spill instances, with 25, 50, and 100 spills, and vary the number of ASV duos from one to ten. These scales reflect and exceed realistic large-scale spill fragmentation observed in major incidents [30]. For each configuration, experiments are run with a 300-second time limit on a laptop equipped with an Intel Core Ultra 9 185H CPU and 64 GB RAM. As a MILP solver, we rely on Gurobi [31], where configuration parameters are left at default settings. Solver Comparison. Results for a representative set of scenarios are given in Fig. 5. We evaluate four solution approaches: (i) Greedy corresponds to our assignment-only heuristic (H1) without solving a MILP; (i) Greedy + DP + ILS denotes the full heuristic (Sec. IV-B), again without invoking the MILP solver afterwards; (i) MILP denotes a pure MILP solution without using any warm-start heuristic; (iv) finally, MILP + warm-start initializes the MILP solver with the full heuristic solution (H1-3). Lower bounds cor- respond to the best dual bounds obtained by the MILP (iv) branch-and-bound solver within the 300-second time limit. For moderate agent counts (5 and 10 agents), the MILP + warm-start solver reaches optimal solutions that match the computed lower bounds, or near optimality for 100 spills. Interestingly, the combinatorial difficulty decreases as the number of agents increases. Conversely, for small agent counts, each agent must service more spills, increasing both assignment complexity and per-agent routing complexity. For 25-spill instances, (near-)optimal solutions are ob- tained for two agents. For the 25-spill single-agent case, the DP solver (H2) recovers the optimal solution directly within the allotted time, while the MILP lower bound did not fully converge. For more challenging configurations (e.g., 50 or 100 spills with 1–2 agents), the benefit of the warm- start heuristic becomes pronounced. Initializing MILP with the DP-augmented heuristic significantly improves solution quality within the time limit compared to MILP alone, with accumulated damage of routing without warm start exceed- ing twice the damage with warm-start in the worst case. Moreover, the standalone heuristic consistently produces high-quality solutions, often within a minute of computation. Our greedy priority queue heuristic (H1), consistently provides high-quality solutions in very short runtimes. The complete heuristic, applying single-agent DP and ILS re- finement, can improve solution quality by up to 15% with negligible additional runtime. The MILP solver, guided by this warm start, yet guarantees a small optimality gap in limited time, and in some cases improves the standalone heuristic solution by over 10%. Overall, this demonstrates that we can effectively deal with realistic problem instances. B. Controller Evaluation We evaluate the PID and FBL controllers (Sec. VI-B) for path tracking under varying conditions. Across the tested scenarios, both controllers produce stable tracking along the planned path—indicating that the FBL design, while requiring only half of the tuning parameters, captures the benefits of PID in this setup. That said, the PID consistently had an edge: it typically showed slightly lower steady-state errors. FBL matches PID ’s behavior well but pays a modest performance penalty. Controller setup. We tuned the PID gains using an iterative trial-and-error procedure. Starting from a baseline (low K p with K i = K d = 0), we increased K p until the response became sufficiently fast without sustaining oscillations. Next, we introduced K i to eliminate steady-state error, and finally added a small K d (with derivative filtering for causality) to reduce overshoot. We stopped tuning once further adjust- ments produced only marginal improvements in rise/settling Fig. 6: RMSE maps for cross-track (top) and heading errors (bottom) as functions of curvature radius and velocity for both vessels using feedback linearization control. time, overshoot, and steady-state error across the tested scenarios, while maintaining stable and consistent behavior. For the FBL controller, the parameters β u and β ω are se- lected to achieve a desired phase margin (through an explicit relation [32]) to improve stability. The crossover frequencies Ω c u and Ω c ω are chosen (through trial and error) to decrease the settling time of the transient response.Specifically, the parameter β ω is chosen to satisfy a phase margin of φ lead ω = 60 ◦ . In a stable plant of a double integrator (Eq. (11)), according to the Nyquist stability criterion [32], this is accomplished by the relation β ω = 1 + sin(φ lead ω ) 1− sin(φ lead ω ) .(21) For a plant of a single integrator (Eq. (12)), according to the Nyquist stability criterion, we already have a phase margin of 90 ◦ , thus we choose a phase margin of φ lead u = 20 ◦ , to slightly reduce overshoot. The value β u is derived similarly to Eq. (21). Simulation model and parameters. An evaluation was con- ducted in simulation using the coupled vessel-boom model in Eq. (5), where the boom enters as a measurable tow-point load f l . The boom is simulated as an articulated chain of n = 40 rigid links, with total length L = 40[m]. The vessel parameters were selected to approximate a small ASV and were assumed to be identical for both vessels. For ease of implementation, the boom is initialized as a straight segment of length L, between the two boats. To allow the boom to actually enclose and retain oil, the two towing vessels are commanded to follow offset paths whose separation is strictly smaller than L (see Fig. 2). If the inter-vessel distance were kept equal to the boom length, the boom would remain perpendicular to the direction of motion, effectively reducing to a taut straight line with negligible lateral containment; moreover, the resulting drag loads would be excessive and could practically damage the actuators. Specifically, the mass and yaw inertia were set to m = 600 [kg] and I = 500 [kg m 2 ], the propulsor offset (lever arm) to r = 2 [m], and the drag coefficients to κ ℓ = 100 [ N·s 2 m ], κ t = 10,000[ N·s 2 m ], and κ r = 1,000 [ Nm·s 2 rad ]. We assume calm-water conditions to isolate the track- ing performance attributable to the controllers themselves. The model neglects heave/roll/pitch, added-mass and wave- radiation effects, and spatially varying currents [26]; Thus, the reported results should be interpreted as a baseline comparison under a benign sea state, while incorporating wind/current/wave disturbances is left for future work. Test scenario. Each scenario consists of tracking a smooth reference path for an ASV duo. Both vessels start at the same initial pose and orientation, separated by a fixed distance. They are commanded to follow the same path with syn- chronized setpoints (see Sec. VI-A). The generated reference paths are Dubins paths [33] between (x 0 ,y 0 ,θ 0 ) = (0, 0, 0) and (x f ,y f ,θ f ) = (100, 65,π), as in Fig. 2. We also sweep the curvature radius ρ ∈ [10, 20][m] and reference speed v ref ∈ [5, 15][ m s ] to span from mild to aggressive maneuvers. The chosen initial and end points produce a trajectory long enough to exhibit both transient and steady-state behavior while remaining representative of local maneuvers used during the approach and encirclement of a spill. For the results, vessel 1 is the left vessel and vessel 2 is the right vessel. Results. For each (ρ, v ref ) pair, Fig. 6 reports the trajectory RMSE as our tracking metric. For cross track error (vessel 1), both controllers degrade with sharper turns (smaller ρ) and higher speeds. Perfor- mance is mainly curvature-driven. For FBL, this is apparent when the velocities are 12− 15 m/s and the curvature radius is 10− 13 m, and for PID with velocities 13− 14 m/s and curvature 10 m. Apart from the simulation where the turns are the sharpest and the velocities are high, PID outperforms FBL in curvature radii of up to 14 m (error of up to ≈ 6.5m vs an error of up to ≈ 7m). In the range of 15 − 16m, FBL outperforms PID (error of ≈ 7m vs ≈ 7.8m). From curvature radii of 17m, both controllers perform similarly with an error of ≈ 7.75m. Similar trends are observed for vessel 2, although with higher error values due to a longer reference path. For heading error, both controllers worsen at high speed and tight curvature. For FBL, the ”rough dynamics” region is substantially larger for both vessels (≈ 20 experiments), with high errors of up to ≈ 15 ◦ , while for PID, this region is much smaller (3 experiments) with high errors of up to ≈ 12 ◦ . When the dynamics are less demanding, FBL achieves heading errors of ≈ 5 ◦ for both vessels, while PID achieves heading errors of ≈ 3 ◦ for both vessels. Overall, both controllers incur larger cross-track and head- ing errors at higher speeds and tighter turns, reflecting the increased difficulty of tracking aggressive trajectories. These errors remain small relative to the spatial scale of typical oil slicks and containment operations. Theoretically, FBL guarantees zero steady-state error, un- like PID. In simulations, both tolerated different dynam- ics and completed all runs, but PID consistently achieved slightly better tracking. Practically, PID is harder to tune (eight gains per vessel), whereas FBL requires only four parameters. VIII. CONCLUSION AND FUTURE WORK In this work, we initiated the study of the problem of multiple oil-spill cleanup using boom-towing ASV duos. Our contribution lies in a careful modeling of the problem, and developing effective solution approaches for the correspond- ing routing and path tracking problems. In the future, we plan to validate the framework on real oil-spill data sets and evaluate our controllers under varying sea conditions (e.g., winds and currents), which we have neglected so far. From an algorithmic perspective, one issue that our current modeling of the routing problem overlooks is that once a spill is cleaned, other ASVs can pass through it, rather than around it, as we currently do. This might require a significant overhaul of our MILP approach, which we hope to explore in the future. Acknowledgments. The AI system ChatGPT was used for light editing and grammar enhancement, as well as a prelim- inary literature review. REFERENCES [1] P. F. Kingston, “Long-term environmental impact of oil spills,” Spill Science & Technology Bulletin, vol. 7, no. 1, p. 53–61, 2002. [2] L. Jr, M. Smith, and P. Ashcroft, “Analysis of environmental and economic damages from british petroleum’s deepwater horizon oil spill,” SSRN Electronic Journal, 2010. [3] I. Saadoun, “Impact of oil spills on marine life,” Emerging pollutants in the environment-current and further implications, 2015. [4] ITOPF, “Oil tanker spill statistics 2024,” ITOPF Annual Report, 2024. [5] Transportation Research Board and National Research Council, Oil in the Sea I: Inputs, Fates, and Effects. Washington, DC: The National Academies Press, 2003. [6] M. Fingas, Ed., Oil Spill Science & Technology, 2nd ed.Oxford, UK: Gulf Professional Publishing, 2016. [7] F. J. Pereda, H. G. De Marina, J. M. Giron-Sierra, and J. Jimenez, “Towards automatic oil spill confinement with autonomous marine surface vehicles,” in IEEE OCEANS, 2011. [8] A. Blum, P. Chalasani, D. Coppersmith, B. Pulleyblank, P. Raghavan, and M. Sudan, “The minimum latency problem,” in ACM Symposium on Theory of Computing, 1994, p. 163–171. [9] I. Leifer, W. Lehr, D. Simecek-Beatty, E. Bradley, R. Clark, P. Den- nison, Y. Hu, S. Matheson, C. Jones, B. Holt, M. Reif, D. Roberts, J. Svejkovsky, G. Swayze, and J. Wozencraft, “State of the art satellite and airborne marine oil spill remote sensing: Application to the BP Deepwater Horizon oil spill,” Remote Sensing of Environment, vol. 124, 2012. [10] S. V. Kumar, R. Jayaparvathy, and B. Priyanka, “Efficient path planning of auvs for container ship oil spill detection in coastal areas,” Ocean Engineering, vol. 217, p. 107932, 2020. [11] S. Dong, J. Feng, Z. Gu, K. Yin, and Y. Long, “A review of artificial intelligence and remote sensing for marine oil spill detection, classification, and thickness estimation,” Remote Sensing, 2025. [12] H. C. Kuzmenko, D. Ehevich, and O. Gal, “Autonomous oil spill response through liquid neural trajectory modeling and coordinated marine robotics,” Applied Ocean Research, 2026. [13] J. Giron-Sierra, A. Gheorghita, G. Angulo, and J. Jim ́ enez, “Preparing the automatic spill recovery by two unmanned boats towing a boom: Development with scale experiments,” Ocean Eng., 2015. [14] F. Arrichiello, H. Heidarsson, S. Chiaverini, and G. Sukhatme, “Co- operative caging using autonomous aquatic surface vehicles,” in IEEE International Conference on Robotics and Automation, 2010. [15] O. Elmakis and A. Degani, “USV port oil spill cleanup using hybrid multi-destination RL-CPP,” IEEE Access, 2023. [16] Y. Xu, L. Zhang, P. Zheng, G. Liu, and D. Zhao, “An emergency scheduling model for oil containment boom in dynamically changing marine oil spills: Integrating economic and ecological considerations,” Systems, vol. 13, no. 8, 2025. [17] L. Zhang, J. Lu, and Z. Yang, “Optimal scheduling of emergency resources for major maritime oil spills considering time-varying de- mand and transportation networks,” European Journal of Operational Research, vol. 293, no. 2, p. 529–546, 2021. [18] L. Zhang and J. Lu, “Optimizing oil spill emergency logistics: a time- varying multi-resource collaborative scheduling model,” Environmen- tal Science and Pollution Research, vol. 31, 2023. [19] E. DeCola, A. Hall, and M. Popovich, “Assessment of demulsifi- cation and separation technologies for use in offshore oil recovery operations,” U.S. Department of the Interior, Bureau of Safety and Environmental Enforcement, Tech. Rep., 2018. [20] S. L. Ross Environmental Research, LTD., “Transfer of decanting technology research to oil spill response organizations and regulators,” U.S. Minerals Management Service, Tech. Rep., 2005. [21] U.S. Environmental Protection Agency, “Development document for effluent guidelines and standards for the centralized waste treatment category,” U.S. Environmental Protection Agency, Tech. Rep., 2000. [22] A. E. F. Muritiba, T. O. Bonates, S. O. Da Silva, and M. Iori, “Branch- and-cut and iterated local search for the weighted k-traveling repair- man problem: an application to the maintenance of speed cameras,” Transportation Science, vol. 55, no. 1, p. 139–159, 2021. [23] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, Introduction to Algorithms, 3rd ed. Cambridge, MA, USA: The MIT Press, 2009. [24] M. Held and R. M. Karp, “A dynamic programming approach to sequencing problems,” Journal of the Society for Industrial and Applied mathematics, vol. 10, no. 1, p. 196–210, 1962. [25] H. R. Lourenc ̧o, O. C. Martin, and T. St ̈ utzle, “Iterated local search,” in Handbook of metaheuristics. Springer, 2003, p. 320–353. [26] T. I. Fossen, Handbook of Marine Craft Hydrodynamics and Motion Control. Wiley, 2011. [27] H. K. Khalil, Nonlinear Systems, 3rd ed. Upper Saddle River, N. J.: Prentice Hall, 2002. [28] E. J. Routh, A Treatise on the Stability of a Given State of Motion, Particularly Steady Motion.Macmillan and co., 1877, classical reference for the Routh stability criterion. [29] A. Hurwitz, “Ueber die bedingungen, unter welchen eine gleichung nur wurzeln mit negativen reellen theilen besitzt,” Mathematische Annalen, vol. 46, p. 273–284, 1895. [30] S. Sun, C. Hu, L. Feng, G. Swayze, J. Holmes, G. Graettinger, I. Macdonald, O. Garcia, and I. Leifer, “Oil slick morphology derived from aviris measurements of the deepwater horizon oil spill: Implica- tions for spatial resolution requirements of remote sensors,” Marine Pollution Bulletin, vol. 103, 2015. [31] L. Gurobi Optimization, “Gurobi Optimizer Reference Manual,” 2024. [32] K. Ogata, Modern Control Engineering, 5th ed. Upper Saddle River, NJ: Pearson, 2010. [33] L. E. Dubins, “On curves of minimal length with a constraint on average curvature, and with prescribed initial and terminal positions and tangents,” American Journal of Mathematics, vol. 79, no. 3, p. 497–516, 1957. APPENDIX A. Normalized feedback implementation We describe a normalized-feedback variant of FBL. The normalized lead controller from Eq. (13), wherein the static gain is equal to one, is placed in the feedback to reduce the overshoot of the system, while its static gain remains in the main branch. This is depicted in Fig. 7 and detailed below. This modification diminishes the effect of the lead- controller’s zero in the complementary sensitivity function. Furthermore, we apply a first-order low-pass filter to the reference signal to avoid discontinuous jumps in the control signals applied to the vessel. This is especially important when switching the reference from zero (when waiting for the other vessel to reach its setpoint) to a non-zero constant value (when both vessels can move again). Next, we explain how the revised normalized-feedback variant is derived, and show that, when coupled with the plants in Eq. (12) and Eq. (11), the closed-loop system is stable, while improving the the transient response. To make this more explicit, we compare the complemen- tary sensitivity function (CSF) for two architectures: (i) the standard series lead controller with unity feedback, and (i) the normalized lead placed in the feedback branch, with the static gain kept in the main branch (Fig. 7). The CSF maps the reference input Ref(t) to the closed-loop output y(t). In the Laplace domain, with Y (s) = Ly(t) and Ref(s) = LRef(t), the closed-loop transfer function is T (s) = Y (s) Ref(s) . The CSF T (s) provides direct insight into the closed- loop static tracking gain (T (0)), which indicates the ratio between a steady-state output and a constant reference input, and the effective bandwidth (how fast the system can track changes in the reference) [32]. By analyzing its structure for both architectures, we can show that the stability guarantees (Claim 1) extend to the normalized setting. For each architecture (i) and (i), we first establish stability, which is a prerequisite for discussing static gain and transient response. We keep the notation from the stability analysis. For surge we use (β u , Ω c u ,γ u ) with K p u = Ω c u γ u . For yaw we use (β ω , Ω c ω ,γ ω ) with K p ω = Ω 2 c ω γ ω . In the normalized topology, the main-branch gain is K p √ β , and the normalized lead H(s) = βs + √ βΩ c s + √ βΩ c is placed in the feedback path. Next, we consider the individual surge and yaw dynamics (respectively). Surge dynamics (single integrator plant). Consider the surge plant P u (s) = γ u s . With the standard series lead C u (s) = K p u √ β u s + Ω c u s + √ β u Ω c u , the complementary sensitivity becomes T std u (s) = Ω c u ( √ β u s + Ω c u ) s 2 + 2 √ β u Ω c u s + Ω 2 c u . With the normalized lead in feedback C main u (s) = K p u √ β u , H u (s) = β u s + √ β u Ω c u s + √ β u Ω c u , we obtain T norm u (s) = Ω c u √ β u s + Ω 2 c u s 2 + 2 √ β u Ω c u s + Ω 2 c u . Both complementary sensitivity configurations share the characteristic polynomial s 2 + 2 p β u Ω c u s + Ω 2 c u , hence the same stability condition (Claim 1). Moreover, T std u (0) = T norm u (0) = 1, so constant references are tracked without steady-state error. The only structural difference is the zero location: z std u =− Ω c u √ β u , z norm u =− p β u Ω c u . Since β u > 1, the normalized zero lies further left in the complex plane. Consequently, its influence on the dominant closed-loop poles is reduced, yielding smaller overshoot and less transient peaking. Yaw dynamics (double integrator plant). Consider the yaw plant P ω (s) = γ ω s 2 . With the standard series lead C ω (s) = K p ω √ β ω s + Ω c ω s + √ β ω Ω c ω , the complementary sensitivity is T std ω (s) = Ω 2 c ω ( √ β ω s + Ω c ω ) s 3 + √ β ω Ω c ω s 2 + √ β ω Ω 2 c ω s + Ω 3 c ω . With the normalized lead in feedback C main ω (s) = K p ω √ β ω , H ω (s) = β ω s + √ β ω Ω c ω s + √ β ω Ω c ω , we obtain T norm ω (s) = Ω 2 c ω √ β ω s + Ω 3 c ω s 3 + √ β ω Ω c ω s 2 + √ β ω Ω 2 c ω s + Ω 3 c ω . Piecewise constant Fig. 7: Lead controller in the normalized variant: the controller lies in the feedback while its static gain sits in the main branch along with the gain K P . The reference signal is a piecewise constant signal, filtered with a low-pass filter to avoid discontinuous jumps in the control signal. Here, again, both realizations share the same characteris- tic denominator and therefore the same stability condition (Claim 1), as well as unity DC gain: T ω (0) = 1. The difference lies exclusively in the zero: z std ω =− Ω c ω √ β ω , z norm ω =− p β ω Ω c ω . As in the surge case, the normalized implementation pushes the zero further left. This attenuates overshoot and reduces high-frequency amplification without altering the stability proofs of Sec. VI-C. This proves that the normalized feedback realization pre- serves the closed-loop characteristic polynomial and unity steady-state gain, while repositioning the controller zero deeper in the left half-plane. The result is improved transient behavior (reduced overshoot and peaking) without modifying the stability guarantees.