UNCERTAIN ENVIRONMENTS

RELATED APPLICATIONS

The present application claims the benefit under 35 U.S.C. §119(e) of U.S.

Provisional Patent Application Serial No. 62/689,011 entitled“OPTIMAL SHIP ROUTING IN STRONG, DYNAMIC, AND UNCERTAIN OCEAN CURRENTS AND WAVES,” filed June 22, 2018, which is incorporated herein by reference in its entirety.

FIELD OF INVENTION

The techniques described herein relate to the field of automatic determination of routes in dynamic and uncertain environments, and more particularly to stochastic techniques for determining routes that optimize one or more performance criteria.

BACKGROUND

Planning optimal routes for vehicles is important in many applications. For example, shipping routes can be optimized to reduce time and/or risk or optimized to reduce cost. Planning optimal routes can increase profits associated with a shipping business or ensure the saftety of a vehicle.

SUMMARY

Some aspects of the present application relate to a method for use in automatically determining an optimized route for a vehicle. The method includes, using at least one computer hardware processor to perform obtaining a target state, a fixed initial position of the vehicle, and dynamic flow information; and determining an optimized route from the fixed initial position to the target state using the dynamic flow information

Some aspects of the present application relate to at least one non-transitory computer- readable storage medium storing processor executable instructions that, when executed by at least one computer hardware processor, cause the at least one computer hardware processor to perform a method for use in automatically determining an optimized route for a vehicle. The method includes, using at least one computer hardware processor to perform obtaining a target state, a fixed initial position of the vehicle, and dynamic flow information; and determining an optimized route from the fixed initial position to the target state using the dynamic flow information

Some aspects of the present application relate to a A system that includes at least one computer hardware processor; and at least one non-transitory computer-readable storage medium storing processor executable instructions that, when executed by the at least one computer hardware processor, cause the at least one computer hardware processor to perform a method for use in automatically determining an optimized route for a vehicle. The method includes using at least one computer hardware processor to perform obtaining a target state, a fixed initial position of the vehicle, and dynamic flow information; and determining an optimized route from the fixed initial position to the target state using the dynamic flow information

The foregoing is a non-limiting summary of the invention, which is defined by the attached claims.

BRIEF DESCRIPTION OF DRAWINGS

Various aspects and embodiments of the disclosure provided herein are described below with reference to the following figures. It should be appreciated that the figures are not necessarily drawn to scale. Items appearing in multiple figures are indicated by the same or a similar reference number in all the figures in which they appear.

FIG. 1-1 is a block diagram of a route determination system, according to some embodiments.

FIG. 1-2 is a block diagram of a computer system, according to some embodiments.

FIG. 1-3 is a flowchart of a method of determining optimized routes, according to some embodiments.

FIG. 2-1 is a plot of a maximum speed profile produced by wind, according to some embodiments.

FIG. 2-2 is a plot of an optimal tacking path for a sailboat, according to some embodiments.

FIG. 2-3A is a map of a time evolution of optimal trajectories for two gliders after five days, according to some embodiments

FIG. 2-3B is a map of a time evolution of the optimal trajectories for two gliders after ten days, according to some embodiments FIG. 2-3C is a map of a time evolution of the optimal trajectories for two gliders after twenty days, according to some embodiments

FIG. 2-3D is a map of a time evolution of the optimal trajectories for two gliders after twenty-five days, according to some embodiments

FIG. 2-4A is a magnified view of the time evolution of the optimal trajectory for the first glider, according to some embodiments, and corresponds to the region in the lower dashed box of FIG. 2-3 B.

FIG. 2-4B is a magnified view of the time evolution of the optimal trajectory for the second glider, according to some embodiments, and corresponds to the region in the upper dashed box of FIG. 2-3 B.

FIG. 2-4A is a plot of minimum arrival-times for the first glider, according to some embodiments.

FIG. 2-5B is a plot of minimum arrival-times for the second glider, according to some embodiments.

FIG. 3-1 is a schematic of a setup for determining a stochastic time-optimal path, according to some embodiments.

FIG. 3-2A is a plot of the mean of a coefficient for a dynamically orthogonal decomposition of a velocity field, according to some embodiments.

FIG. 3-2B is a plot of the mode of a coefficient for a dynamically orthogonal decomposition of a velocity field, according to some embodiments.

FIG. 3-2C is a plot of the probability density function (PDF) of a coefficient for a dynamically orthogonal (DO) decomposition of a velocity field, according to some embodiments.

FIG. 3-3 is a cumulative histogram of relative error in arrival-time determined between the DO decomposition approach and a Monte Carlo (MC) approach, according to some embodiments.

FIG. 3-4A is a Frechet distance (normalized) between reachability fronts computed using the MC approach and the DO approach at time T=25, according to some embodiments.

FIG. 3-4B is a Frechet distance (normalized) between reachability fronts computed using the MC approach and the DO approach at time T=50, according to some embodiments.

FIG. 3-4C is a Frechet distance (normalized) between reachability fronts computed using the MC approach and the DO approach at time T=75, according to some embodiments. FIG. 3-4D is a Frechet distance (normalized) between reachability fronts computed using the MC approach and the DO approach at time T=100, according to some

embodiments.

FIG. 3-5A is a plot of the spatial distribution of reachability fronts at time t=0.05, according to some embodiments.

FIG. 3-5B is a plot of the spatial distribution of reachability fronts at time t=25, according to some embodiments.

FIG. 3-5C is a plot of the spatial distribution of reachability fronts at time t=50, according to some embodiments.

FIG. 3-5D is a plot of the spatial distribution of reachability fronts at time t=75, according to some embodiments.

FIG. 3-5E is a plot of the spatial distribution of reachability fronts at time t=l00, according to some embodiments.

FIG. 3-5F is a plot of the spatial distribution of reachability fronts at time t=l25, according to some embodiments.

FIG. 3-6 is a plot of stochastic time-optimal paths for three targets, according to some embodiments.

FIG. 3-7 includes plots showing the mean field, variance of the DO coefficients, the DO mode fields for five mode fields, and the marginal PDF of the corresponding DO coefficient of a stochastic double-gyre flow field at a beginning time, t=0 days, and an end time, t=l3.5 days, according to some embodiments.

FIG. 3-8 includes plots showing the stochastic flow field for a first realization with a negative first coefficient and a second realization with a positive first coefficient at a beginning time, t=0 days, and an end time, t=13.5 days, according to some embodiments.

FIG. 3-9 includes plots showing stochastic reachability fronts over time, according to some embodiments.

FIG. 3-10A is a plot of the stochastic time-optimal paths colored with the velocity DO coefficient 1, according to some embodiments.

FIG. 3-10B is a plot of the stochastic time-optimal paths colored with the arrival time at the target, according to some embodiments.

FIG. 3-11 A includes plots of the DO mean flow (row A) and the DO mean flow for mode 1 and mode 2 fields (row B), according to some embodiments. FIG. 3-11B includes plots of the marginal PDFs of coefficients 1 and 2 (row C), according to some embodiments..

FIG. 3-11C includes plots of the variance of the first eight modes (row D), according to some embodiments.

FIG. 3-12 includes plots of two realization of the stochastic flow fields at different times, according to some embodiments.

FIG. 3-13 includes plots of stochastic reachability fronts overlaid on the velocity streamlines of mode 1 (column A) and mode 2 (column B).

FIG. 3-14 includes plots of the stochastic time-optimal paths colored by velocity coefficient 1 (column A) and coefficient 2 (column B) to target 2, target 5, and four other targets, according to some embodiments.

FIG. 3-15 includes plots of arrival time distributions at each of six targets, according to some embodiments.

FIG. 4-1 is a schematic of a minimum-risk time-optimal path planning setup, according to some embodiments.

FIG. 4-2A is a plot of the domain of flow strength for a stochastic simulated front crossing, according to some embodiments.

FIG. 4-2B is a plot of the PDF flow strength for a stochastic simulated front crossing, according to some embodiments.

FIG. 4-3A is a plot of stochastic reachability fronts at time T=2, according to some embodiments.

FIG. 4-3B is a plot of stochastic reachability fronts at time T=2.5, according to some embodiments.

FIG. 4-3C is a plot of stochastic reachability fronts at time T=3.5, according to some embodiments.

FIG. 4-3D is a plot of stochastic reachability fronts at time T=6, according to some embodiments.

FIG. 4-3E is a plot of stochastic reachability fronts at time T=9, according to some embodiments.

FIG. 4-3F is a plot of a time-optimal path distribution for a vehicle navigating in a stochastic steady front with uncertain flow strength, according to some embodiments. FIG. 4-4 includes plots of cost matrices (row 1), risks associated with each path (row 2), and the risk-optimal paths (row 3) for risk-seeking behavior (column a), risk-neutral behavior (column b), and risk- averse behavior (column c), according to some embodiments.

FIG. 4-5A is a plot of risk-optimal waypoint objectives for a stochastic steady front crossing, according to some embodiments.

FIG. 4-5B is a plot of PDF error due to following a risk-optimal path for a stochastic steady front crossing, according to some embodiments.

FIG. 4-6A is a plot of the Frechet distance between the time-optimal path and the risk- optimal, risk-seeking path, according to some embodiments.

FIG. 4-6B is a plot of the Frechet distance between the time-optimal path and the risk- optimal, risk-neutral path, according to some embodiments.

FIG. 4-6C is a plot of the Frechet distance between the time-optimal path and the risk- optimal, risk-averse path, according to some embodiments.

FIG. 4-7A is a plot of risk-optimal heading objectives for a stochastic steady front crossing, according to some embodiments.

FIG. 4-7B is a plot of the PDF of errors due to following the risk-optimal heading objectives for a stochastic steady front crossing, according to some embodiments.

FIG. 4-8A is a plot of a realized risk-seeking path corresponding to a particular flow realization and is colored by the discrete Frechet distance between this path and the true- optimal path for the realized environmental flow, according to some embodiments.

FIG. 4-8B is a plot of a realized risk-neutral path corresponding to a particular flow realization and is colored by the discrete Frechet distance between this path and the true- optimal path for the realized environmental flow, according to some embodiments.

FIG. 4-8C is a plot of a realized risk-averse path corresponding to a particular flow realization and is colored by the discrete Frechet distance between this path and the true- optimal path for the realized environmental flow, according to some embodiments.

FIG. 4-9 includes plots of the mean (first row), standard deviation (second row), skewness (third row), perturbation from the DO mean using a first realization (fourth row), and perturbation from the DO mean using a second realization (fifth row) for a stochastic wind-driven double-gyre flow at times T=0 days (first column), T=6.75 days (second column), and T=13.5 days (third column), according to some embodiments. FIG. 4-10 includes plots of cost matrices (row 1), risks associated with each path (row 2), and the risk-optimal paths (row 3) for risk-seeking behavior (column a), risk-neutral behavior (column b), and risk- averse behavior (column c), according to some embodiments.

FIG. 4-11 A is a plot of risk-optimal waypoint objectives for a stochastic double gyre flow field, according to some embodiments.

FIG. 4-11B is a plot of PDF error due to following a risk-optimal path in a stochastic double gyre flow field, according to some embodiments.

FIG. 4-12A is a plot of the Frechet distance between the time-optimal path and the risk-optimal, risk-seeking path, according to some embodiments.

FIG. 4-12B is a plot of the Frechet distance between the time-optimal path and the risk-optimal, risk-neutral path, according to some embodiments.

FIG. 4-12C is a plot of the Frechet distance between the time-optimal path and the risk-optimal, risk-averse path, according to some embodiments.

FIG. 4-13A is a plot of risk-optimal heading objectives for a stochastic double gyre flow field, according to some embodiments.

FIG. 4-13B is a plot of the PDF of errors due to following the risk-optimal heading objectives for a stochastic double gyre flow field, according to some embodiments.

FIG. 4- 14 A is a plot of a realized risk-seeking path corresponding to a particular flow realization and is colored by the discrete Frechet distance between this path and the true- optimal path for the realized environmental flow, according to some embodiments.

FIG. 4-14B is a plot of a realized risk- neutral path corresponding to a particular flow realization and is colored by the discrete Frechet distance between this path and the true- optimal path for the realized environmental flow, according to some embodiments

FIG. 4-14C is a plot of a realized risk-averse path corresponding to a particular flow realization and is colored by the discrete Frechet distance between this path and the true- optimal path for the realized environmental flow, according to some embodiments

FIG. 4-15A includes plots of the mean (first row), mode 1 (second row), mode 2 (third row), and mode 3 (fourth row) for a stochastic flow exiting a strait at times T=0 minutes (first column), T= 50 minutes (second column), and T= 100 minutes (third column), according to some embodiments.

FIG. 4-15B includes plots of the marginal PDF for the coefficient of mode 1 (first row), the marginal PDF for the coefficient of mode 2 (second row), and the marginal PDF for the coefficient of mode 3 (third row) for a stochastic flow exiting a strait at times T=0 minutes (first column), T= 50 minutes (second column), and T= 100 minutes (third column), according to some embodiments.

FIG. 4-16 includes plots of the standard deviation (first row), the skewness (second row), the kurtosis (third row), a first representative realization (fourth row), and a second representative realization (fifth row) of a velocity field for a stochastic flow exiting a strait, according to some embodiments.

FIG. 4-17 includes plots of stochastic reachability fronts colored by DO velocity coefficient #1 at different times, according to some embodiments.

FIG. 4-18 includes plots of time-optimal path distributions colored by arrival time, according to some embodiments.

FIG. 4-19 includes plots of risk-optimal paths (for risk-seeking, risk-neutral, and risk- averse choices) with waypoint objectives in a stochastic flow exiting a strait for five different target locations (column 1) and corresponding PDF of errors due to following the risk- optimal paths (column 2), according to some embodiments.

FIG. 4-20 includes plots of risk-optimal heading objectives for risk-seeking, risk- neutral, and risk-averse choices in a stochastic flow exiting a strait for five different target locations (column 1) and corresponding PDF of errors due to following the risk-optimal heading objectives (column 2), according to some embodiments.

FIG. 4-21 includes plots of trajectories obtained by following risk-optimal headings (for risk-seeking (first column), risk-neutral (second column), and risk-averse (third column) choices) with waypoint objectives in a stochastic flow exiting a strait for five different target locations (rows a-e), according to some embodiments.

DETAILED DESCRIPTION

Determining optimized paths has applications for exploratory, surveillance, and imaging missions using EiABs and/or UAVs as well as in determining optimized shipping routes for container ships. In many application, water and/or air flows encountered by the vehicles may be comparable in magnitude to the nominal speeds of the vehicles.

The inventors have recognized and appreciated that rigorous and efficient techniques for determining optimized routes are needed. Any number of performance criteria may be optimized, such as travel time, energy consumed, or risk encountered.

The inventors have recognized and appreciated that a deductive forward computation of optimized trajectories of vehicles operating in dynamic and uncertain flows may be determined with a computational efficiency that is one or two orders of magnitude better than conventional techniques. In some embodiments, a time-optimized route to a target state is determined by starting with a fixed initial position and determining a forward reachable set of fronts using an unsteady Hamilton- Jacobi (HJ) equation.

The inventors have further recognized and appreciated combining decision theory with stochastic path-planning can result in a new partial differential equation-based scheme for risk optimal route planning in uncertain and dynamic flows. By combining a principled risk optimality criterion grounded in decision theory with stochastic dynamicallt orthogonal level-set equations, an efficient computational scheme to predict optimized paths from a distribution of stochastic optimized paths is determined.

The inventors have further recognized and appreciated that accounting for uncertainty in optimized route path planning results is important for many applications. To efficiently solve the partial differential equations involved, some embodiments use dynamically orthogonal reduced-order projections that result in several orders of magnitude in

computational speed-up relative to Monte Carlo techniques.

FIG. 1 is a schematic diagram of a route determination system 1-100, according to some embodiments. The route determination system 1-100 includes a computer server 1-101, a computer 1-103, a data storage device 1-105, a first vehicle 1-107, and a second vehicle 1- 109, which may be communicatively coupled together by a network 1-110. The network 1- 110 comprises one or more networking devices for transmitting information from one point of the network 1-110 to another. The network 1-110 may include a local area network (LAN), a wide area network (WAN), and/or the internet. The network 1-110 may include connections, such as wired links, wireless communications links, and/or fiber optic cables. The network 1-110 may include wireless access points, switches, routers, gateways, and/or other networking equipment as well as any suitable wired and/or wireless communication medium or media for exchanging data between two or more computers, including the Internet. The wireless connections may be implemented using radio signals, optical communication signals, and/or satellite links. Computer server 1-101, computer 1-103, and data storage device 1-105 are illustrated as connected to network 1-110 with a wired connection (either electrical or optical). The first vehicle 1-107 and a second vehicle 1-109 are illustrated as being wirelessly connected to the network 1-110. However, embodiments are not so limited. For example, any of the components may be connected with a wired connection or a wireless connection.

As depicted, the first vehicle 1-107 is illustrated as a shipping vessel and a second vehicle 1-109 is illustrated as a submarine. In some embodiments, the vehicles may be other watercraft such as tankers, bulk carriers, container vessels, passenger vessels, autonomous underwater vehicles (AUVs), sailboats, underwater gliders, or yachts. In other embodiments, the vehicle need not be watercraft. For example, in some embodiments, the vehicles may include aircraft, such as airplanes, helicopters, drones, or unmanned aerial vehicles (UAVs). In some embodiments, the vehicles may be manned or unmanned. In some embodiments, the vehicles may be autonomous or manually controlled by a human operator.

The data storage device 1-105 may be one or more storage devices, such as hard drives, tape drives, optical drives, and other suitable types of devices. The data storage device 1-105 may be located in a single location or may be distributed in different locations. In some embodiments, the data storage device 1-105 may provide data and other information to the computer server 1-101, the computer 1-103, the first vehicle 1-107, or the second vehicle 1- 109. For example, the data storage device 1-105 may contain information about the environment associated with one or more vehicles, such as dynamic flow information associated with the air or water. For example, the dynamic flow information may include one or more flow fields.

FIG. 1-2 is block diagram of a computing device 2-100 according to some embodiments. The computer server 1-101 or the computer 1-103 may have one or more of the components described in connection with the computing device 2-100. Additionally, the first vehicle 1-107 and/or the second vehicle 1-109 may include a computer with one or more of the components described in connection with the computing device 2-100.

Computing device 2-100 may include at least one computer hardware processor 1- 210, a memory 1-220, a non-volatile storage 1-230, an input/output (I/O) device 1-240, a network adapter 1-250, and/or a display 1-260. Computing device 1-200 may be, for example, a desktop or laptop personal computer, a personal digital assistant (PDA), a smart mobile phone, a tablet computer, a computer server, or any other suitable computing device.

Network adapter 1-250 may be any suitable hardware and/or software to enable the computing device 1-200 to communicate wired and/or wirelessly with any other suitable computing device over any suitable computing network, such as network 1-100. The network adapter 1-250 may be used to obtain information from the data storage device 1-105. For example, dynamic flow information and flow uncertainty information may be obtained from the data storage device 1-105.

The display 1-260 may be any suitable display for displaying to a user a visual representation of the optimized route determination results. For example, the display 1-260 may be a computer monitor, a LCD display, or a touchscreen display. The results displayed may include a map, a list of headings, and/or instructions for following the determined route.

The non-volatile storage 1-230 may be adapted to store data to be processed and/or instructions to be executed by processor 602. For example, the non-volatile storage 1-230 may include at least one non-transitory computer-readable storage medium storing processor executable instructions that, when executed by at least one computer hardware processor, such as the processor 1-210, cause the at least one computer hardware processor to perform a method for use in automatically determining an optimized route for a vehicle.

The memory 1-220 may be a volatile memory device such as random- access memory (RAM) that may be controlled by the processor 1-210.

Computer hardware processor 1-210 enables processing of data and execution of instructions. The processor 1-210 may cause computer executable instructions stored in the non-volatile storage 1-230 to be loaded into the memory 1-220. The processor 1-210 may then the instructions to perform a method for use in automatically determining an optimized route for a vehicle. The data and instructions stored on the non-volatile storage 1-230 may comprise computer-executable instructions implementing techniques which operate according to the principles described herein.

While not illustrated in FIG. 1-2, a computing device may additionally have one or more components and peripherals, including input and output devices. These devices can be used, among other things, to present a user interface. Examples of output devices that can be used to provide a user interface include printers or display screens for visual presentation of output and speakers or other sound generating devices for audible presentation of output. Examples of input devices that can be used for a user interface include keyboards, and pointing devices, such as mice, touch pads, and digitizing tablets. As another example, a computing device may receive input information through speech recognition or in other audible format.

FIG. 1-3 is a flowchart of a method 1-300 for determining an optimized route for a vehicle, according to some embodiments. The method 1-300 may be implemented using all or a portion of the route determination system 1-100. For example, a computer on the first vehicle 1-107 or the second vehicle 1-109 may perform all of the actions in the method 1- 300. Alternatively, the server 1-101 and/or the computer 1-103 may perform a portion of the method 1-300, while other actions of the method 1-300 may be executed by the first vehicle 1-107 or the second vehicle 1-109. For example, the server 1-101 may determine an optimized route for the first vehicle 1-109 using flow information from data storage device 1- 105 and then send information detailing the optimized route to the first vehicle 1-109 where a computer on the first vehicle 1-109 may control the first vehicle 1-109 to implement the optimized route.

At act 1-301, the method includes obtaining a target state, a fixed initial position, and dynamic flow information. In some embodiments, one or more computers performing the method 1-300 may perform the act of obtaining 1-301. For example, any of the computers or servers described above may perform the act of obtaining. In some embodiments, the target state includes a desired destination for the vehicle. In some embodiments, the fixed initial position may include the current position of the vehicle. For example, radar or GPS may be used to determine the initial position of the vehicle. The dynamic flow information may include, for example, a flow field. The dynamic flow information may be obtained from a third party, such as a weather service. The dynamic flow information may also include uncertainty information pertaining to the certainty that the flow field is accurate. The dynamic flow information may, for example, be stored on the data storage device 1-105 and obtained via the network 1-110.

At act 1-303, the method 1-300 includes determining an optimized route from the fixed initial position to the target state using the dynamic flow information. In some embodiments, the act of determining 1-303 may include any number of actions described in detail in Section II, Section III, and/or Section IV, below.

At act 1-305, the method 1-300 includes displaying one or more indication of the optimized route. This act 1-305 may not be performed in every embodiment. In some embodiments, the indications of the optimized route may include a list of speeds and/or heading to be implemented by a human operator of the vehicle. In other embodiments, the indications may include a map of the vehicle and its surrounding.

At act 1-307, the method 1-300 includes automatically controlling the heading of the vehicle to follow the optimized route. This act 1-307 may not be performed in every embodiment. In some embodiments, the act 1-307 is performed for an autonomous, unmanned vehicle. In other embodiments, the act 1-307 may be performed even when there is a human operator of the vehicle if, for example, the human operator is operating the vehicle in an autopilot mode.

Some aspects of the technology described herein may be understood further based on the disclosure and illustrative non-limiting embodiments provided below in Sections II, III, and IV. Each of Section II, III, and IV are self-contained such that equations referenced in each section only refer to the equation in the respective section. Similarly, sub-sections and appendices referenced in each section refer only to the sub- sections and appendices associated with the respective section. Similarly, references cited in each section refer only to the references associated with the respective section. II. Minimum-Time Path Planning In Strong Dynamic Flows Using A Forward Reachability Equation

1, Introduction. This paper provides a theoretical synthesis for time optimal control of anisotropic robotic units operating under the influence of strong and dynamic flow -fields. The deduction is based on a formulation of forward-reachability in terms of the viscosity solution to an unsteady Hamilton- Jacobi (ILl ) equation. This problem is of great interest in several areas of engineering. For example, it· is applicable to Unmanned Aerial Vehicles (UAVs), for exploratory, surveillance or imaging missions, transportation and delivery, as well as environmental and climate research.

It is also pertinent to Autonomous Underwater Vehicles (AUVs) such as ocean gliders and surface crafts, for diverse missions including ocean mapping and sampling, naval reconnaissance and harbor protection. Flows encountered by those vehicles can often be comparable in magnitude to their nominal speeds (speed relative to the flows). For both air and sea vehicles, it is imperative to develop rigorous and efficient techniques for optimal routing, minimizing either the energy consumed, the travel time, or other performance criteria [24] Additionally, certain anisotropic vehicles” are characterized by direction-dependent motion constraints, which must be accounted for. In this paper, we focus on the deductive forward computation of minimum-time trajectories of anisotropic vehicles operating in complex time and space dependent flows.

Minimum-time problems have been well studied in the control theory and operations research literature. The first results are likely due to blalton who coined and employed Isochrones for time-optimal ship routing problems [18, 391. Since then, the most common approach to solve such problems uses the Dynamic Programming Principle (DPP) to derive an unsteady Hamilton -Jacobi— Bellman (HJB) equation for the globally optimal‘time to go ^’ . which is the minimum time required to drive the system to a ftsced target state [2, 10). The viscosity solution )13| of the governing

HJB equation provides an optimal feedback control law and the resulting optimal system trajectory. Several methods have been proposed to solve thi HJB equation. One of them is the extremal field approach [41, 42] , in which a family of initial value Euler Lagrange equations are solved to determine extremal (locally optimal j trajectories. Alternatively, one can compute the time-evolving‘backward-reachable’ set (also known as capture basin), the set- of points from which there exist controls that drive the system to the target state in a specific time |3, lb] .

The computation of backward-reachable sets is a fundamental tool in the safety verification of embedded and hybrid control systems. It enables one to ascertain the set of states from which a system may be driven (perhaps inadvertently) to potentially dangerous configurations. For instance, in air traffic control, an unsafe configuration occurs when the distance between any two aircrafts is lesser than a required mini mum |37, 35] Similar examples occur in oceanic applications, wherein the operating agents are AUVs and/or ships [14, IS] As a result, the design of efficient optimal control methods to compute the backward reachable set is receiving much attention. Specifically the DPP can he used to arrive at an unsteady HJ equation that governs the evolution of the backward reachable set, as explored in |35, 33j and references therein. It Is important to note that this HJ equation generally does not hold in a classical sense, i.e., the value function is not differentiable everywhere. In fact, the optimal time-to-go may not even he continuous everywhere. The formal link between backward-reachability and the solution to the HJ equation is established through the framework of viscosity solutions, which may be computed numerically using level - set methods [38| for low -dimensional problems and the more recent semi-Lagrangian schemes |8] for higher dimensional ones. The unsteady HJ equation offers an expedient and robust characterization oi the backward -reach able set, even for systems with complex non-autonomous dynamics [(3j and state constraints (i.e., requirement that system trajectories remain inside a given set at all times) 171, We refer to [34, [ _< 2.3] and 133: for detailed discussions of other algorithms for computing the backward-reachable sets in a variety oi applications and dynamics, both discrete and continuous.

The optimal time-to-go is equivalently characterized by a static HJ equation, when the system is locally controllable jib, 3, 50, Si. In our application, the local controllability property holds only when the speed of the external flow is lower than that of the operating vehicle. Moreover in this case, the optimal time-to-go is a uniformly continuous function [3] Similarly, when the external flow is steady, the time-derivative term of the unsteady HJ equation drops out and we arrive at a static HJ equation for the optimal time-to-go. In both of these cases. Fast Marching Meth ods 144], Ordered Upwind Methods (OUMs) [45, 431 or Fast Sweeping Methods [23] can be then used (see aim [501 ] For other references to path- lanning in dynamic Mows (not necessarily optimal), the reader is referred to [29, 42]

A complementary view-point to backward reachability and control is the concept of "forward reachability’. The forward reachable set is the set of all states to which a system can be driven at a given time, starting from a fixed initial state (recall that in the definition of backward reachability, it is the target state that is fixed). While not extensively studied-as it leads to an open-loop control law-, forward reachability has potential advantages in several applications. Examples include predictive pa.th - phmning and safe guidance of vehicles operating in dynamic fluid environments using data- assimilative forecasts of the environmental fields. In such examples, which are a main emphasis of this paper, i is critical to accurately forecast the regions that can be explored or reached by vehicles in a given time, prior to deployment. In ocean applications, as gliders surface only once every few hours to communicate, their continuous feedback control is not strictly possible. Moreover, predictions of currents (which directly affect AIJVs) are periodically updated by assimilating sensor observations and therefore, optimal AUY trajectories need to be frequently recomputed to account for such updates. Similar comments can be made for the planning of specific air vehicles (e.g, drones) based on data-assimilative wind forecasts. Finally, an open-loop control that is quick to recompute serves equally well as a feedback control. Due to all of the above, we focus on forward reachability in this paper. cations, as gliders surface only once every few hours to communicate, their continuous feedback control is not strictly possible. Moreover, predictions of currents (which di rectly affect AUVs) are periodically updated by assimilating sensor observations and therefore, optimal AUV trajectories need to be .frequently recomputed to account for such updates. Similar comments can be made for the planning of specific air vehicles e.g, drones) based on data.-assimilative wind forecasts. Finally, an open-loop control that is quick to recompute serves equally well as a feedback control. Due to all of the above we focus on forward reachability in this paper.

We determine the time-optimal control by computing the optimal‘arrival timed which is the minimum time required by a system to reach a target state, starting fro a fixed initial state. Using various concepts rooted in the theory of non -smooth calculus we demonstrate that the forward reachable set of anisotropic vehicles Is governed by the viscosity solution of an unsteady ELI equation, even in cases when the arrival time is discontinuous. This result provides a foundation for extending the use of numerical tools such as the level-set method to solve for the forward reachable set and the optimal trajectories. For completeness, we also show that when the system is locally controllable, the arrival time is continuous and satisfies a static If J equation. All the above results are the main contributions of this paper, providing a synthesis for forward minimum-time path-planning in dynamic flows, An exposition of the corresponding methodology for strictly isotropic vehicles and its application to a variety of ocean-specific path-planning scenarios is provided in [SI, 30]

The paper is structured as follows. §2 formally defines the minimum arrival time problem an the notation used. The main theoretical results and discussions are presented in §3, In §4, these results are used to illustrate properties of minimum-time trajectories and arrival time fields for a sailboat operating in a uniform wind field, and autonomous gliders in the Sulu Archipelago, where the flow-field estimates are obtained from a realistic daia-assimilative ocean modeling system. Conclusions are listed in §5,

2. Problem Statement, Let Q be an open set and U be the set of unit vectors « in K". Let F (&,£) : U x {0, oo) ® be a given bounded time- dependent function. Consider a vehicle {P} moving in il under the influence of a dynamic (non- homogeneous, unsteady) flow-field, V(x, i) : W x [th oe) --- M ^s , We wish to predict a control for P that mmiimzes its travel tmie between given start and end points, denoted by y _s and y _f respectively Let a general continuous trajectory from y _s to yg he denoted as Xp(t; y _s ). Any vehicle’s motion, being the sum of nominal motion due to steering and adveetion due to the flow-field, is governed by the kinematic relation

(2.1 ) ^ = U (Xp(i;y _g ), i) = F _P (h(i)fr) h(i) + V (X _P (f; y _s ) , where ( h(i),£j is the cfrrection--dependent speed of the vehicle relative to the flow, satisfying 0 < Fp | hit), t ) h h(i) is the unit vector in any vehicle’s

heading (steering) direction at time i; and U(Xp(:i; y _s ), i) is the total vehicle velocity. Let T(y) : W—* £ denote the‘first arrival time ⁵ function for a given Fp, h(i); i.e., T(y) is the first time the vehicle reaches any point y . starting from y _s for the given Fp. Clearly, T(y _« ) = 0, The limiting conditions on Xp(i; y _s ) are

We aim to predict the optimal controls for h (t) and Fp that minimize

T(y _f ) subject to constraints (2.1) and (2.2). Let the optimal travel time to reach f be T*( yr) and the corresponding optimal trajectory be X _f> (t; y _s ). We assume that V (x. t) is exactly known. In realistic ocean applications, forecast flow-fields are always associated with some levels of uncertainty 126, 271. V (x, £) can correspond to, for exam le, the mode or mean of the predicted flow-field, up to the predictability limit. We consider cases where the distance traveled by the vehicle is much larger than its dimensions and thereby assume the interaction between the vehicle and the How-field to be purely kinematic. The notation * j in this paper denotes the L-2 norm of *. We assume that F ^ h(f),t^ and h(i) are Lipsehitz continuous and that V (x, i ) is bounded and Lipsdiitz continuous in both x and t, i.e. d C, Cy > 0 such that

and

(2.4) Ov (|xi— C Ϊ 4- life

The forward reachable set 'R j: y _s of points at time i > 0 is defined as the set of all points y e ί I such that there exists a trajectory Cr(t; y _g ) satisfying (2.1), with Xp(0; y _s ) = y _s and Kp{t y _s ) = y. The boundar of this set is called the re&ehabiiUy front and is denoted by y _s ). Recall that the subset of trajectories Xp($; y*) that reach _f was denoted by X (f;y _s ). Just as p(i; y _S ) represents any feasible trajectory starting at y _s , T{y ) denotes the first arrival time (not optimal) at any point y. We refer to Appendix A for relevant definitions and theorems that will be used iu the following section, which synthesizes the main theoretical results. Next, we also omit the time dependency on several quantifies where no confusion is expected. 3 Main Theoretical Results. Let f be the viscosity solution of the HJ equa tion: with initial conditions

(3.2) _< i(x _f 0)=x(x),

where v : W ^. is Lipselutx continuous, and the boundary conditions on f are open

(see Appendi B). In the following lemma, we prove a monotoriicity result related to f. We show that the generalized gradient of f is non- positive on trajectories

Xpif; y _s ), along the direction 1 ) for t > 0. This lemma is then used in

a theorem that establishes the relationship between forward reachable sets and the viscosity solution of an unsteady HJ equation, which can also be viewed as a modified level-set equation. For the definitions, lemmas and theorems used in the proofs here, we refer to Appendix A.

when f ⁹ denotes the gener

Proof The function

Lipsehitz in all three arguments and convex in the third argument, by Claim 1. The viscosity solution to (3.1) is unique, by Remark 3 and locally Lipsdiitz, by Lemma G.

Let us fix £ > 0 and let (q p ⁾ € &F ( Xp(i; y _s ), t . We first derive a useful inequality. From Theorem 4 , Consequently,

where the first inequality follows from i j Since (3.5) must

From CoroEarv 1,

where O _cj C ll (0, oo) is the set of points at which ø is not differentiable. We now show by contradiction that the .RHS of (3.7) is noiwpositive.

Since f is locally Lipsehitsq there exist e _x , _t > 0 such that for all (xq i ¹ } f Q _K * and j '— X (t;y _s )| < e _x , W— t\ < <¾ the derivatives’vd(xh fi) and are bounded in magnitude. Let. us assume that the RHS of (3.7) equals k > 0. Then, there exists a sequence ((xf, b _wi loi ^® ^Xp{*i * ) suc that

By the BolKano-Weierstruss theorem, the sequence {(x), i;() J^ admits a subsequence {(x . i ) f _{; i} such that the sequence converges. Let this

limi be (q ^w _s ^f ). From (3.8), Taking limits as i oo yields

by construction. Therefore

from (3,7),

Since the choice of i > 0 was arbitrary, we obtain (3.3),□

Proof (I Let > 0, From Remark 3, the

viscosity solution to (3,10) is unique and by Lemma· 6. it is locally Lipsehitz. We now argue that d _p is locally Lipsehitz. Observe that :=

^Xp(£;y _s },if. Since p{t) is continuously differentiable in (0, oo)

-fp- 1 ] and f° is locally Lipsehitz, <j> _p l) is also locally Lipsehitz in (0, oo) by the rule stated in Theorem 6,

Let ti > 0 be fixed. Since f _r is loc-ally Lipsehitz, there exists an open interval around in which f° _r is Lipsehitz. Thus, for any > ti in this interval, Lebourg’s Mean Value Theorem (Theorem 5) implies that there exist £ ₃ € and s 6 d<p _p (t _$ ) such that

(3.14} f {ί ₂ ) ( ) = s a (t _{2 ~} ti) ,

Using the chain rule of Theorem & again {* denotes the adjoint),

$fr(%) C (g (i ₃ )}" &F (gp(fo))

Hence, for any s e. ΰf _> (f _$ ), there exists (q, ) such that

implying that for any a V cfo _p ifo), a < 0. Using this result in (314) yields _r {1 ₂ ) < f _r {H) for all t ¾ Since y _r is locally Lipsehitz in (0, 00}, we conclude that d (f) is non-increasing an (0, oc). Moreover, since f _r is continuous on

(from (3J1)} and non-increasing in (ϋ,oo), we have $ {t) f ^a ( p{i:y _s } _;i | < 0 for all t > 0.

(2), It as shown in part (1) that f° all t > 0, for any trajectory

Xp(f;y _s ) that satisfies {2.1} and the initial conditions X (0;y _e ) = y _s . Therefore, for a trajectory Xp|i;y _s ) that reaches a given end point y€ W at time T(y) (not necessarily optimal),

Since this inequality holds for any arbitrary arrival time T{ ), it will also hold for the optimal arrival time T°i y). implying

(3.18) ^ ^® (y,T°(y)) < 0 for all y€ W .

For y = y _S (3.12) holds trivially. For any y y y _s , fr ^o (y, 0) > 0 by (3.11). The continuity of y and (3.18) together then yield

(3.19) T°(y) > inf (t : r(y, t) = 0} .

.}>0

Final]} , any point on the zero level-set of f ⁰ belongs to a characteristics of (3,10) emanating from y _s , since y _s is the only point in W where f ⁹ is initially zero. There fore, when the zero level-set reaches y for the first time, it implies the existence of a trajectory X|»(t; y _s ) with Xf» U; y _s ) = y _s that satisfies (21). Dor this trajectory,

(3.19) holds with an equality, thereby establishing (3.12). Physically, this means that the fastest arrival time at any end point y S W occurs when the zero level-set of f° reaches y for the first time, and equivalently that the reachability front (boundary of the forward reachable set) coincides with the zero level-set of f ^> at time t.

(3). Under the assumption mu¾ F iihfr) > sup _x€$ * f 1 V (x, £} I b the start point y _s belongs to the interior of the forward reachable set R.{i; y _s ) for all t > 0, i.e, for airy t > 0, there exists e _t > 0 such that all points x' satisfying jy _s— xfi < _t are members of fo(i; y _s ). This is equivalent to the "Small Time Local Controllability ^? condition discussed in [31 as a· result of which, T ^® is continuous in SI. See [31 for a formal proof of this statement Let us fix y G SI. By definition T°{y) satisfies

(3.20) T°{y) = inf (T°(y) + k] ,

where y y il is a point such that there exists a trajectory Xp((;y _s ) satisfying (2.1) and the limiting conditions

(3-21) X _P (T°(y); y _s ) = y, Xp(T°(y) y hi y _s ) = . In order to show that T° is a viscosity solution to (3.13), we show that it is both a viscosity suhsolution and a supersohition to (3,13).

Viscosity Sufoselutioan Rro Definition 4 and Remark 1, T° C(O) is a viscosity subsolution to (3.13) if at every y 4 W arid for every f? ¹ function r _e : S --- K such that c _s (y ^’ ) = T ^i> (y) and r _s ^. T° has a local minima at y,

Since T _S > T° in neighborhood of y, we obtain for h > 0 small enough.

⁷ - _* (y) My) r°(y) ^Te y)

Moreover tor this choice of h and the resulting y, (‘d.'20) implies

T°(y) £ T ^u (y ) + h

Combining the above two inequalities yields

(3,23) r _* (y) - My) < T y) - Tfoy) < h ,

Since t _« is differentiable at y, Taylor's theorem may be used to expand r ₉ {y) near y,

to rfotsfo I V

Vr, (y) · ^(T°(y); y _fi ) = F (h, . T°(y)) h* Vfo(y } - V{y _s T°( )) - Vr* (y) < 1 , thereby establishing {3,22). Therefore, T° is a viscosity subsolution to (3.13).

Viscosity Superso!ution. T° is a· viscosity supersolution to (3.13) if at any y€ W and for every C ¹ function r* : W— » K such that r*(y) = T°(y) and T T° has a local maxima at y.

(3 26) max { f ( h, T°( y)) h · Vr ⁴ (y)} 4 V(y, T°(y)} VT ^S ( ) - 1 > 0 .

For any 0 < h < °(y), there exists y C-I O satisfying T°(y) 4 h = T°(y ) and a trajectory X (t;y ₈ ) satisfying (2.1 ) and the limiting conditions

Of course, the optimal trajectory leading to y is a valid choice for Xp(i; y _s ). For h > 0 small enough,

(3.28) h = T°( y) ^. T°(y) < T*(y) - r ^s (y) .

As in the earlier section, WB may use Taylor’s theorem to expand T*(y) near y to obtain

As in the earlier section, we may use Taylors theorem to expand r* ^' (y) near y to obtain

(

Inserting (3.20) in (3.28) and dividing by h,

Observe that from (2.1 ),

width roves that T ^® a visc sity onporsoiotion of (3 13) Therefo e, T ^® is a viscosity Solution to (3.13).

THEOREM 2, Let the assumptions of Theorem 1 hold.

1. In addition, if f is regular at ( (t; y _s ), t) for all t > 0 ami Xf _> satisfies

Proo (1), When the trajectory X _p (t;y _si is regular, he. when #' is regular at points (.X _p (i; y _s ), t) for all t > 0, Lemma 4 implies t) =

¾6° X _p (b y _s ), £) for all t > 0. Since ° is the viscosity solution to (310), we obtain from Theorem 4 that for any

for all {«¾,£>}€ d-^ C _r (ί; y _g ), t) = tifCKffit y _s ) _f i). Specifically, for the member (q ^s ,]R) of fth ⁰ (X (h _s ),tS _> tha satisfies (3.31), the definition of generalized gradient (A 9) implies

(333) = 0.

Combining this result with (3.3), we obtain

Since y° is regular at {XifU: y _s ), t) by assumption, we then also have For any h > 0, the definition implies

(3.35)

Since f ^u is locally Lipsehiiz, d C > 0 such that for h > 0 small enough,

(3.36)

where o(/< _· }€ I denotes a vector whose individual components are £>{ }. Adding and subtracting -f h ) from the numerator of (3 35) a d

using the triangle inequality, we obtain

From (A.?), the first term on the right converges

as h j 0 and by (3.34), its value is aero, The second term uniformly converges to aero as h I 0, by definition. This implies

and consequently that (3.37)

Since (3.37) holds for all i > 0, f _r is right differentiable in {0, oc) and the value of fire right-derivative is zero for all f > 0. This implies tha f _> is constant in (0, oo). Since d* _p (G} = 0, we obtain > _p (t) = f ^a (X (f ; y ₈ ), i) = 0 for all t > 0. Therefore, trajectories X (£; y _s ) that are regular and satisfy (3.31) always remain on the zero level- set of h (2). For the second part of the proof, we fix y* £ il. From Theorem 1, the opti mal trajectory to y _f satisfies °(X _p (i; y _s ), f) = 0 for all t > 0. Hence, $ _p {t) := i) equals zero for all D < t < T*(y _f ). Let us fix a time t such that f ^i> is differentiable at y„), f). The usual chain rule then yields

(3.38)

where the derivatives of f° are evaluated at (X* (i; y _s }, (). Since y ^s is assumed to he differentiable at this point, (3.10) holds in the classical sense and

Substituting this in (3.38) gives

equality holding iff F = F and h* . Using

this result in (3.39) yields

which completes the proof, P

3.1 Remarks.

Relationship between forward reachability and ": Theorem 1 lays down the math ematical foundation for computing forward reachable sets of vehicles operating in dynamic flow-fields, for any given anisotropic, time -dependent maximum relative speed function F As shown the forward reachable set coincides with the

set of non positive i.e.. for any t > C),

(3 40) U{t _s ) = (y€ W : o'ly . i) < 0} . where f° is the continuous viscosity solution of (3.10). Therefore, the problem of forward reachability reduces to that of accurately computing the viscosity solution to (3.10). This unsteady formulation also has several merits over the stationary one given by (3.13),

Advantages of the unsteady formulation: As mentioned in «1, T° can become discontinuous when jV(x, f)i when local controllability does not hold.

Accurate treatment of such discontinuities requires expensive and often ambiguous grid refinements. All of this is avoided by the unsteady formulation (3.10); by con struction, it governs a continuous f ^a . Note that (3.10) contains the extra dimension of time t when compared to (3.13), This higher dimensional embedding is what allows us to continuously track the forward reachable set (and reachability front). This may seem to add cost, but powerful solvers for (3.10) can be obtained by leveraging schemes (e.g. level set solvers) dedicated to such unsteady 113 equations and allowing subgrid-scale accuracy. Moreover, one is typically not as interested in tire actual value of as much as in its sign. In other words, since TZit y^) depends only on the sign of ffofx, t), it suffices to track the zero level set of h Narrow -band level set methods _< an be used for this purpose, achieving significant reduction in computa tional cost [31], when compared to the regular level set method

unit vector normal to the contours of <f>P.

Novelty: Our approach to reachability differs from common techniques used to derive Hi equations for backward reachable sets [7, 33, 16|. Such schemes first characterize the target set using a scalar function, and consider a control problem that optimizes this function for various admissible trajectories, Using the DPP and a constructive definition of the value function, they show that it satisfies an unsteady HJ equation. Presently, using non - smooth calculus, we showed that all admissible trajectories lie in the region where is non-positive. This result was then used to show, in a deductive fashion that d°(y3) < 0 is equivalent to the forward reachable set. We note that this was obtained without invoking convex reachability [36] . 4. Applications. We now illustrate our minimum-time path— lanning results using two examples. The first example Is that of a sailboat navigating in still water, downstream of a uniform wind field. Here, the maximum relative speed of the sailboat (F) is anisotropic, since it depends on the angle of the sail .relative to the wind. In the second example, we use our methodology to compute the time-optimal trajectories of two isotropic autonomous underwater gliders operating in complex and realistic ocean flows of the Sulu Archipelago in the Philippines.

4.1. Minimum-time tacking of a sailboat in a uniform wind. Mobile agents such as sailboats cannot be steered at the same speed in all directions. They rely on wind drag for thrust, and their speeds vary with the magnitude of the wind and angle of the sail. Anisotropy is also exhibited by surface crafts, kayaks and AUVs that encounter wave effects. This anisotropy poses an additional challenge for the computation of minimum-time trajectories, with many racing and historical ramifications 118! . Through this example, we illustrate how the methodology developed in this work can be utilised to perform time optimal path planning, rigorously accounting for such anisotropic motion constraints. Specifically, we consider a sailboatoperating in a steady and spatially uniform wind field and compute its time optimal trajectory to a point located directly downwind from its starting location.

Fig. 2-1 depicts a polar diagram of the sailboat's maximum relative speed (F) produced by wind blowing in the horizontal direction.

Fig. 2-1 : Maximum speed profile (h, t, i) produced by wind blowing in the positive-x direction Ϊ.

This idealized speed profile F that we employ represents the effect of wind drag on the sailboat and is given by the analytical expression

(4.1) (fo b: li fe)

for h such that 0 < $(h;T ) < Here. $(h; ίc„, ( : U ® (0, 2 K ) denotes the angle between It and the wind direction 1 . For h such that TT < #{h; h _a .) < 2?r, we replace d(h; h,,,) by 2w— #(h; h _H, ) in (4.1) so as obtain a symmetric polar diagram (see Fig 2-1 }. We see that for this speed profile, the sailboat cannot be steered directly upwind, and that the maximum speed is achieved when the sailboat is steered at an angle of approximately 60° to the wind. This analytical form for F is chosen so that tire resulting; polar diagram qualitatively matches the speed profile of a typical sailboat |40j. Nonetheless, the path planning methodology is also valid for other sailboats with varied speed profiles and those operating in unsteady and spatially variable winds. The sailboat stmts at a point with coordinates y _s = (0 2.0.5) and must navigate to the point yr (0 8, 0.2) in minimum time. These points are depicted in Fig. 2-2 as a shaded circle and star, respectively. The sailboat experiences a steady and uniform wind field in the positive-x direction, i.e., h _¾ = i. The minimum-time trajectory of the sailboat is determined using a two-step numerical algorithm, extending |32, 311 to anisotropic motions. In the first step, the forward reachable set 'Ri J y _s ) is computed by numerically solving the unsteady Hi equation (3.10), using a narrow-band level-set method, 7?,(i; y _s ) is evolved until time

Since F does not vary with time, the optimization in (3,10) is performed offline, prior to the start of the tune-marching loop. We note that F(h i;h _w )h nf° =

{ Fi t; h _¾; )h - h jV i where n = is the unit vector in the direction of Y °

Prior to the time-marching loop, we populate a lookup table of the optimal values of (h, t; h _w )h h (optimization done over ίi) for various choices of n. At each time- step, the optimal \aluo f F(h t] h _m )h fi at every grid point is simply read-off from the lookup t&bh . hatod the unit vector n at that grid point, The boun ary of Tb(dys (ie, dRi _* 11 «xfr&cfce from the scrnfor Held F ⁹ nod stored for later use. The sec n step of the algorithm toes deter ine' fht pt l tfoeeoto by s lving the ordinary ditTetentoil o t nation photo backyar hi titar. uning n P w s > Recall that 13.32) k velid only m. points a hee k" differenncshLe tone* ;M k kwallv Lipsehliz (from Lemma fig Is k itif tiabh' hm* eu r>whi re b B.nh^n *rbt rV theorem. Thereto», the Leeks nokmg equati n ¾2 fto ab _{<* %} nlid al ost even when·. At trajectory points win re to s not dish > variable. he smsnerhvd ,» io iu i k»n to f yields one of the ener lis'd sp.uhd gra ient _* . oik to W A con t queue _*, ,»r.t one of the optimal trajectorie from lm fondly p ^' L. ^' H i k o toeved, The repaired optimal heading directions e e oouquutvl sing tfo dkcieie contours fofoe ¾\ _l e t were snood during the first step, W- e h..Mre wo n- '

F that does not. vary with tone, the algorithm h valid < wn In the ease -svlmre f ^' k time-dependent, The unumrka! toiunc vex to solv «-. s ,3 iff» and hboto asv summarized in Appendix B .

The optimal trajectory of the sailboat is shown in Fig 2-2.

Fig.2-2: Optimal sailboat tacking paxli superimposed on intermediate zero level -sets of f°.

Also shown are the intermediate reachability fronts &RAi: y _s ). uniformly spaced in time. We observe that since the sailboat can move faster sideways of the wind direction, the intermediate reachability fronts are elongated normal to the wind direction. Secondly, the speeds in the upwind direction (negative-x) are much lower than those downwind. As a consequence, the reachability fronts are more closely spaced in the upwind direction. An optimal trajectory to the end point is the inverted‘V ⁵ shape; the sailboat tacks at an optimal angle to the wind for the first half of its journey and at the symmetrically opposite angle for the second hall This numerical result agrees well with intuition and basic sailboat racing [e.g. the sailboat can be navigated faster at an angle to the wind than downwind). Of course, for boats that can be steered fastest downwind, the time-optimal trajectory is the straight line connecting the start and end points. Such tests have been completed [221. but are not shown here. A key feature ofonr path- planning methodology is that it applies to any physical speed profiles and wind-fields, time-dependent and inhomogeneous. In ail cases, it determines the time-optimal trajectory in a rigorous but also computationally efficient maimer.

4.2. Minimum-t ime pat h— planning in the Siiht Archipelago. We now illustrate minimum-time path-planning results for the case of two autonomous gliders operating in the Sulu Archipelago in the Philippines. The geometry is complex, with numerous islands and passages of varied sizes. The ocean flows are also complex and of a multi-scale nature, due to the geometry and internal dynamics, but also due to the variable forcing from the larger-scale open ocean, tides, and atmosphere. Therefore, the flow fields are highly unsteady, with many dynamic features such as strong tidal currents, intermittent eddies, coastal currents and larger-scale streams.

In tills example, we consider two underwater gliders that are deployed in the Sulu sea at coordinates 6.84° N 120.71° E and 7.61° N 121.51° E. These two start points are depicted by full black circles on Fig, 3-Fig 5, They are to be driven, in minimum time, to a common target end point in the Sulawesi sea, with coordinates ye. 4,26° N 120,42° E, shown as a full black star. The gliders are assumed to be isotropic and capable of speeds up to 0.25 m/s (T 0.5 knots) relative to the flow, i.e., (h,f) = 0.25 m/s. Since F is isotropic, the term max _} -

(3.10) reduces to F \Vm \. Hence, in this case, no optimizations are required at each time-step to solve (3.10)

The multi-scale ocean flow dynamics are modeled using the MIT Multidisciplinary Simulation, Estimation and Assimilation System (MSEAS) (21]. The flow-field estimates are obtained from data-assimilative simulations during February 5-March 24, 2009 as a part of the Philippine Straits Dynamics Experiment (PhilEx, 119]). Details regarding the numerical schemes used to solve the governing ocean equations may be found in [28 21] As they execute their mission, the gliders are assumed to follow a. yo-yo pattern in the vertical (from the near surface to 400 in depth or near the seabed if shallower), and the horizontal flow they encounter during this process is integrated over the yo-yo pattern. The resultant horizontal velocity is assumed to force the simulated gliders, i.e. this horizontal velocity is V (x, i). For a spatial discretization of 3 k x 3 km resolution and a time step of 5 min are employed. Fig. 2-3 depicts the time-evolution of the optimal trajectories of the two gliders, overlade on the flow-field, colored by its magnitude (cm/s).

Fig. 2-3: Time-evolution of the optimal trajectories of two gliders released in the Sul si sea (start locations numbered 1 and 2, and marked by full black circles). The target end state in the Sulawesi sea is marked by a full black star. Glider positions at intermediate times are grey shaded circles, ^ Streamlines of the ocean currents are overlaid on a colored map of the flow magnitude (cm/s ^' ) .

Intermediate glider positions are show.» as grey shaded circles. In parte! (a), we .see that Glider 2 encounters a stron cyclonic gyre, frequent features hi the Suite sea. As a result if is forced to move northward away from the end point. Glider I, in contrast, does not ex rieiK : this e dy aid is able to directly head south. Along the way, it encounters strong tidal flows as it traverses through passages between islands, as Shown in panel (b) and in the corresponding magnified view in Fig. 2-4a.

0 Fig. 2-4: Magnified views of optimal glider trajectories, 10 days after their release. These magnified areas correspond to the regions marked in Fig. 2-3b.

These strong and oscillatory flows cause the glider to take a zig zag baek-and-forth path as it crosses the islands. Glider 2. initially forced offtrack by the strong eddy, turns and resumes its journey towards the end point after roughly 10 days of deployment, completing a rem rkable loop in the process (Fig. 3b, Fig. 4b). Panels (c) and (d) of Fig. 3 display the remaining portions of the optimal trajectories of both gliders. Overall, Glider 1 takes roughly 15 days to reach the end point, while Glider 2, though deployed at the same time, arrives at the target in 25 days. These are striking differences in optimal paths and shortest arrival times, The results also show that due to the strong currents (e.g, tides and gyres), certain points along the optimal trajectory may be visited multiple times, indicating discontinuities in the field of the optimal arrival time T°(y).

Fig. 2-5 shows the fields of optimal arrival times T°{y) for each of the two deployment locations.

Fig, 2-S: Minimum arrival-time fields T”(v } . corresponding to the deployment location of each glider.

Such maps are generated by recording the first time each point in the domain becomes a member of the forward reachable set y _s ). Various ridges in Figs 2 -5a and 2-5b dearly indicate the discontinuities in T Ay). This example, therefore, illustrates the applicability of the unsteady formulation (3 10) for forward reachability and the ability to implicitly handle strong flows, without requiring expensive and often ambiguous grid refinements near the boundary of the reachable set. S. Conclusions. We presented a theoretical synthesis of forward reachability for minimum-time control of anisotropic autonomous vehicles navigating in dynamic flow -fields A moti vation of this work is to optimally operate such vehicles by utilizing the numerical predictions of environmental flow fields. The approach is based on exactly evolving a forward reachable set, which is the set of all states to which the system can be driven, at a given time, An open-loop control is derived for the minimum-time trajectories, taking into consideration the infrequent communication between the vehicles and the central controller. We showed that the forward reachable set is governed by the viscosity solution of an unsteady Hamilton— Jacobi equation. This equation is shown to be valid even when environmental Sow speeds exceed the maximum relative speed of the vehicles, Le„ when local controllability does not hold. The unsteady H.J formulation is shown to be equivalent to a static one under the restriction of local controllability. The algorithm is applied to determine the minimum-time trajectory of a sailboat moving in a uniform wind-field and the optimal arrival-time fields of two autonomous underwater gliders deployed in the Suiu sea with target end points .in the Sulawesi sea. These applications demonstrate the usability of the algorithm in realistic scenarios, highlighting the general validity of the unsteady HJ formulation, including overcoming the limiations of the static one.

Appendix A. Preliminaries.

in this section, we describe some of the relevant definitions and other theoretical details needed for the main results. Most of the material presented in this section may be found in [3, 12 _» 11, 17, §h hi what follows, we let. n€ fib ί l C M ^¾ be an open set and n ® M.

(A.2) \H(x t q) -- Af(x ₂ ,f ₂ ,q)j < C _s \q\ x (fx _j --- x ₂ j 4- :q -- t ₂ \)

The proof of this claim uses a straight-forward application of the triangle and Cauehy- ScliwarK inequalities, and is omitted.

DEFINITION 1. Let Xo e W. The set of super-differentials of x at Xo is defined as

f x(c) C _. [xn} q · (x xp) ^ \

(A.3) ύ ₊ x{Cc) := q e S” : limsup

c c¾ |X— Xol

Similarly. the set of sub- differentials of f at x _§ is

(A.5) #÷C(xo) = 0-€( ©) = {Vfixo}} -

2. If ti+ iXf) and d-fix Q) are both non-empty then x is differentiable at Q and (A.5) kokh.

DEFINITION 2. Let x he locally IApsckitz at o. The generalized directional derivative x ^® (co; u) : R ^w ® K, far an u€ M ⁿ is defined os

xί ^c + feu) ---- £{x)

(A.6) iim sup— ^: - - - h

hiii

The directional derivative, defined ns i fe

The set of generalized gradients of x at x¾ is the nan-empty set

(A .12) <¾^(x ₀ ) U fi-C xs) C i¾(x _f} ) .

REMARK 3, ( Uniqueness. } The Cauchy problem {3,1} - (3.2) admits at most me bounded , uniformly mnUnuom msmsitg solution (sec Corollary 5.5 of (9j).

LEMMA tl. Let y be the unique viscosity solution to {3.1} and (3.2). Then f is locally Lipsckitz, s.e. for every (XQ, ¾) f (0, oo), them exists a neighborhood of (cb, ί o ) in -which F is Lipschiiz continuous. The proof of this lemma uses convexity of H(x,t, q) = ma¾ {F (u, f) u - q} + V(x, t) · q in the third argument and may be found in [48] (Proposition 1.18} or [4] (Proposition 3.20) or [11] (Theorem 7 2.3) .

THEOREM 4. A locally Lipsckitz function f : Q x (0, oo) - ¾ is a viscosity sahMcm to {3.1} iff for ever x, t)€ ii x (U, oe),

The proof of this theorem may he found in [17|.

(A.19) /(¾) - f(x) = g x (y - r) .

for some e df{z ) , where- s = Xx ÷ (1 ^. A)y

THEOREM: 6. ( Chain Rule.) Lei Q C and :¾ C h two open sets with m. n it N. Lei g : if ® S¾ be continuously diffemr&iable near x. y ifr , and let F : ffr ® M be LipschUz near g{x). Then / := F o g s Lipsckiiz near x and

(A, 20) #/(x) c (g f i)) ^* dF (gfx)) , where ^* denotes the adjoint The proofs of Theorems 5 and 0 may be found in 11 ] . Appendix B. Numerical schemes.

We now summarize the numerical schemes utilized to discretize and solve (3.10) and (3 32) . These equations are solved using a Finite Volume framework. The term ¥f ^> in (3,10) is discretized using either a first order |29f or a higher-order upwind scheme [511 ; V (x _t i ) - V^ ⁵ is discretized using a second-order Total Variation Diminishing scheme on a staggered C-grid f4§| .

Evolution of the forwani reachable set We discretize [3.10) in time as follows:

(B.I) is solved only in the interior cells of the domain. For boundaries that are open inlets outlets or side walls (Le. not interior obstacles), a zero normal flux can be used on f?. We have also used weaker open boundary conditions such as zero second- order derivatives, so as to allow free advection of non-zero gradients, e.g [47, 4tij. In general, more complex open boundary conditions can be used, as done in regionalocean modeling {25. 20] . Obstacle in the domain are masked. For points adjacent to the mask, the necessary spatial gradients are evaluated using neighboring no es that do not lie under the mask,

The reachability frost d ^' R is extracted from the narrow band f° held at every time step using a contour algorithm. We also note that instead of this extraction, we Backi cking ; (3,32) is discretized using first order 129] or higher-order j;51] time integration schemes. Ideally, it suffices to solve (3.10) until the level set front first reaches y _# . However, due to the discrete time steps, a more convenient stopping criterion is the first time, T, when , ] < 0 Due to this, yj does not lie on the final contour <¾T(T; y _s , i _s ) exactly, Thus, we first project y _f onto d ^" R T y _s , i _s ). Starting from this projected point, the discretized form of (3.32),

(B.2)

is inarched backward in time until we reach a point on the first saved contour and this process generates a discrete representation of X _p (t;y ₈ ). Along the way, we project each newly computed trajectory point, X. _p {t— Ai; y _s ) onto the corresponding intermediate level-set contour (see |29]}. Par F that does not vary with time, the optimal heading hd f i) in (B.2) is read off from the lookup table (described in §4 1) based on the direction of the gradient $° (X (i;y _s ), i).

Further algorithmic and implementation details may be found in |29j.

REFERENCES

it] DAVID A DALSTKINSSON AND JA ES A . SETHJAN, A fast level set method for propagating Interfaces, Journal of { omputationai Physics, 118 (1995 , pp. 261* - "277.

i2|i M. ATHANS AND P.L. RAΪ,B, Optimal Control: L» Ini dmtion to the’Fheory And Its Apph- c&iiorm, Dover Books OK Engineerin Series. Dover Publications, 2006.

|3| M- BA.RD1 AM> I, CAPD KZO-DOLCETTA. Optimal Control and Visaosdy Solutions of Eamiltim- Jaeobi-Melknun Equations Modem Birkhauser Classics, Birkiiauaer Boston, 2008.

4j; S. BΪLNOHΪ AND B. TOR ON, SBV Regularity for Hamilton- Jacobi equations with Hamiltonian depending an SIAM Journal on Mathematical Analysis, 44 (2012), pp. 2178-2203.

IS;! A.M. BLOCH, Nonhohmormc Mechanics and Control, Interdisc: ipihi ary Applied Mathematics,

Springer, 2003,

16:1 O, BOKANOWSK!, A. BBΪANE AND H, ZIDAM, Minimum time control problems for non- autonomous differ tial quation . Systems L· Control Letters, 58 (2000), pp. 742 - 746.

!7j O, BOKANOWSKL i» . Po.RCADEL, AND H. ZTDANL Roachahihly and minimal times for stale constrained r.onltnear problem s >.vdhoui any controtiahilif y assumption, SIAM Journal of Optirna] UontroL 48 {20101, pp. 4292-4316.

!¾] OLiVIER BOKANOCOiKI, JOCKEN GARCKE. MiCHAEL CiUEBEL, AND IRENE KLOMPMAKER, An adaptive sparse grid simii-hignmgian scheme far first order hamiUon-jacoM bellman equa tions, Journal of Scientific Computing. 55 12013). pp. 575-60 ^” .

|9| A'JWiRTO BfiEsSSAN, Viscosity solutions of Hamilton- Jacahi Equations and Optimal Control Problems, An Illustrated Tutorial, Penn State University, Spring 2011.

fiOl A . li, B YSON AND Y . CL HO, Applied Optimal Control: Optimisation, Estimation and Control, John Wiley and Sons, 1975

fll| P. C ^' AXNARSiA AND C. StffiSTRA&l, Sem xmcave Functions, Hamilton Jacobi Equations, one (optimal Control, Progress in NonBnear Differential Equations and Their Applications, Birkhauser Boston, 2004. [121 F. H. CLARKE, YLT. S. LEDYAEV, R. 3. STEB.N. AND P . R WOLENSKJ, Nonsiocvth analysis anti control th ry, Springer- Verl&g New York, Inc., Sec&ueus, J. USA, 1998.

[ 13] M G. CRANDALL A D P. L. LIONS, Viscosity solutions o h milicm-jacobi aq ation·! TVansao fir s of the American Mathematical Society, 2 ' 77 i 1983). pp. 1 43.

[14] T. B. CURTIN, D. M. CBIMMIN-S. -P Cuncio. M. BENHMIN, AND C. ROPER, Antonamims underwater vektdcs: Trends and transformations, Marine Technology Society Journal, 33 ( 2065-09-01 T 00 : 00 : 00 } pp.65-75.

[IS! R, E, DAVIS, N. E. LEONARD, AND D. M. FRATANTONL Routing strategies for underwater gliders, Deep Sea R&seareh Part II: Topical Studies in Oceanography, 56 (2009), pp. 173 - 187.

[16] M. FALCONE, T. GiORGi, AND P. LORETI, Lewd sets g f wsscosHy ^oh insnw some applb at»ms to fronts and rendezvous problems, SLAM Journal oil Applied Mathematics, 51 < 1934) '

dients. Journal of Mathematical Analysis and .Applications, 141 ( 1989), pp 21 - 2€.

[18! FRANCIS GALTON. OR ike employment of meieoroiegicsl statistics in deirn mining the best course for a ship whose sailing qualities are known, Proceedings of the Royal Society of London, 21 (1872), pp. 263 274 ^'

[19] ARNOLD I.. G ORDON A D CESAR L. VILLANOY Oceanography. Special issue an the Philippine Stmits Dynamics Experiment, vol. 24 The Oceanography Society, 2011.

[20! P, J. HALEY JR. AND F. F. J, LER USIAUX, Mtdtiscoh two-way embedding schemes for free- surface primitive equations in the midtidis iplinsry simulation, estimation and assimilation systems, Ocean Dynamics, 80 (2010), pp. 1497 1587.

[21] PATRICK J HAL Y JR. AND PIERRE F. J LERMESIAEX, Multiscale two-way embetldmg schemes for free-surfaee primitive equations in ike muliidiscipimary simulation, estimation nd assimilation system. Ocean Dynamics, 60 {2010), pp. 1497 1537.

[22] BENJAMIN HESSELS, Time-optimal path planning o sea-surface vehicles under the. effects of strong cwrmts and winds, Bachelor’s Thesis, Department of Mechanical Engineering, Massachusetts Institute of Technology, June 2014.

[23] GHJU-YBN KAO, .STANLEY OSHER. AND YEN-HSI TSAL Fast sweeping methods for static hamilton-jacohi squahons, SIAM journal on numerical analysi , 42 ( 005 J, pp. 2612 2832

[24] NAOMI Emiicii LEONARD. DEREK PALKY, FRANCOIS LBKIEN. EODOLRHE SEPULCHRE. DAV ID

PSATAKTONi, AND Russ DAVIS. Collective motion, sensor networks, and ocean sampling, Proceedings of the lE-FE. special issue on the emerging technology of neiAvorketi control systems, 95 { 2007), pp. 48 74.

[25! P. F. J . LEKMtSSAUX, Error Suhspece Data Assimilation Methods far Ocosn FteM Esiimoiion:

Theory. Validation and Applications, Ph thesis. Harvard University, May 1997.

[26] P'. F . J. LER LBIAUX, Uncet &mty estimation and prediction far interdisciplinary ocean dynamics, Journal of Computational Physics, 217 (2006), pp. 178-199.

[27] P. F. J. LBBMVSTAHX, C-8. Cmu. 0; G . GA A KiEwtcz, P ABBOT, A . E. ROBINSON , R. .

MILLER, P. J HALEY, W. G. LESLIE, S. J MAJU DAR, A PANG, AND F LEKIEN, (Quantifying uncertainties in ocean predictions, Oceanography, 19 (2006), pp. 90-103.

[28] F. F. LBRMUSIAUX, P. i . HALEY, W. G LESLIE, A . AGAUWAL, O. LOGUTOV, AND L.J.

BURTON, Aiu!tiscale physical and biological dynamics in the Philip ines archipelago: Predictions and processes,, Oceanography, 2-4 (2011), pp 7(1-89.

[29] . LOIXA . Path planning in time dependent flows using lend set methods, master’s thesis,

Department of Mechanical Engineering, Massachusetts Institute of Technology, September 2012

[30] T. HOLLA , F. J . HALLY JR , AND P. F. J LERMUSIAUX, Time-optimal path planning in dynamic flows using level set equations: realistic applications, Ocean Dynamics, 64 (2014), pp. 1899- 1417. ;G3ί] T LOLLA. R . F. J LBKMUSIAUX, M. P. IIECKERMANN , AND P. J. HALEY J R. , Tima- optimal path planning in dynamic flows using level set equations: theory ¹ and schemes.. Ocean Dynamics, 6-4 ( 2014), pp. 137 -4 97.

[321 T. LOLLA. M. P, U ECKERMANN , K. YIGIT, P , J . HALEY Jg,, AND P. F. J , LERMUSSALX, Path planning in time dependent flow fieM Is using level set methods, in Proceedings of IEEE International Coisferenoe on Robotics and Automation, 2012, pp. 166 173,

[335 JOHN LYGEROS , On reachability and minimum cost optimal control. Automatic», 40 ( 20i ), pp. 917 - 927.

[34) IAN MITCHELL, Application of Level Set Methods to Control an lieaekabilit-y Problems In Continuous and Hybrid Systems. PhD thesis. Scientific Computing and Computational Mathematics Program, Stanford University. August 2002.

|35) IAN MITCHELL AND C’LAIRE TOMLIN, Oeerappr xsma ng reachable sets by hamMt&n-jacebi projections, Journal of Scientific Computing, 19 (2903). pp. 323-346,

[38] I. MITCHELL AND C . 3. TOMMN. Level set meth ds for computation hyl id systems, in Hybrid Systems: Computation and Control, .Nancy Lync and BruceH, Krogh. «ids., vol. 171 1 of Lecture Notes in Computer Science, Springer Berlin Heidelberg. 2000, pp. 310-323,

[37] L M MncHELL, A . M flAYEN. AND C, J . TOMLIN. .4 t -m-dcye enl Hamilton- Jacobi formu lation of rnnchahle sets fox- continuous dynamic games, IEEE Tramad ions on Automatic Control, 50 (2005), pp. 947 - 957.

[38 [ STANLEY OSHER. AND JAMES A SETHTAM, fronts propagating no h curvature-dependent speed:

Algorithms based on Hamilton- Jacobi formulations, Journal of Computational Physics, 79 (1988), pp. 12 49.

[39; _: A PitLPOTT AND G . LEYTAND, Rawing to barbados, in Algorithms for Optimization with Incomplete Inforenation. Susaime Albers. Rolf H. Mbkring, Georg C¾. Pfiug, and Rudiger Seim It», eds., no. 05031 in Dngst nhf Seminas· Proceedings, Dags iihl, Germany 2005. Internationales Regegimngs- und Farschungsssentram fiir Informat ik (IBFI), Seliioss DagstuM, Germany.

[441] Y OAV RAK, Sailboat speed os. udnd speed, Catalyst, Journal of the Ara&tenr Yacht Research Society, {April, 2000), pp. 263-274.

[411 B. RHOADS, 1. MEZiC, AND A POJE, Mininmm time feedback control of autonomous underwater vehicles, in Decision and Control (CDC), 201U 49* ^L IEEE Conference on, Dec, 2010, pp. 5828 -5834.

[421; BLAND RHOADS, IGOR MEKKT AND ANDREW G, POIE. Minimum time heading control of underpowered nekictm in time-varying ocean currents, Ocean Engineering, 66 (20131. pp. 12 - 31.

[43] J SETiRAN AND A . V LA lMiViBK ¥ . Ordered upwmd ethods for static Hamilton- Jacobi equa tions: Theory a d algorithms, SIAM Journal on Numerical Analysis, 41 (2003], pp. 325 - 363.

[44] J . A . SETHIAN , Fast marching methods, SIAM Rev., 41 ( 1999), pp. 199-235

145; J. A . SETHIAN AND A , VLADIMIRSKY , Ordered upwmd methods for static hamilton fjacobi cgua- Hans, Proceedings of the National Academy of Sciences, 98 ( 2001 ), pp. 11069 -11074.

[48] D . N , SUB RAMAN! A D P, P. J. LBRMUSIAUX, Energy-based path planning by stochastic dynamically orthogonal hv d-sei ap mi tion, Ocean Modeling. (2015) . in preparation.

[47] D. N, SDBRAMANI. T. LOLLA . P, J. HALEY in \L , P. F. J, LERMUSJALX, A stochastic optimisation method for energy-based path planning, in Dyaxnie Data-driven Environmental Systems Science Conference, Lecture Notes m Co ute Science (cds. Ravela, S. and Sandu, A.}, Springer international Publishing, 2015. in press.

=48:; DANIELA TONON, Regularity Results For Hamilton- Jacobi Equations, PhD thesis, Scuola Inter B&rionaie Super iore di Stu<li Avansati, September 2011.

[49] M. p. ITSCKERMANN AND P. F. J. LER DSIAHX. 2.39 Finite Volume MAT LAB Framework

Documentation, teoh. report, MSEAS Group. Massachusetts Institute of Technology, Cambridge, MΆ USA, ^* 2011 ^‘ ISO] ALEXANDER V LAD1 IHSKY, iftalic. pdes far ma-dependsnt canirol problems . interfaces an Fret; I osmdaries. 8 (2006), p. 281.

[Slj K. Y SGIT, Path planning methods for autonomous u derwater vehicles, master’s thesis, Department:· of Mechanical Engineering, Massachusetts Insiitsifee oi Technology, June 2011.

III. Stochastic Time-Optimal Path Planning In Uncertain, Strong, Dynamic

Fields

L Introduction

Planting optimal paths of aut nomous platforms in dynamic environments such as the ocean and atmosphere is hnpmtaat for maximall utilizing the platforms’ capabilities. in the ocean, commonly used autonomous vehicles - underwater gliders, propelled underwater vehicles and surface crafts - often undertake complex missions such as oceanographic data collect too. search and rescue operations, oil and gas discovery, and acoustic surveillance and security tasks | l-3|. Path planning is the task of predicting paths for these vehicles to navigate between any two points white optimizing some or all operational parameters such as time, energy, data collected. and safety. A related concept is dynamic reachability forecasts. the task of predicting die dynamic set of all the locations that can he reached by these vehicles.

One major challenge in dynamic reachability forecasting and optimal path planning for realistic ocean conditions is that current forecasts are uncertain. There cotikl be uncertainties in the initial conditions. hotmtlaiy conditions, parameters and even terms in dte equations themseives {4J. In fee present paper, owr objective is to develop fundamental ibid efficient stochastic equations and methodology for computing the teachability fronts and time- optimal paths of vefeicies navigating m uncertain, strong, and dynamic flow gelds. A« we will «cc. such an approach also helps us in quantifying the sensitivity at optimal paths fo errors in How predictions.

A traditional foens of path pSaiiiiing has been on autonomous agents in cort^vicv hot static environments. These traditional methods can be broadly grouped into graph-based search methods, nonlinear optimization meth s me hiding those with evolutionary algorithms, and dynamics-based approaches including Lagrangian Coherent Structures. Fvr a review ee. c.g . j> -- an references thetein. Recently. fnixkimeniai ievel-«ct Partial Differentia) quations ⁽ PDJEsl that govern reachability have been developed for lime -optimal (?.9.l0j and energy-optimal f i t ] path planning of autonomous swarms in strong and dynamic deterministic flows. This fevei-set. PDE-based methodology overcomes tbs limitations of traditional path planning methods designed for robotic motion in static environments. It can indeed directly utilize a prediction of the dynamic environment to plan time- or energy-optimal paths that inSeihgeody utilize favorable flows and avoid adverse currents, it also directly avoids physical obstacles ie g., islands) and forbidden regions due to operational constraints <e.g.. minimum water depth) thereby ensuring vehicle safety. It feus been employed to solve pursuit-evasion games in dynamic Hows 1 12,13). its use with realistic deterministic ocean reanalyses was demonstrated several ocean regions | I4J 5{ and in real-time sea exercises wrth real AUVs [lb]. However, the corresponding PDEs and methodology need to be extended to dynatnic stochastic environments.

Progress has been mad? on path planning for autonomous robots in uncertain environments Wfellman et al. f f?1 extended graph-based search methods to account for uncertain edge costs with applications to stochastic bus netw rks Kewiani et al. j iS] extend gr.iph-hns.ed algorithms tc explicitly consider uncertainty in the mobditv and terrain parametenzations using stochastic surface response methods. Potential field methods have been used in a Monte Carlo approach to solve path planning problems in uncertain Sows by Batxaquand and Laiomhe [Iv], Rathbun et al. :[2(!| report. the use of evoiuiiorwry algorithms to plan paths of Unmanned Aerial Vehicles IUAV) m an airspace with uncertain obstacles, Wang et al. 1 f . I perform path planning of autonomous vehicles to dy namic and uncertain ocean currents nsing an ensemble apprt h The authors solve the deterministic Boundary Value Problem fBVP) for each member of the ensemble to compote the statistics of optimal trajectories. Considering all t these approaches, on the ne hand, the extension of graph- based methods is not suitable for dynamic environments as iheii cost Increases exponentially and there is & lack of optimality guarantee. Gn the other hand, the Monte Carlo and ensemble approaches can have slow convergence and prohibitive costs that increase with the complexity of uncertainty. For optimal path planning in dynamic stochastic environments, fundamental and efficient stochastic PDEs {S PDEs> are needed

In what follows, we first outline die problem statement, present oar exact planning level-set S-PDEs and new efficient stochastic DO level-set PDEs (Section 2). Next, we showcase applications that validate and demonstrate the capabilities of our equations (Section 3). Finally, we provide some concluding remarks (Section 4).

2, Stochastic time-optimal planning PDEs

2. /. Prehistn statement

Consider an open set 'D e 'R ^{x å} ' {hc· = 2 or 3 for 243 or 343 in space), lime t e [0, co} and a probability space B P.i, where W Is the sample space, E is the o-algebra associated with ff, and V is a probability measure. Starting atr = 0, let a vehicle {P) navigate in 2> from x. _t to x ,· with a nominal speed F(#) > 0 under the influence of a stochastic dynamic flow-field vfx, , <a) : V x [0, oo), where w e W is a random event. Let j>(x _s , f; w) be a general continuous trajectoiy from x^to x A(¾, 3- i _, l.

The main notation is provided in Table 1. Our goal is to predict the stochastic reachability fronts and time-optimal paths.

2,2, Stochastic level -set partial iffemmsi equations

We approach the above stochastic time -optimal path planning problem starring with the exact level-set PDEs for deterministic path planning [?!· Hero, the source of stochasticity is the uncertain dynamic How v(x, t; to) (in comparison. the vehicle speed was artificially made stochastic lor energy optimization in |ll }>. The stochastic reachability from for P is thus governed by the stochastic Hamiilon-Jacobi (HJ) level- set equation with the Initial condition f(c, 0: a>) = jx - x _sj : and, where needed, open boundary conditions such as j = 0. whom n is the outward normal to the boundary d2>. m is a random event, and F is the re3ebabihty-fnxit-iracktng scalar level set field (e.g.. signed distant'.' Junction. For every w. the opti al arrival-time T {% ·. mt at x _f is obtained by integrating Eq. t : ; until the first time i such that fί /. r; &> < 0 The corresponding optimal srajtvt t X pt s,. f; e» is then given by the particle backtracking equation (where f is differentiable}

Soh ing the S-PDE\ (D and ;.N would fulfill ur good. Eq. 1 can be hiiegrate by direct Mottle Cark· methods; however, the computational cost of such a metho is prohibitive for practical applications. Classic reduced-order stochastic methods could be utilized [2 -25J. However due to the variable propulsion and advectien terms, and the highly dynamic flows, they can sometimes diverge or be melfrcient. We thus employ a dynamic stochastic order reduction. specifically the Dynamically Orthogonal iDOi field equations |2b.2?] that weie shown to be m instantaneously -optimal reduction using differential geometry arguments s j. Nest, we obtain these efficient stochastic DO fcvel-set equations.

2.3. Stochastic Dynamically Orthogonal level-set equations for uncertain flow fields

A brief overview of the DO theory is provided in Appendix A. As mentioned previously, in [11] we obtained the

stochastic DO level-set equations when the vehicle-speed Fit) in Eq. ! ) was made stochastic (i.e.. Fit) ® Fir. w)).

Here. She difference is that we have a deterministk F(t ) and fee stochastieity comes from the uncertain fio -ikid v{\, < . _<> · _> . This stochastic ilow-Hdd can be computed before solving the stochastic levei-set: equations and considered cx .m ίl!a·!>hί ·h :hV to the 5-PDE · i

Lei us introduce DO decompositions to the uncertain 8ow~fteid v and the stochastic fevel-sei f as

where m<{h wϊ, V j = 3 . . , « _{s v ,} and Y t; <o), Yi = 1 . . , are zero-mean stochastic processes that represent the generatiy complex probability density functions of the veioeity and level-set fields respectively.

Henceforth, we drop the spaiiai and temporal variables in parenthesis for brevity of notation. We use repeated indkes to indicate a summation over », _, * velocity DO modes oi- n _s a level- set .DO modes, as the case may be. ¾¾ also define y º \Vip

Substituting Eqs. (3) and (4) in Hq. ( I }, we obtain

We then derive the mean equation by applying tbs expectation operator to Eq. (¾. the coefficient equations by projecting Eq. (5) onto the modes </>;, and the mode equations by multiplying Eq. (5) with i ^' ,: and then applying the expectation operator. We provide the detailed derivation in Wei [29 J. The results are as Ibtlows:

isi'r C denotes the covariance matrix iietween the two stochastic variables indicated in the respective subscripts (Tabs* i }.

Here we do not consider a separate DO decomposition for the non-polynomial nonlinearity g bat choose to hiffidfe only that te:rmreaifeation-by-realizatk>R and compute the required statistical quantit'es ¾ y] and å{« y j. Other approaches tor handling the non-polynomial nonlinearity, viz. , KL-Gamma and Tay hr- Gamma j i I] could also be utilized (not shown here).

The solution of DO E s. f§) provides reachability fr nts for ail realizations >¾. The minimum arrival time for each realization is tint first time t for which <fex _f . t: w) < 0. Subsequently. the stochastic particle backtracking Eq. 2 can be solved to obtain foe time-optimal trajectories X (x„ t w . The solution also provides the tone-series of time-optima i headings h* tr; to), again for all w. in Section 2, we validate (he solution obtained by the DO level-set E , <ht by comparing it with that obtained from direct Monte Carlo (MCs methods. We also showcase applications for pmdicting the stochastic reachability frcmts. ads. and time-optimal trajectories for vehicles nav igating in uncertain flow- fields.

2.4, Cmtp ma} easts

The computational cost of solving the level-set PDE ( I i for w fixed is commonly dominated by the cost of the adveetion term. Thus, the cost of a MC method scales with the number of realizations and grid points as 0<n _r n _? ). The cost of solving the stochastic DO level-set equations is also often dominated by the adveetion terms (b the velocity mean and modes). This cost scales as £>(jX, _v + T i)¾) and is independent of n _r . .«,· only affects the cost of solving the coefficient Eq. (fib), which is simply an ODE independent of»,,. For typical values of»*,, », , B, _¥ . and «,· _, ,· of realistic applications, tlx: computational speed-up can be several orders of magnitude.

3. Applications

We apply the efficient stochastic DO level-set PDEs to compute stochastic teachability fronts and time-optima] paths in three tost cares. They correspond to different stochastic Ro configurations: (i) a simple canonical s ¾¾¾teady front w here only the strength of file flow is uncertain, (ii ) a stochastic dynamic barotropic quasi-geostrophie double- gyre circulation;, anti (iii) a stochastic Row past a circitiar island. We employ the first test ease to verify the solution of the DO equations by comparing it: with its cotrespondhig Monte Carlo (MC> solution, in all test cases, we describe and study the variability in reachability fronts and time-optimal paths. The numerical parameters are listed in Table 2.

3.1. Test Case l; Stochastk steady-from with uncartainftaw .« tgift

The first test ease is a canonical Row— a steady zona! jet with uncertain flow strength, in addition to being a good candidate for developing and testing numerical schemes, this flo scenario is commonly encountered in the ocean. For example, the crossing of a sheUbreak front [is] or of channels [16] caa he idealized by this canonical flow. Hence,

Taiffc 2

.MkΐΉΰt 'Lΐί p5f i sr siii ί¾t·· fcesi iistf « isiuw dciffisr ia T&Mί> Ί }

Tfeftt Case !!, fri /> i is. ft :.v

m 0 t t 5L--2 2000 i ? too too o.h! aot ts-5 sm s

3 240 00 at»? 0.067 !e-2 «3009 S Studying the properties of time optimal paths in such flows is a valuable firs step for planning in realistic ocean flows later on.

The domain, flow configuration (mean, modes, and coefficients), and start and end points are provided in Fig.

3.2.

Fig. 3.2. Mean, Mode, and PDF of coefficient for ilie DO decomposition of she velocity field used in Test Case 1: The zonal jet is t« m West to East between / = 40 ami >' = 60, and has tin uncertain strength of uniform densiiy tlisirihihion with lower limit 0.5 and upper limit f .5. The velocity is 0 cist «'here it; the domain. The start point is mariied with a cireuiar marfcer and these end points with a star marker.

The front is a steady jet flow f ram w est t o e ast, t hat t s s onbitcd an d c nstam tn th e re gion 40 < y 60. The strength of the flow is a random variable uniformly distributed from 0.5 to 1.5 with the jet proper. The flow is zero elsewhere in Ihe domain. To represent tire stochastic steady trout w ith a uuitor ly distributed uncertain strength, we need only one DO mode for velocity, i.e., n,. _v = f . We utihzc a. ~ 2 Ottfl realizations to re resent the uniform PDF of flow strength ⁽ see Ft». 3 -2a}. We consider missions with a vehicle moving at a non-d imensional nominal speed of Fit) — 1 , V i— (0, oe). The start point is (150.20) and the three Targets are 1 0 8tn. 1 1 0,80), and (210,80). We compute the reachability fronts and sets by solving ihe stochastic DO level-set Fqs "·' and time-optimal paths by solving the. backtracking Eq, (2 i

First, we verify the solution of Eqs. ⁽ 6) by comparing it to the solution of Eq, ⁽ 1 ) computed by the MC method. Then, we study the properties of the reachability fronts, sets and time-optimal paths.

3. LI. Verification of ike DO solution

To verify the DO solution, we first compute t he corresponding M€ solution and then compare the two. The MC solution is obtained by solving the level-set PDE for each of the n, realizations of the velocity field separately. We look at the differences in arrival time at the target (210,80) and in the reachability front of all realizations.

Pig. 3- i shows the histogram of relati ve error in arrival tune at ihe fatget (210.80).

The arrival times computed by the DO method are nearly identical to the corresponding: times computed by the MC method. 26.55% of ihe realizations have zero error and. 82.8% of realizations have less than 0.1% error in arrival times. The maximum error is 024%, which is 0.15 non-dimensional time or three discrete time-steps. Overall the etror in arrival-time is negligible, especially considering that the DO method is 4 orders of ma ni ude faster than the MC tneihod.

To quantify the difference in the reachability front, we utilize the discrete Ftechei distance [1 1,50,.'· ! ] which measures tire maximum distance between two discrete curves in 2-D. Fig:, 34 shows (he discrete- Frechet distance normalized by the grid spacing for all realizations, plotted at four non-dimensional times. All Frechet distances are less than the grid spacing, which indicates that the reachability front computed by the DO and MC methods are identical and the errors are less than or of the order of the spatial discretization. We also compared qualitatively by plotting the reachability fronts computed by both methods for some characteristic realizations and observed that DO and MC compute identical reachability fronts (not shown here, see }29]). We thus verified that the reachability front and arrival times computed using the DO level-set equations are· accurate for our applications. All Freehet distances are less than the grid spacing, winch indicates that She reachability front computed by tire DO and MC methods ate identical and the errors ate less than or of the order of the spatial discretization. We also compared qualitatively by plotting, the reachability fronts complie b both methods for some characteristic realizations and observed that DO and MC compute identical reachability fronts (sot shown here, see |29]}. We thus verified that the teachability front and anival times computed «sing the DO level-set equations are accurate for our applications.

3.1.2. Analysis of stochastic reachability fnmls ami time-optimal paths

Fig. 3-5 shows the distribution of She stochastic reachability· fronts computed by solving the stochastic DO level- set Eqs.

Each reachability front is colored with the flow· strength of that realization. The six panels show snapshots of the spatial distribution of the reach ability fronts at six non-di ettsionttl times.

There is no uncertainty in the flow field outside the jet proper. Hence, (he reachability fronts for all realizations are identical until the time when they first experience the uncertain flow. As shown in Fig. 3-5b. the part of the reachability fronts in the jet proper spreads out while the part that has not yet reached the jet remains identical for all realizations. As time progresses, the spread in the reachability' fronts increases as each flow realization has a different flow strength.

(Fig 3-5c-f). Most notably, there is a locus of points where pairs of reachability fronts corresponding to different flow realizations crossover in Fig. 3- 5c and d. For target points downstream of this locus of crossover points, the higher the flow strength, the faster the vehicle arrives at the target. On the other hand, for target points upstream of this locus, fire lower the flow strength of the stochastic· jet, tire faster the arrival-time at the target. In tire former situation, the flow aids the motion of the vehicle to the target and thus stronger flows are favorable; however, in the iattei _' dluutiof . li e flow hinders tire motion of the vehicle to tire target and thus weaker flows are favorable. By t = 75 , ; * V \ u'h ies .» all flow realizations have arrived at the target (150,80 !. By i— 100 (Fig. 3-5e), vehicles in all flow re ill utn ns have crossed the targe! (210,80 > also and the reachability front of the weakest flow strength is just arnvmg at tire target (90.801. By t— 125 f Fis. 5 -5f>, vehicles of almost all flow' realizations have crossed (90,80).

Fig. 3- shows the stochastic time -optimal paths for ail three targets. As before, each path is colored with the strength of the flow realization corresponding to that path.

F¾.3-6- Stochastic timeoptimal patitsprTm Case 1: Ail tin -optim&l paths ate colored with the flow strength of the corresponding flow realization. The variability of the titne-o _j ximal paths is greatest for the target point upstream of :he start point

The spread of tire paths increases progressively as the target moves upstream front the start point. The spatial distribution of the time-optimal paths to target (210,80) is narrow, in other words, the sen itivity of the time- optima! paths to error in the flow field prediction is low. On the other hand, Sie spatial distribution of the time- optimal paths to (90,80) is wide. The sensitivity of the paths to errors in the flow field prediction is higher than the other two targets. For this test case, the more upstream the end point is, the larger the effect of the flow strength uncertainty is on the time-optimal paths, non-dimensional stochastic baroiiropk quasi-geosteiphie model (e.g.. [32,331). written as conservation of momentum in the Lange v in form,

where Re is the Reynold’s number, / is the Coriolis coeilcteist. and a is the strength of the wind stress. For die Coriolis coefficient. we employ a f -pki approximation f = f _r . ÷ bg with / _e = 0. We utilize a horizontal length scale L— 10 ^s m. vertical tengili scale D— UP m, velocity settle U— 1.9S ah/s, lime scale T— 1.6 yrs. density p = 1025 kg/m' stress scale r« = 0.16 MPa, eddy viscosity si _¾ ,· = 19.77 nr/s, and p _Q = L977¾ 10 ^{~ : f} ms. This leads to non-dimensional numbers Re = y÷ — 1000. a— , ^i: . = 1 , 000, a i B— ¾E ^~ /U— i , 000. A deterministic non-dimensional steady zonal wind stress, t— . forces the flow in the basin. The DO equations corresponding to the stochastic barotropic quasi-geostrophic dynamics Eq, < 7; are provided in A pen ix & .2.

We sirnuiate tlie flow in a basin of size 1000 km x 1000 km by solving Eq. (A.3) and using the relevant numerical parameters in Table 2. The barotropic zonal and meridional velocities are initialized (Eq. (7e s) from a spatial correlation kernel with a length scale of 5 ⁽ X ⁾ k and uniform variance. Five velocity DO modes and 5000 DO realizations are utilized (Table 2), and the simulation is span-up for 1 non-dimensionai time, i.e.. 1.6 yrs. For path planning, we utilize the next 13.5 days of stochastic flow. We consider a vehicle with a nominal speed of 40 cm/s (|34|, e.g. ). Fig. 3-7 illustrates the DO mean, DO modes, and marginal PDF of the DO coefficients for Ore flow field at the beginning (i = 0 days) and end if = 13.5 days) of the planning horizon.

Over the 13,5 day period, the tnerin flow is a strong zonal jet from east to west with a cyclonic gyre to the north of the jet and an anti-cyclonic gyre to the south (Fig 3-7 Mean). From the variance of the coefficients (Fig. 3-7 Variance of Coeffs.), we see that Mode 1 has tire most stochastic energy and explains the broad features of tire uncertainty in the flow. Mode 1 has two zonal jets and three gyres (Fig. 3-7 Mode l ). The two jets are located north and south of the zonal jet of the mean. The stochastic coefficient corresponding to this mode has a bimodal distribution. For positive coefficients, tire southern jet (0,2 < y < 0.5} contributes a flow from west to east and die northern jet (0.5 < y < 0,8) from east to west. For negative coefficients, the direction of this contribution is reversed. Tire northern and southern jets can Ire locally considered to be similar to the stochastic front crossing example (Test Case ! }. Higher positive (negative) coefficients correspond to greater variation from the mean m the (opposite) direction of the modes. Modes 2—5 and their coefficients contribute to noil-symmetric flow realizations (Fig. 3-7 Modes 2-5), such as tire non-sym errie wind-driven gyres seen in the North Atlantic (e.g.. [35.36] ). Two realizations with most negative Coefl. 1 (Realiz. # 1 ) and with most positive Coeff. I (Realiz #5000) are shown in Fig 3-3.

Pig. 3-8, Tm> rwUmbms of she swchas c flaw field for Test Case 2; Sireaiilioes of sire flow lire overlaid oa a color plot <·( the flow jusgitUude. Reaiiz. #1 corresponds io the roost negative Coeff. 1 sad realiz #5000 to tits roost positive Coeff. i Alt such realizations ate integrated in tiine bv one DO flow field siffiuiaskni.

Both these and 4998 other realizations are all simulated y one DO flow field simulation (Appendix A).

3.2.1. Analysis of stochastic reachability fronts and time-optimal paths

Now we plan stochastic time -optimal paths from a start point (0.2, 0.2) to the target (0.4, 0.8). Fig. 3-9 illustrates the evolution of stochastic reachability fronts (i.e., zero level-set contours) over the planning horizon, showing snapshots at nine discrete times. Fig, 3-9. Stochastic reachability fronts for Test Case 2: The reachability from for each of the flow resiizsfions is colored with its corresponding velocity DO coefficient 1, m m) Time ; is in toys, iiiid.v - an y - av.a, are nr thousands of fcm.

Each flow realization has a reachability front and all such fronts are computed by one DO simulation. Overall, the reachability front gets advected by the mean flow, but individual realizations experience the effect of the modes and coefficients of the flow. As such, the reachability fronts spread, leading to a density distribution (Fig. 3-9). The mean double-gyre and clockwise in the cyclonic northern gyre. To study the spread, we look at Fig. 3-9. Here, the reachability fronts are overlade on streamlines of velocity DO mode

1, v~j x, r), and colored with velocity DO coefficient 1 , m _¾ (/; w). The initial spreading of the reachability fronts (Fig. 3-9b-d) in the southern part of the domain (>· < 0.5) is dne to the flow contribution by the southern jet of velocity DO mode 1 . The reachability fronts corresponding to positive coefficients are advected to the east more than those corresponding to negative coefficients. The spreading of the reachability fronts in the northern part of the domain (y > 0.5 i (Fig. Ffli-i} is due to the flow contribution by the northern jet of velocity DO mode f . Here, the effect on the reachability fronts is opposite to the effect due to the southern jet. i.e., more positive coefficients are advected to the west more than the realizations corresponding to more negative coefficients. The spread of the reachability fronts in the strong mean zonal jet {approximately 0.45 y < 0.55 ⁾ is not affected much, as the contribution to the flow by all velocity modes in this region is weak (Fig 3-9c-f). The finer details in the shape of the reachability fronts are due to contributions of the velocity modes with lower stochastic energy (modes 2-5).

Next, we study the spatial distribution of the time-optimal paths (F;g >· sis)

Fig. 3- id. Siockaslie !iim-op!itnai paths for Test Case 2: The lime-optima! p ths are colored with (a? the velocity ΪCΊ coefficient i. /se e op and id) the artivai time (in days) id the hit get shown ;t- and y- asms aie in thousands of tan

In Bg o-iOa, each path is colored with the corresponding flow realization's velocity DO coefficient 1 , mrί; aή, la Fig. iOb, each path i colored with the arrival-lime for that path. Ail paths cross the strong mean zonal jet by riding the anti-cyclonic southern gyre first and then the cyclonic northern gyre. The spread in the paths is a result of the variability in die reachability fronts as described above. The most negative coefficient /; _¾() ; &;) defines the western edge of the spatial distribution of paths and the most positive mjit; to) defines the eastern edge as seen in fig. .V K-a. When paths are colored by the arrival time (Fig. Mtib), one sees that the arrival time for paths that are very' different spatially (and corresponding to different flow realizations) can be nonetheless similar. This is due to the complex PDF of the stochastic velocity' field, in this test case, the arrival time range varies from 12.1 to 13.6 days. For reference, under no flow conditions, a vehicle with a speed of 40 cm/s would take 18.3 days to travel between the start and target points considered in this test case 2.

The examples so far also show that the time-optimal path can be very sensitive to the uncertainty in the flow field A strength of ihe stochastic DO level-set equations is that they allow the rigorous prediction of this sensitivity in a computationally efficient manner. The examples also raise the need of computing risk-optimal paths in uncertain flows, a topic addressed in [3? j.

3.3 Test Case 3: Stochastic flow past a circular island

In the third test case, we consider time-optimal path planning in a stochastic flow behind a circular island in a channel. Fig. 3-1 ! shows ihe domain and flow configuration (vonkity overlaid with streamlines). The flow is governed by the stochastic barotropic quasi-geostrophic Eq. f?c), but with/= 0 and t— [0 , 0] ^J .

Fig. 3-11. Flow field for Ten Case 3: R«v A stows the DO «teas slow, with streantllass overlaid on colored voriietiy. Row B shows the same for riii ifc ! ami sank· 3 feiiis. Row C stows she marginal PDT of eiefSAcieuls i and 2. ow D shows the deca in v riance, of the first eight modes Each column Is a snapshot, at the time specified at die top. The starting point of the propelled vehicle is indicated by a ratlar marker and the 6 targets by star markers We model a channel 16 km x 6 km with a eficuhtr island of diameter l km. Here, we utilize a turbulent eddty viscosity A _jf = 10 m¾, length scale L = 1 ton, velocity scale V = I m/s, and time scale T = 1000 s, leading to a non-dimensional Reynold’ s number of 100. We consider a deterministic barotropk inlet into site channel of v = [1 , 0]·' m/s on ihe western boundary. The northern and southern boundaries have free-siip wall conditions and the eastern boundary has open conditions. The barotropk initial conditions arc generated from a covariance kernel with a spatial correlation length scale of 5 kni in Ihe zonal direction and 2 km in the meridional direction. 8 velocity DO modes and 10,000 DO realizations are employed here {Table 2). The siniuia.tioa is sptm-up until stochastic eddies are generated downstream from the island Then, the flow for Ϊ .5 b (non-dimensional time of 5.4} is considered for path planning.

Fig. .>· ! i shows ti¾e velocity mean i row A>. the first iwo velocity DO modes crow B ), marginal PDF of the corresponding DO coefficients i row C >. and the variance of ilrst eight DO modes trow D o Each column of Ftg. 3 i : s ow s I he above quantities at t — 0. 2700, and 5400 s. From die dee.n in c aiianec i r w D). w e see th.tt most of the flow uncertainty is captured by die DO modes 1 and 2, and hence other modes and coefficients are not shown i but we used 8 modes i. The mean flow accelerates to the «outii an·! north of the inland as « i s confined b\ the channel and has almost zero magnitude in the lee just behind the island. The stochastic initial conditions result m uncertainty in the eddy strength, shedding frequency, and whether the eddies are first shed to the north or south of I he island. The modes I and 2 have eddies downstream from the island which together with the coefficients explain (he uncertainly in Site How downstream

Fig, 2 shows two realizations corresponding lo the most negative (Realiz. # ! ) and positive (Realiz. #10,000} coeff. I

We see that in Realiz. #1 , eddies shed to the north of the island while for realiz, #10,000, they shed to the south of (he island.

3.3. /. Analysis of stochastic reachability fronts and tttne-oprimal paths

Now we plan stochastic time-optimal paths for navigating a propelled vehicle with a nominal speed of 1 m/s from a start point (2.3) upstream of the island to six possible target points (see Figs. 3-1 1 -1 : the star! point is indicated by circular markers and targets by star markers). Fig. 3 : 3 shows the evolution of the stochastic reachability froirt overlaid on streamlines of velocity modes i (column A) and 2 (column B).

In each plot, we color the reachability fronts with fheir respective coefficients. The level-sets do no! grow west of die start point as the vehicle speed and the mean flow there are equal and opposite, making the reachability front stationary at the start point. Until t— 1800 s, the reachability fronts for all flow realizations are almost identical as the flow uncertainty in the reachable sets is insignificant i Pis;. 5-OA.i, B.i). The shape of the reachability fronts is set by the mean flow between t— 1800 and 2400 s. The reachability from: grows normal to itself at the nominal vehicle speed in the dose-to-zero mean How region behind the cylinder (Fig, M3A.fi, B.ii). After t — 1800 s, the spreading of the reachability fronts due to the effect of velocity DO modes start. The mean flow advects all the reachability fronts to the east. From t = 3000 s, the reachability fronts start spliting inio two groups roughly corresponding to positive and negative values of coefficient i . At t = 3600 s, the group of reachability fronts with positive coefficient 1 have a kink near g - 3.5 and those with negative coefficient 1 have a kink near y— 2 5 corresponding to flic eddy direction i ;¾. 3 - ;3A.ϊn). At t — 4500 s, the reachability tronis have two groups corresponding to positiv e and negative coefficient 2 (Fig. 3 3A.v). At t 5400 s, the bchav tor E similar to t— 3600 s, with two distinct groups of reachability fronts corresponding to positive and negative v alues of coefficient 1. Next, we study the spatial distribution of the time-optimal paths (Fig. 3- 14).

Fig.3-14. Stochastic pmhs or Ttst Case 3: Cede*»» A ίb» fcas paths colored by velocity i irffitatrt i t ^’ 2). Row ; shows aii paths to target 2, row « to target 5, and row iii ·o bk other 4 targets. The v- and y- axes ate at kins.

In Bfs. 1· 14, column A has all paths colored with velocity coefficient 1. and column B has all paths colored with velocity coefficient . Row i has all paths to target 2, row if to target 5, and row in to all other targets. Paths to targets 1. 3. 4. and 6 (row Hi) have little variance, since these paths are mostly outside the regions of flow uncertainty'. Paths ro target 5 (row n) have the maximum variance followed by paths ro target 2 (row i), as these are· affected by the flow uncertainty. There are two seis of paths: one· set that goes north of the island and another that goes south. The combination of velocity modes and coefficients creates flow realizations that lead to northern paths or southern paths. This bifurcation mirrors the uncertainty in eddy shedding downstream from the island. An error in estimating the direction of eddy shedding will result in large errors in the time-optimal path predictions to targets 2 and 5 (e.g., north vs. south) for this test case. However, predictions to targets 1, 3, 4. and 6 are largely unaffected by such an error. This illustration quantitatively confirms the fact that sensitivity' of path predictions depends not only on the flow itself but also on the target locations, and all such paths and sensitivities can be one reachability front simulation is required for one start location and start time. The paths to multiple targets are obtained by solving the backtracking stochastic ODE Eq. (2), which is very inexpensive to solve.

fig. 3-15 shows the distribution of arrival times at the six targets.

As expected from fee distributions of time -optimal paths, the distributions of anrival times at targets 1 , 3, 4, and 6 are tighter compared to other two targets. Target 5 indeed has fee most variance In arrival time, and the distribution is bi odal similar to fee bitnodal distribution of velocity DO coefficients I and 2,

4. onclusion

We obtained and applied fundamental equations far ti sue-optimal path planning is uncertain, dynamic and '4RH¾> Sows. We first presented fee stochastic level-set PDEs that govern fee exact evolution of reachability fi -vds ^' < : vehicles navigating in uncertain, strong, and dynamic flow fields. ¾ then developed efficient stochastic dynamically orthogonal level-set PDEs that solve the above equations in a reduced dynamic suhspaee, providing several ortfeos of magnitude computational speed-up when compared to direct Monte Carlo methods. To compute stochastic time- optimal paths, we «di ed the governing stochastic particle backtracking ODE. We then applied fee equations to compute stochastic teachability fronts and time-optimal paths in three different scenarios: a steady-front f uncertain strength, a stochastic douhte-gyre circulation, and a stochastic Sow past a erretdar island. We utilized the first test cave to verify that t>ur DO level -set equations can compute fee distribution of stochastic reachability fronts and tune- optimal paths as accurately as traditional Monte Carlo methods. In the second case, we quantitatively explained she effect of fee first DO velocity inode (with the most stochastic energy I on setting fee shape and distribution of the reacfaahtiity lions and utne-opumal paths, We also showed feat fee paths corresponding fo different flow realizations can have large spatial differences hut similar arrival times, doe to fee complex PDF of the stochastic velocity fields. In the third case, we described how the variance of time-optimal paths and arrival times depended mostly on the first two DO velocity modes. We showed that the variability could be large due to the uncertainty in estimating the direction of eddy shedding for targets directly downstream of an island. Impertanriy, all these paths are computed by one simulation of fee stochastic DO level-set equations. Overall, our analysis offers insights into the behavior of vehicles navigating in canonical uncertain flows often encountered in coastal ocean regions. In addition to computing stochastic reachability fronts and time-optimal paths, the new equations provide a rigorous framework for quantifying the sensitivity of time-cptim&i paths to variability and errors in how field predictions. In the future, our methodology can he applied with realistic uncertainty flow predictions (e.g., [16,38,39]}, The stochastic DO level-set equations can fee augmented with decision theory to compute risk-optimal paths for vehicles navigating in uncertain flow fields. Even though we emphasized (he navigation of autonomous vehicles in ocean flows, the methodology and equations are general: they can be utilized far serial vehicles, land robots, or ships (e.g. , [401).

APPENDIX A

In this appendix, we first outline generic Dynamically Orthogonal field equations and then provide the barrotropic quasi-geostrophic DO equations. We refer to Table 1 for the notation.

A /. Generic DO equations

Prior to describing the Dynamically Orthogonal held equations 126,27,41], we fust provide some definitions. For any two fields *> and ¾ x t », the spatial inner product is defined as

dx , (A. i t

the orthogonal com onent operator as

an the ΰc bΰ&i oh o or&i r as

Consider the time evolution of a general stochastic continuous field y(c . ί; w) (e g. , ¾9(x, r; a»)

level-set scalar field and ^( ^' v . t ; co) = r(x, r; re) for the velocity vector field), described by a S-PDE

where C is a general nonlinear operator. Introducing a generalized dynamic Karhanen-Loeve (KL) decomposition (a DO deeomp·: t ition)

into (Eq. iA 4\s. appl ing the expectation ifcq. (AJ)) and spatial inner-product (Eq. (A 1)) operators, and an orthogonality condition on the evolution of the stochastic: subspace.

we obtain (he DO equations for the statistical mean y, stochastic coefficients T . and deterministic modes y; as

where we haw dropped (x, r) from the iraier prc liici femis for brevity of notation. Here, the n _s ^ motfes . aii stochastic coelftcients 1',· are dynamic and governed by different iai equations, distinguishing them from classic empirical orthogonal functions and proper orthogonal decomposition approaches. These dynamic nwtfes describe an evolving snbspace, and file coefiicients evolve the unc^tai ty within that dynamic sobspace. For nonlinear dynamical systems, the intrinsic rn Unearitles isEq. (A 4) ansi the correspcntding dynamics retained in Eqs. i _, A.7&HA,?e) often enable tmneating the number of » _« . modes an coeffickints to a mnch smaller iiitsiiher than the dimensions (spatial and stochastic) of the original dis retizarios of Eq (A 4). It is this dynannes and the a l sp X'o _¾ , fiiat allow die DO equations to accurately capture most of the stocbasticity of the original variables [27.42,43],

42 Dynamical f? orthogonal quasi- gewwophic equations

la this appendix, we provide the Dynamically Orthogonal barosropfc Quasi-Geostrophie equations that we utilised to generate stochastic Sow fields. Detailed description of similar equations for Boossinesq dynamics, their derivation and numerical schemes are provided in J41] and the implementation is provided in [44 j. The equations for the DO mean (Eqs. ¾ A,8&), iA S ), coefficients (Eq. (A cO, and modes (Eqs. (AM), <A.¾e>) of the S-PDE (T> are

where we have dropped the parenthesis tx . t; <») for brevity of notation. We solved Eqs. (A.S? to obtain the stochastic. Sow fields utili zed in Test Cases 2 and 3 of the present paper (Section 3).

References

IV. Risk- Optimal Path Planning In Stochastic Dynamic Environments

1. Introduction

The present work is based on fundamental differential equations that govern the evolution of the reachability and time-optima! paths in strong and dynamic currents (Lolla et ai. 2ul4b), andis related to energy- optimal paths (Sitbramani et a!., 2015: Subramani and Lermusiaux, 2016). These equations were used to compute optimal paths, both in realistic data-driven simulations (Loll t ah, 20i4a; Sobramani et a!,, 2017a) and with real vehicles (Subramani et a!., 2017b) They were also employed to solve pursuit evasion problems tSun et aL. 20 l ⁷ a.b ^' t. The input probabilistic predictions of flow fields are obtained from the variance-optimal reduced-order stochastic dynamically orthogonal equations (Sapsts and LemniMaux, 2009; Ueekermann et a!., 2013; Feppon and Lermusiaux, 2018 >. The resulting stochastic PDEs (S-PDEs) for time-optimal path planning have some advantages (Subramani. et aL, 2018); (i) for a given stochastic environmental flow prediction, the tochastic time-optimal paths are exact; (ii) the computed paths naturally avoid stationary and dynamic obstacles; and (Hi) the probability of a location being reachable {or non- reachable) are directly predicted.

From the probability distribution of time-optimal paths, we have to assess the risks and make a decision of risk-optimal path. This subject of decision making under uncertainty has been well studied in the fields of economics and management. One widely used model is the expected utility theory and its several variants as reviewed in Sehoemaker (1982 ⁾ . The key ingredients of this model are the evaluation of the -utility cost of the outcome due to a decision and the proba bility of that outcome. The expected utility theory can be utilized in a prescriptive or normative framework to inform optimum decision making under complex decision scenarios and can be customized to fit the risk behavior of users (Sehoemaker, 1982; Epstein, 1992; Von Neumann and Morgenstem, 2007). Specific utility functions are also available for different risk behaviors (e.g. Arrow, 1958; Fishbura. 1988; LiCaizi and Sorato, 2006). In literature, the combination of expected utility theory ⁷ and robot motion planning in dynamic environments has been limited. For example, Rapidly-exploring Random Trees (RRTs) with utility-based random trees have been used for single-query robot planners (Bry and Roy. 201 1). An expectation driven iterativ refinement approach with a heu-ristic has been proposed for robotic path planning problems (Boddy and Dean, 1989). The A ^* search approach has been combined with expected utility (Burns and Brock, 2007). A maximum utility based central arbiter for distributed architecture for mobile navigation has been developed (Rosenblatt, 2000). A risk-aware planner based on a Markov Decision Process has also been used to minimize the expected risk of surfacing of gliders along paths (Pereira et al. , 2013).

Our cu! S here i to combine a principled nsk siptimaiif criterion grounded in decision theory \\ ith tnrr stochastic dynamically orthogonal level-set equations to develop efficient compofaik al schemes lo predict risk-opthna! paths from the distribution of stochastic time-optimal paths. We also seek lo apply the ne« schemes lo several stodantic How scenarios, analyze the effects of ilitlerem risk metric and enteric. and iscu s tire p>opcnics of nsk optimal paths.

The remainder of this section is organized as follows. Next, we provnle a brief review of prior path planning results. In Sec. 2, we develop the theory and scltentes for risk-optimal path planning. in See. 3, we apply the new schemes to compute risk optimal paths for a variety of stochastic flow scenarios. In Sec 4, we conclude and provide futu research directions.

1.1. Previous Progress in Optimal Path Planning

Traditionally, path planning has focused on land based robots in stationary environments (e.g., Hwang and Ahuja, 1992: La Valle, 2006; Latombe, 2012). However, a major challenge for

marine and aerial platforms is the effects of the uncertain, strong and dynamic currents/winds. Several authors have extended many of the algorithms for static environments to plan paths of autonomous vehicles in dynamic environments ' for reviews, see e.g,, Lolla et al, (2014b); Pe reira et al. (2013)}, These applications include graph based search methods such as the modified Dijkstra’s algorithm (Mannarini et al., 2016). A '-search fGarau et al. , 2005), RRTs (Rao and Williams. 2009), kinematic tree-based navigation (Chakrabarty and Langelaan, 2013 ), stochastic planners with uncertain edge weights (Wellman et al, 1995), and stochastic surface response methods (Kewlani et al., 2009), Other techniques such as nonlinear optimization methods (Kru ger et al. 2007; Witt and Dunbabin, 20081, sequential quadratic programming (Bey!kin. 2008), evolutionary algorithms (Alvarez et al, 2(10 ! . Aghahaba, 2012), fast marching methods (Sethian, 1999; Petres et al,, 2007), wave front expansion (Sotilignac et aL. 2009; Thompson et al., 2010} have been used with varying degrees of success. Monte Carlo methods to account for uncertain ties and compute statistics of optimal trajectories have been used with potential field methods (Batraquand and Latombe. 1990) and two point boundary value problems (Wang et al. 2016). Partially Observable Markov Decision Processes have been utilized for computing policies of ro bots moving in uncertain environments (Amato et al,, 2016). We refer to Leraiuskm.s et al. 2016, 2017) and references therein for a detailed review of ail these results, and to Blythe ( 1999) for an early review and outlook of decision-theoretic planning.

Briefly, the main issue with several of the above methods is that they are inaccurate, or require application specific heuristics, or are computationally intractable strong and dynamic flows. Several methods do not rigorously utilize environmental predictions for planning. Additionally, accurate accounting of uncertainties requires large ensemble sizes for convergence of Monte Carlo methods, making them computationally very expensive for real-time use.

2. Theory and Schemes

2.1. Problem Statement

The risk-optima) path planning problem in uncertain flows can be formulated as follows (Fig. 4- 1)

Figure 4-1: Schematic of minimum-risk- time-optimal path pksmmg setup: The goal is to compute the time -optimal path with minimum risk under uncertainty for vehicles navigating from \ _s to Xf in an uncertain (low field V0r, f;<u) . The effective· velocity, U experienced by the vehicle is the vector sum of the vehicle's forward motion Fifihit} and the background 8«w V.

Consider a domain Ό with a spatial index x, temporal index t, probabilistic sample space W, random event w e W, and probability distribution function !¾(«}. For a vehicle Q navigating from start x _s to a target X _f in a stochastic flow v(x, t; *), we denote the time-optimal path distribution by X _jj t _f , t ·) and the time-optimal heading distribu-tion by Mr, ·}. We define risk-optimal paths as the time-optima! paths that minimize the expected utility cast of following a path that is possibly sub-optimal in the uncertain environment. Next, we develop the theory to compute such risk-optimal paths

2.2. Theory of Risk-Optimal Paths

Risk-optimal paths can be obtained in three steps. First, v(x, n ·) is computed by solving the DO stochastic environmental How equations (Appendix A.! ; Ueckermann et al, (2013); Suhra- mani (2018), Second, C _d (c _4· , i\ ·) and hit; ·} are computed by solving the stochastic dynamically orthogonal level -set Hamilton Jacobi equatio s (eq. A.7) followed by the stochastic backtracking equation (eq, A.5). We developed and applied these equations in Subramani et a!. (2018), Third, a risk metric is selected, the risk of all path choices Xg(x _s, t; ·} is evaluated, and the risk-optimal path XgiX , i) computed by minimizing the risk. Here, we focus on all aspects of this third step.

To evaluate and minimize the risk of following a wrong path from among the distributi on of time-optimal paths, we employ the expected utility hypothesis (e.g.. Arrow, 1958; Sehoe aker, 1982) The risk corresponding to a path choice is the expected utility cost of that choice being sub-optimal due to the uncertain environment The utility cost depends on the error incurred due to the sub-optimality and the cost a decision maker assigns to this error, in accordance with their risk profile or chosen risk metric. Mathematically, we can formulate the risk evaluation and minimization problem as follows.

The time-optimal path Cr(c _L r, a>j) or heading time-series h(t; w ) corresponding to a random event w _{/ €} W is exact for the environmental realization v(x _> ;<¾¾); however, its use in another environmental realization v(x./; o½) would be sub-optimal, and potentially infeasible. Let the trajectory achieved by following the waypoints or headings of the optimal path for &¾ in w„, be C ΐn, /; oyjw _f e ). We define the error incurred due to this path as e(Xg(s _s , t\ & _¾ !& _½) ) and the cost f\\ ₍ « i; - f(eiX _#{ x _s . r. to _/ j&½))) to quant i !y the effect of this sub-optimality on the decixu n maker. Here / is a utility cost function corresponding to the risk profile of the decision maker. Thus, the risk of utilizing the time-optimal path Xg(x,. /: uy) is the expected cost, i.e,

Then, the risk-optimal path is computed by optimization as

Xg(X _j , r) - arg min RiXgi ^ r; ny)) . (2)

Of course, there cart be many time -optimal paths that have the same risk RI(Cr( ^. ί: wi)) In eq (.1 ). Similarly, the risk-optimal path defined in eq. (2) can also be multi-valued.

2,2. Schemes for Risk Evaluation and Minimization

The specific choice of the error incurred due to sub-optimality and the form of utility cost function determine the risk evaluation and risk-optimal path selection. For brevity of notation, we use hereafter i _m for the error «(Cr(c.„ t w _/ !a _¼ )) and for the cost C(Xg(x _f , 1 wiIw _ίh )) There can he several error metrics corresponding to the operational parameters of interest and multiple analy tical forms of the utility cost function / corresponding to different risk tolerance profiles of different decision maker, e.g,. the pilot or mission designer (Schoemaker, 1982; LiCaizi and Sorato, 2006). Next, we present the schemes we consider.

2.3.1. Error Metric e _m -

The purpose of the error metric in our problem is to quantify the difference between following a prescribed path (specified as waypoints or a heading time-series) in an environment for which

 it is sab-optimal and following the true optimal path for that environment. We consider vehicle operations where the objectives are either achieving the optimal set of provided waypoints or the optimal heading time-series, in both operational modes, the optimal path can always be followed in the environment flow realization for which it is truly optimal, i.e., for which it reaches the target. X/ in fastest time. For a vehicle programmed to achieve waypoints, the ability to traverse a non-optimal path depends on whether the currents are always weaker than the vehicle's forward thrust or not. If yes , all prescribed paths are then feasible in finite time and di ffer from the optimal path mainly in total travel time. If not, then the vehicle does not have local control lability and may get adverted away from the desired path; it may not be able to reach xy at all and thus need to be aborted. For a vehicle programmed to achieve the optimal heading time- eries for a given environment, if it is operated in an environment for which these headings are sub-optimal., the vehicle would in general end up at u different end location than its intended target.

To quantify the physical dissimilarity between the paths Cr(c _> r, ie _¾ j < _½) and XgfXi, /; & _{½ }} , we employ the discrete Frechet distance (Eiter and annila, 1994; Alt and Godau, 1995) as the error metric. A zero value of the discrete Frechet distance means that the paths are coincident and higher values progressively imply more dissimilarity.

If the vehicle is programmed with waypoint objectives and is controllable, then Xg(x*, t. ) is simply X _Q (X _s . t\i _}) and the physical dissimilarity error metric can be directly computed from she distribution of optima! paths.

Alternative error metrics are possible for vehicles programmed with waypoint objectives or heading objectives. For example, for a vehicle programmed with waypoint objectives, the sub-optimal trajectory reaches the target in time 7( y;cy!c½) instead of T f.xy; c½). Thus, the error metric could he [77A7 << ₎ - T{x , a½)l. For a vehicle programmed with a heading ob _j ective, an appropriate error metric could be the error in arrival location. Specifically, follo wing \ _{0{ s} , t ά. _/ ί« _{ί )} , the vehicle ends up at xj- inst ad of x _f - and the error metric is the distance

It the vehicle is not controllable as it gets caught in an unexpected current or wind field stronger than the vehicle speed, in practice the ission may quickly get aborted. Here, the error metric could be set its maximum if the mission is aborted. If that current or wind field larger than the vehicle speed was forecast as possible, than the above metrics would continue to be employed.

All the above metrics have the property that if the path followed was the exact time-optimal path for the realized environment then the error is zero. For non-optimal paths, the error gets progressively higher with the degree of sub-optimality.

We note that, the above metrics have been chosen based on our experience with real time i sion with REMUS 600 All Vs fSubramani et ai., 2017b; Mirabito et a!., 2017; Edwards et ai., 2017). Of course, other error metrics are also possible. Examples include metrics based on the maximum background flow encountered along the path, or minimum bathymetry clearance along the path, or the number of times a vehicle has to perform a particular maneuver such as surfacing. Specific error metrics for individual missions and vehicles can be formulated based on an accurate performance characterization of the vehicle utilized, its navigation capabilities, fin configurations, and thresholds for satisfying navigation objectives. In the present paper, our goal is to provide a genera! theor and a recipe for risk-optimal path design and prediction for any acceptable function ti _m . If a specific metric is well adapted to a given situation, one would simply need to update the definition of the rest of our schemes would still apply. 2.3.2. Analytical Forms of the Cost Function f

The analytical form of the cost function / that translates the error to a cost is important to characterize the risk tolerance of vehicle operators. A concave cost function implies risk-seeking preferences for choices within the range of concavity, a convex cost function implies risk - seeking preferences for choices and a cost function with a constant slope is a risk-neutral preference for choices (Fishburn, 1979; Schoemaker, 1982; LiCalzi and Sorato, 2006; Von Neumann and Morgenstera, 2007). In our applications, we will employ /(<¾,} - logs 1 + t?; _m } for a risk seeking decision maker, = e _im for a risk neutral behavior, and/i¾) - - 1 ÷ exp(<¾) for a risk averse behavior.

24. Algorithm

The algorithm to implenient the theory and schemes developed in Sec. 2.2-2.3 has four major steps as summarized in Table 1. See Table B.2 for our notation. A subscript DO is used to indicate that a DO for KL) expansion and truncation has been applied to the corresponding stochastic vari able.

Table I: Ris -OpiifoaJ Path Planning: Algorithm

I. Probabilistic Flow Prediction

5. Predict the probability of the velocity field Vpoix, s; r) for r = 1 n _r by solving (he discrete stochastic DO harotropie quasi-geostrophk equations (etjs _< A.3) or discrete stochastic DO Primitive Equations (Chapter 2).

If, Stochastic T ie-Optmial Path Planning

1 Compute stochastic maximum icadiahiiifv fronts r) for r = 1 · · · n _T all at once by solving the stochastic

DO level -set equations (eq. A.7) with v»o(x,f; > }.

2. Compute discrete time-ctptiwEi l paths Xr>(Xs _< ri r) and time -optim ! headiogs hit; ?) for all r = I n _r using the backtracking equation: icq. A .5).

III. Risk Evaluation and Optimization

1. SinmkiSe the trajectories Cb(Ci, n htti) for waypoint objective Xg(x _s .r ) or heading objective ft(t; r) for till l = I · · ·«,. Under six- conditions of complete controllability nd for waypoint objectives. Xg(x. _f , P,i|tn) =

2. Compute the error metric matrix ¾« as the discrete Frecfsei distance between Xg(x _s ,¾;/i»i> and Xy(x _{( ,} /: «?) for all l, m— I · · · or one of the alternative error metrics defined in Sec. 2.3.1,

3. Compute the cost matrix C _>m = fi t¾) for all L m = I - · ¾ ,

4. Compute the risk of choosing the path XQ<X _s , t I) as B;— C*».

5. Find l that minimizes f¾. Then. foe risk optimal path XQ( _S ) = C _b (c _{L ,} f:l).

IV. Iterate

1. If ihete tv e\ n -nce that the number of DO realizations r = i · · n _r or adaptive sizes of Use DO envitonntentai or level-set · itbs .tces, ft . · or rtv.e (Appendix A ), ore not sufficient to describe the probability density of v(x, t *) or of X _jj t . r,·» ,ugment the B umber o realizations or subspaee sizes, and iterate. 3. Applications

In this section, we showcase the above theory and schemes by applying them to three environ mental flows of increasing complexity. In the first, we consider an autonomous vehicle crossing a simulated stochastic steady front tin the ocean or atmosphere) similar to the one we have used to previously tLolla et ai.. 2014b; Subramani and Lermusiaux, 2016; Subramani et a!.. 2018), but now for stochastic flows and risk-optimal paths. Here, we consider a hi -modal Gaussian Mix ture Model (GMM) PDF for the uncertain flow strength. In the second, we consider a stochastic quasi -geostrophic wind-driven double-gyre flow, as in Subramani et !. (2018) but now add the risk evaluation and optimization. In the third, we consider an autonomous vehicle navigating in simulated flow exiting a strait/estuary (in the coastal ocean) or wind blowing through a widening constriction of as urban canopy (such from a narrow street onto a wide street both lined with tall buildings).

3.1. Stochastic Steady Front with Uncertain Flow Strength

The first stochastic flow scenario we consider - a steady front with an uncertain flow strength described by a Gaussian Mixture Model PDF - is an idealization of missions where an AUV or glider crosses a front or channel. Fig. 4-2a shows the domain and Fig. 4-2b shows the PDF of the flow strength.

Figures 4-2: Domain and ike PDF of flow strength for the stochastic simulated from crossing test case: (a) In a s uare domain of nor.-dintcnsional side lengths 100 x [03, an idealized stochastic steady front is modeled as a zonal- jet with uncertain strength, flowing from west to east between y = 40 and y = 60. (b) The PDF of the flow strength is a Gaussi n Mixture Model with two Gaussians with non-dinwnsionai mean, standard deviation and mixture weight of ( H), 3, 0,65.» and <20, 1, 0.35) respectively.

The domain is a non-dimensional square basin of size 100 x 100 with an uncertain jet flowing from east to west and confined between 40 < y < 60. The direction of flow is considered known and only the strength is uncertain. The PDF of the flow strength is a Gaussian Mixture Model of two Gaussians components with non-dimensional mean, standard deviation and mixture weights of ( 10,3,0.65) and (20,1.0.35) respectively. A single _v = 1 mode is sufficient to describe _· the variability in the flow. We utilize ,<?,·,, = 5, 000 realizations sampled from the GMM distribution of the stochastic coefficient corresponding to mode 1. Here, we consider a non- symmetric PDF for the flow strength in order to illustrate the key aspects of our risk-optimal planning. This completes step I of Table. 1.

Next, we simulate the stochastic reachability and time-optimal path distribution for a vehicle with non-dimensional speed 20 starting at. (50,20) and heading to a target at (50,80). The start point is depicted by a circular marker and the target by a star marker in all our figures. Eq. A.7 is solved to obtain the distribution of stochastic reachability fronts. In Fig. 4-3a-e, five discrete snapshots during the evolution of the reachability front are shown. Fig. 4-3: (a)-(e) Stochastic reachability front evolution at five discrete times and (t) time-oplitnai path distribution for a vehicle navigating in a stochastic steady front with uncertain flow strength.

Each realization of the rea-chability from is colored by the jet strength of the corresponding flow realization. Fig. 4-3f shows the distribution of the time-optimal paths computed by solving eq. A.5, This completes step If of Table. 1 , Next, we compute the risk-optimal paths by completing step 111 of Table. 1.

3.1.1. Risk - Optimal Paths

Waypoint objective. For vehicles programmed with a waypoint objective and under conditions of complete controllability, no new trajectory simulation is required. We simply compute the error metric matrix as the dissimilarity matrix of all time-optimal path realizations. This t¾ matrix is of size n _r,4 x n _r and contains the pair-wise discrete Frechet distance between pairs (/, m) of paths for all /. m ~ 1 - n _r, . The utility cost matrix C*,, for the required cost-functions can then be computed from e _m . Fig, 4-4 row I shows the utility cost matrix for a risk-seeking cost function (column a), risk-neutral cost function (column b) and risk-averse cost function (column e).

Fig. 4-4: Computation of risk-optimal paths with waypoint objective for stochastic steady front crossing: Rows 1 ,2 and 3 correspond to step Hi.3, 111.4 sad Iii.5 respectively of Table, i . Columns a,b,e correspond to risk-seeking, risk-neutral and risk-averse behavior. To facilitate visualization, the waypoint objective choices l have been sorted, in rows 1 and 2. by the strength of the flow realization for which it is the exact time-optimal path.

The cost functions are computed by normalizing the errors to lie between 0 and 128 and applying the concave function fie ) = tog ₂ (t + e) to model risk-seeking behavior, the constant slope function / ( _<? > - e to mode! risk-neutral behavior and the convex functio f(e) - 2* to mode! risk-averse behavior. In the present range, the above functions exhibit concavity and convexity enough to model risk-behaviors (see, e.g., Schoemaker, .1 82). The /-axis of the cost matrix corresponds to the waypoint objective choice and the -axis corresponds to the realizations ol the flow field. The risk of the waypoint objective choices is then computed by marginalizing ou ^v i /« and obtaining the risk curve shown in row 2. The red point is the minimum risk choice. In row 3. each waypoint objective choice is colored by its risk, and the risk-optimal waypoint objective is shown in black.

Fig. 4-5a shows the three risk-optimal waypoint objective choices corresponding to different risk behaviors. Fig. 4-5b shows the PDF of the frechet distance between the risk-optimal waypoint objectives and the true-time optimal path in all realizations of the stochastic flow environment.

Fig. 4-5: Risk-optimal waypoint objectives far stochastic steady front crossing: (a) The risk-seeking, risk -neutral and risk-averse waypoint objective chokes (b) The PDF of error (quantity jag the degree of suh-optitna!ity) due to fallowing the rtsfc-optirnal paths. The ensx ts the discrete Frechet distance between the risk-optima! choice and the true time- opti ad potit ««responding to Use realized envifoartieni. The PDF of the error describes the properties of the risk-optimal paths succinctly. By following the risk- seeking path, we have a high probability of having low errors, taut also a high probability of having high errors. On the other hand, by following the risk-averse path we have hi itet certainty of making medium errors. The risk-neutral path has error characteristic in between the risk-seeking and risk-averse choices. Hence, depending on the risk appetite of the operator, one of the risk-optimal choices will be appropriate. An aggressive operator may want to bet. with the risk-seeking choice with the understanding that, for the realized environment, this choice has a high probability of being close to the true time-optimal path but also a high probability of being far from the true lime-optima! path. Similarly a conservative operator may choose the risk- averse choice and have lower probability of large errors but also for small errors. Ail ambivalent operator may choose the risk-neutral choice.

Gh Fig. 4-6, the errors due to following the risk-optimal choice are visualized on the other time-optimal paths, if they were the tree time-optimal path.

Fig. 4-6. Error visualization for risk-optimal waypoint objective choices: Each time-optima! path is colored by the Fnechet distance between it and the risk-optimal choice (a), (b) and (c) correspond to fee risk-seeking. risk-neutral and risk-averse choices. Tire RDF of these errors is shown in Fig, 4-5b.

This is the physical visualization of the errors presented as a PDF in Fig. 4-5h.

Heading objective. For a heading objective, we simulate the trajectories obtained by following a panic-alar heading objective choice in all the realizations of the Slow. T he e rror m atrix is then computed by the discrete frechet distance between the realized trajectory and the true time- optima! path for that flow realization. Matrices and risk curves similar to that shown in Fig 4-4 can be _· obtained for computing risk-optimal heading objective choices. We only show the risk-optimal heading choice and the error characteristics for these choices here. Fig. 4-7a shows the compute risk-optimal heading choices for our three risk behaviors Fig. 4-7b shows the PDF of the error due to following the risk-optimal heading choice.

Fig, 4-7: Risk-optimal heading objectives pr stochastic steady front crossing; (a) The risk-seeking. risk-nentrai and risk-averse heading objective choices, (b) The PDF of errors dire I» following the risk-optimal heading objectives. The error is quantised as the discrete Frechet distance between the path obtained by following the risk-optimal c hoice and fee true time-optima! path corresponding to that realized environmental flow.

The error is quantified as the discrete Frechet distance between the trajectory obtained due to following the risk-optimal heading in all flow realizations and the true time-optimal path for that flow realization. Fig. 4-8 shows this error visualized on the trajectories obtained by following the risk-optimal paths in all the flow realizations.

Fig. 4-8; Path distribution due to foi!a mg tke risk-optimal heading objectives in ail environmental fictw realizations: Each realize path corresponds to a particular flow realization and 1 s colored by the d iscrete Frechet distance between this path and the true time-optimal path for that realized environmental flow, (a), (b) and (c) show the errors due to following risk-seeking, risk-neutral and risk-averse choices. Fig. 4-?b shows fee PDF of these errors. From Fig. 4-8, we see that by following the risk-optimal heading choice in a flow for which it is sub-optimal results in the vehicle missing its target. The risk-seeking heading objective choice results in paths that are close to the true time-optimal path in more flow realizations that the risk-neutral and risk-averse choices. On the flip side, the risk-seeking choice also leads to paths that are tar fro the true-time optimal paths in more flow realizations than the other two choices. As before, the risk averse choice leads to more certain medium errors.

3.2. Stochastic Wind-Driven Double Gyre

As our second illustration, we consider a stochastic wind-driven double gym. This is the same flow fields that were presented in Subramani et al. (2018). This choice of test case is done to ensure continuity of our theoretical development. They are obtained by solving the quasi- geostrophie DO equations A, 3 with h _ί = 5 DO modes and - 5, 000 DO realizations. The domain is 1 ,000 km X LOGO km discretized into a 100 X 100 regular finite olume g rid. A deterministic zonal wind forcing is applied to drive the double gyre dynamics initialized wit uncertain barotropic velocity components. The flow field for 13.5 days for a mission with a vehicle moving at 40 cm/s .is considered. Fig. 4-9 shows the mean, standard deviation and skewness of the velocity fields at three discrete limes T = 0, 6.75 and 13,5 days.

Fig.4-9: Si ciuudc wind-driven double gyre flow: The mean, siaodatd deviation and skewness of ihe velocity fields ate shown ia the first three t ows. Peiturbstioa from the DO ateaaoftwo representative realisations are show» «» the last two rows. AO fields are shown for three discrete times at the he-gitming. middle and end of the planning horizon in each column.

The perturbations from the DO mean of Iwo extreme realizations are also shown. These two and the other 4,998 realizations are all simulated by one DO simulations with DO mean, mode and coefficient equations (see Subramani et al. (2018)). From the higher order statistics of the flow fields, we see that the flow exhibits non-Gaussian behavior. The extreme realizations have perturbations in the opposite direction with similar magnitudes showing that flow realizations with very different flow proprties can be simulated by a single DO realization.

3.2.1. Risk Optimal Paths

As for the previous illustration, here we compute the risk-optimal paths for both waypoint objectives and heading objectives.

kiypoim objective. Under conditions of complete controllability, we can compute the utility cost matrix directl from the distribution of time-optimal paths. Results from each step of the computation are presented in Fig. 4-10.

Fig. 4-10: Computation of risk-optimal paths jar waypoint abjeab'e in the stochastic double gyre: Rows 1 ,2 and 3 correspond to step III.3, 111.4 and .111 _* 5 respectively of Table. I . Coimnas a.b,c correspond so ride-seeking, risk-neutral and risk-averse behavior. To facilitate visualization, the waypoint objective choices l have been sorted, in rows I and 2, by (he velocity DO Coefficient 1. The rows 1 and 2 show the utility cost matrix and the risk curve. Row 3 shows all waypoint objective choices colored by their respective risk and the risk-optimal choice colored by black. The columns a,b and c correspond to the three risk behaviors: risk-seeking, risk -neutral and risk-averse respectively.

Fig. 4- 11a shows the three risk-optimal paths corresponding to different risk behaviors together. Fig. 4-1 l b shows the PDF of the Freehet distance between the risk-optimal waypoint objectives and the true-time optimal path in all realizations of the stochastic flow environment.

Fig. 4-11. Risk-optional waypoint objectives in the stochastic double gyre flow field: Same as Fig 4-5 but for the mission in stochastic double gyre flow.

The properties of the risk-optimal paths are similar to those in the previous illustration. This shows that our method can be utilized in complex stochastic flow scenarios efficiently. The information in the PDF of the error is visualized on the physical paths in Fig. 4-12.

Fig. 4-12: Error visualization for waypoint objective risk-optimal choices: Each time -optimal path is colored by the Freehet distance between it and ffae risk -optima! choice, (a), (b) and <c> correspond to the risk-seeking, risk-neatral and risk-averse choices. The PDF of these errors is shown in Fig. 4- lib.

Heading objective. The risk evaluation and minimization is completed for vehicles following heading objectives. Fig. 13a shows the risk-optimal heading objectives and Fig. 4- 13b shows the error PDF due to following the risk-optimal headings in all flow realizations.

Fig. 4-13: Risk-optimal heading objectives in the stochastic double gyre flow field: (a) The risk-seeking, risk-neatral and risk-averse heading objective choices (b) The PDF of errors doe to following the risk-optimal heading objectives. The error is quantified as i he discrete Freehet distance between it® path obtained by following the risk-optimal choice and the true rime-optima! path corresponding to that realized environmental slow.

In Fig. 4-14, we show the trajectories obtained by following the risk-optimal headings in all the flow realizations. Here, the paths are colored by the Freehet distance between the obtained trajectory and the true time-optimal path in that flow' realization.

Fig. 4-14: Path distribution dm to following the risk-optimal heading objectives: Each realized path consspcmds to a particular flow realization and i s colored by the discrete Freehet distance between this path and tire tree time -optimal path for that realized environmental flow (a), (b) and (c) show the errors doe to following risk-reeking, risk-neutral aad risk-averse choices. Fig. 4- 13b shows the PDF of these errors. 3.3. Stochastic Flow Exiting a Strait

As our third illustration, we present risk-optimal planning is an idealized stochastic flow scenario encountered in the coastal oceans and in urban environments. In the coastal oceans, a barotropic jet exiting a strait (or an estuary) into a wider channel creates eddies and meanders downstream (Lolla et al, 2015). In urban environments, the wind blowing through narrow con strictions between buildings into an open area also creates similar dynamics. Such flows can be idealized as sudden expansion flows studied extensively i n fluid dynamics (Cherdron et al., 1978; Durst et al., 1974; Feam et al., 1990). We consider a 6 km x 1. kra channel with a narrow constriction of length 1/3 kra on the West. A uniform jet with a velocity 50 cm/s is exiting the constriction into the channel The initial conditions of the barotropic velocity components in the channel are uncertain and sampled from a Gaussian covariance kernel with a decorrelation length scale of .1 km. The simulation is conducted by solving the quasi -geostrophie DO equations A, 3 with /i. _{S ¥} = 10 DC) modes and n _r.v = 10, 000 DO realizations. The stochastic flow is allowed to develop for 1 ,500 mins by which time the flow develops recirculation zones and breaks to either the north or south of the centerline depending on the initial uncertain perturbations. The nonli near dynamics causes the initial gaussian uncertainty to become non-Gaussian by this time. For our risk-optimal planning, we consider the flow for the next 100 mins. Fig. 4-15 shows the mean field, and the first three DO velocity modes and the marginal PDF of the DO velocity coefficients at three diereie times T-0, 50 and 100 mins.

Fig. 4— 15: Stochastic ftmv exiting a strait: Mean. Mode#!.#2,#3, CoefF. #t ,#2,#3 are shown at three discrete times at the beginning, mid le and end of die planning horizon. The mean and modes velocity streamlines are overlaid on a co!orpiot of their magnitude. The coefficients are shown by their marginal PDF.

As can be seen here, the uncertainty in the flow field is highly dynamic and non- Gaussian, but the mean flow is nearly steady. Such flows areoften encountered in the coastal ocean or atmosphere. Fig. 4- 16 shows the standard deviation, skewness and kiutosis of die velocity field at three discrete times.

Fig. 4-16: Statistics and representative realizations of the stochastic flow exiting a strait: The standard deviation, skewness and kur!osis of the velocity lie ids are shown in the first three rows. Two representative realizations are shown in the last two rows. All fields are shown for three discrete limes at tire beginning, middle and end of the planning horizon in each column.

Also shown are two representative realizations with the jet breaking to the south and north. This completes step I of Table, i.

3.3 J. Stochastic Reachability and Time-Optimal Path Distribution

Our path planning problem is to predict the risk optimal paths of a vehicle with a nominal relative speed of 25 cm/s traveling from the start point (depicted by a circle in Fig. 4- 17) to five tar-get locations (depicted by slars in Fig. 4- 17).

Fig.4- 17; Stochastic reachability front evolution colored by DO velocity coeff. # 1. The reachability fronts are computed by one DO simulation by solving the stochastic DO level set equations with the above stochastic DO velocity fields for the stochastic flow exiting a strait. Using the velocity DO mean, modes and coefficients computed above, we solve the stochastic DO level set equations A 7 to obtain the distribution of the stochastic reachability fronts. Fig, 4-17 shows the evolution of the stochastic reachability front at six discrete times, T- 16, 33, 50, 66, 83 and l(K> mins, in each panel, the reachability fronts are colored by the DO velocity coefficient I of the stochastic flow field. The time-optimal path distribution is computed from the stochastic reachability fronts and backtracking eq. A .5. Fig. 4- 18 shows the distribution of the time-optimal paths to the five targets, each colored with the arrival time at the target.

Fig 4- ! 8; Ttme-opiimal path distribution colored by arrival time in mins The distribution of exact time-optima] paths are confuted from Ox stochastic &acfeabiJity frosits using die stochastic back tracking equation.

The paths in panels (a) and (e) have low physical variability but significant arrival time variability. The paths in panels ic) have high physical variability but low variability in the arrival time and the paths in panels (b) and (c) have high variability in both the physical paths and arrival times. This completes step II of Table. I.

3.3.2. Risk Optima] Paths

For each of the live targets, we complete step III of Table i to obtain the risk optimal paths (both waypoint objectives and heading objectives separately).

Waypoint Objectives.. Here, as before, under assumptions of complete controllability, the error metric can be computed directly from the distribution of the tim -optimal paths. Fig. 4-19 column 1 shows the competed risk-optimal paths that minimize risk-seeking, risk-neutral and risk-averse cost functions for all the five targets (rows (a) to (e)>. Fig. 4-19 column 2 shows the PDF of the errors due to following the computed risk-optimal waypoint objectives corresponding to the paths in column 1.

Fig. 4-19: Risk-optimal paths with waypoint objectives in the stochastic jUyw exiting a strait: (a)-(e) correspond to the live target locations. Column 1 shows the risk-seeking, risk-neutral and risk-averse· waypoint objective choices. Column 2 shows die PDF of errors due so following the risk optimal paths. The error is quantified as she discrete Frechet distance between tire risk-optima! choice and the trne lime-optimal path corresponding to the realized environment.

Due to the complex and dynamic nature of the path uncertainty, the different risk-optimal paths show interesting behavior For targets a and e, the risk-neutral and risk-seeking paths are physically very close to each other and have similar error PDFs. However, for targets b, c and d, it is the risk-neutral and risk -averse paths that are physically closer. For all targets, the error PDF of the the risk-seeking and risk-averse paths are similar to the error PDF in other test cases. Crucially, our rigorous probabilistic prediction and risk-optimal path planning framework allows computation of such complex risk-optimal paths very efficiently.

Heading Objectives.. Fig. 4-20 shows the risk-optimal heading objectives (in column 1 ) and the error PDF (column 2) .for the five target locations (rows (a.) to (e)), computed by completing step ill of Table 1. Fig, 4-20: Risk-optimal heading objectives in the. stochastic flow exiting a strait: Rows (ti -(e) correspond to the five targets. Column 1 shows the risk-seeking, risk-neutral and risk-averse heading objective chokes, and column 2 shows the PDF of errors due to following itiv n k-optifna! heading objectives. The error is quantified as the discrete Frechet distance between the path obtained by foliowing the risk-optimal choice and the true time-opthn&! path corresponding to ihsit realized environroenta! flow.

Fig. 4-21 shows the trajectories obtained by following the risk-optimal heading objectives in all the flow realizations. As before, each trajectory is colored by the discrete Frechet distance between that and the true time-optimal path for that flow realization.

Fig. 4-21 : Trajectories obtained b fltilowing risk-optima! headings in the stochastic flow exiting a strait: saHe.) corns- «pond to the five target locations. The cohmias correspond to risk-seeking, risk-neutrai and risk-averse waypoint objective choices. Each trajectoty is colored by the discrete Frechet distance between that trajectory and the exact time-optima! path corresponding to that realization of the envitonineal.

As before, each trajectory is colored by the discrete Frechet distance between that and the true time-optimal path for chat flow realization.

4. Conclusion

We developed the theory for rational risk -optimal path planning by combining decision the ory and our stochastic time-optimal path planning with stochastic DO level-set equations. The schemes and software developed compute risk-optimal paths for vehicles navigating in uncertain, strong and dynamic flows. The path planning proceeds in three steps: (i) obtain predictions of die probability distribution of environmental flows (ii) obtain predictions of the distribution of exact time-optimal paths for the above flow distribution, and {Hi} compute and minimize the risk of following the above time-optimal paths. Three cost functions corresponding to risk-seeki ng, risk-neutral and risk-averse behaviors were utilized to compute and minimize risk. We illustrated foe planning in three stochastic flow scenarios: stochastic steady front crossing, double gyre and flow exiting a strait. The risk-optima! paths minimize the error of following the chosen path, if it is not the exact time-optimal path for that environmental flow realization. Minimizing a risk- seeking cost function results in a pat that has a higher probability of being closer to the exact time-optimal path, but also a higher probability of being far away from the exact time- optimal path. A risk-averse path on the other hand has a high probability of medium error, and a low probability of either extremes (being very close or very far from the exact time optimal path). In complex flow situations, it is difficult to predict the behavior a-priori as shown in the sto chastic flow exiting a strait example. Here, our methodology allows in the characterization and minimization of the risks along time-optimal paths. In the future, risk-optimal planning can be integrated with probabilistic current predictions for real-time missions with real AUVs.

Appendix A. Summary of Stochastic Time-Optimal Path Planning Appendix A.l. Stochastic Dynamically Orthogonal Barotropic Quasi-Geostrophic Equations In this information, we briefly provide the stochastic PDEs that govern the dynamics of uncertain

Barotropic Quasi-Geostophic flows. For Specifics, we refer to Ueckermann et al. (2013). The stochastic Barotropic Quasi-Geostophic equation, written as conservation of momentum in Langevin form is

where Re is the Reynold’s number, / is the Coriolis coefficient, an a the strength of the wind stress. For the Coriolis coefficient, we employ a/CpIane approximation / = /¾ + bn.

Introducing a DO decomposition to n(c,G, ), where El pi o) j ~ 0 V j = i - · we can derive the DO mean (eqs. A Ja,A.3b), coefficient ieq, A.3c) and mode equations (eqs, A 3d.A.3e) of the S-PDE eq A,i as

where we have dr pp' the prrmnffiesfe (x, ft ed for brevity of notation. Here, E is the expectation o erator, C is the cmanufK matrix, M is the third moment, sad is the spatial inner- product operator for any t o iieMs # and ¾A defined as #o _? } * Appendix A.2. Stochastic Dynamically Orthogonal Level- set Equations

Here, we summarize the exact stochastic level-set equations and the efficient stochastic DO level-set equations derived for stochastic time-optimal path planning in Subramani et al. (2018).

The stochastic reachability .front for a vekide-Q (Fig. 1.) is governed by the stochastic Hamiitoa- Ja¾ · c >bi ( H I ; level-set equation

with the initial condition <Mx, 0; w) ~ Ix-x and, as needed, open boundary condition _tVi·}

0, where n is the outward normal to the boundary 0£>. w is a random event, an f is the reachabiltty-from-irad ng scalar level-sel tiekl <e.g . signed distance function) For every w. the optima] arrival-time Tixp.o at x _f is obtained by integrating eq A 4 until the first iirne t such that d>( _f wϊ < 0. The optima! trajectory Xpt ,,/; w) is then given by the particle backtracking equation

To arrive at the DO level-set equations, let us introduce a DO decomposition to the stochastic level-set F as

where Ef fgr; «i)j = 0 V / = 1 - · ^.

1¾e stochastic DO level-set equations are

Appendix B. Notation

The notation for variables used in the present paper is furnished here in Table B.2,

Agbahaba, M, R, 2012, 3d path planning for underwater vehicles using five evolutionary optimization sigcrith-Hs avoi ding static and energetic obstacles, Applied Ocean Research 3&, 48-42.

Alt H., Godau. M„ 1995. Conning the Mchst distance between two polygonal curves International Journal of Com putational Geometry & Applications 5 (Oln02>, 75-91.

Alvarez, A., Caiti, A., Onfcen, R„ 2004 Evolutionary path planning for autonomous underwater vehicles in a variable ocean. IEEE Journal of Oceanic Engineering 29 (2), 418-429

Amato, C, Konidaris, G , Anders, A., Cruz, G , How. 3. R, Ka blmg, L, P., 2016 Policy searjeh for multi-robot eoordi- aalton under uncertainly, The iatematioRal Journal of Robotics Research 35 04), 1764 77S,

Table B,2: Relevant notation.

Scalar

f e 14 Stochastic sabspaeii isdex

Q Vehicle whose risk-op tioiai path is to be determined

F e l Vehicle speed

,*! _¾ · g Dime ask® of ft ihe physics! domain

i , e ?1 Number of discrete grid. points in the ,c -direct ton

Hy € A Number of discrete grid points is the y-dkeciioa

>!.. € A Total number of discrete grid points in

h ,, _{L €} A Dimension of the stochastic subspace of A the level-set scalar field

e ¾ Dimension of the stochastic subspace of v, the velocity vector field

Yi €¾ Random variable describing the PDF of the orthonormal level-set (f) modes in

Pi € JR Random variable describing the PDF of the orthononaal velocity (?) modes, v,· n, e N Number of Monte Carlo realizations of the level-sets (and velocity vector fields) r e N Realization index

Velocity vector field

Mean velocity vector field

DO mode i of ^’ v: Dynamically orthogonal basis for the stochastic sabsp55ce of v

Level-Set field

Mean Level-Set hel

DO mode i of 4: Dynamically ^~ basis for the: stochastic

Arrow, K. 1 . 1958. Utilities, atsiaides, choices: A review m>te. Econometrics: louraal of the Econometric Society. 1-23.

Barraguand. J.. Lsto lse. 3 -C. 1990. A morae-carlo algorithm for path planning with many degrees of freedom. In: Robotics and Automation, 199(1, Proceedings.. 1990 IEEE International Conference on. IEEE, pp. 1712-1717.

Bc lkin. D„ 2098. Path optimkahon for an earth-based demonstration balloon flight. Surf report, Caltech. F¼sadena, CA.

Blythe, j., 1999. An Overview of Planning Untter Uncertainly. S ring Berlin Heideibeqe, Berlin, Heidelberg, pp. 85- 110.

Boddfy, M., Dean. T. L, _» 1 89. Solving time-dependent planning problems. Brown University, Department of Computer Science, Providence, Rf.

Bty. A . Roy, N.. 20.1 1. Raptdly-explonmg random belief trees for motion planning under uncertainty, la: Robotics and Automaton tlCRAi. 203 ί ΪEEE International Conference on. IEEE, pp 723-730.

Burns. B.. Brack, O., 2007, Sampling-based motion planning with sensing uncertainty. in: Rohotks and Automation, 2007 lEEE lnieroationsl Conference on. IEEE, pp. 3313-3318.

Cbakrsbartv. A., Langehtan. J., 2013. UAV flight padi planning in time varying complex wind- ields in: American Control Conference. IEEE, pp. 2568-2574.

Chsrdron. W„ Durst. E. Whilelaw. 1. H., >978. Asymmetric flows and instabilities in symmetric duets with sudden expansions. Journal of Fluid Mechanics 84 (i), 13-31,

Dursi, F„ Meiling, A,, Whilelaw. i. H, . 1974, Low Reynolds number Sow over a plane symmetric sodden expansion. Journal of Fluid Mechanics 64, 1 f i-! 28.

Edwards, I. Smith, l, Girard, A.. W kman. D , Suferamaal D. R, Kulkami. C. 5., Haley, Jr.. P. i„ Mint*». C.. iw, S,, Lentajstaux, R F. J„ Jun. 2017, Data-driven learning and modeling of AUV operational characteitst'cs for optimal path planning. in; Oceans Ί7 MTS/IEEE Conference. Aberdeea

Eiter, T., Mannila. R, 1994, Computing discrete frdehef. distance, Tech, rep,, IbchnEehe Univeisitat Wien,

Epstein. L G, 1992. Behavior undernsk: Recent developments in theory and applications. Advances in economic theory 2 1-61

Fearn.8 M. Muilin. T . Clifte, K. A., 1990. Nonlinear flow phenomena in a symmetric sudden expansion. Journal of Fluid Mechanics 23 1. 595-608,

Fe pon „ Lonnosiaux. R E J. 2018. A seomeSic approach to dynamical model-order redaction. SIAM Journal oa Matrix Analysis sa Applicationsln press.

Fishbura. P. C.. i 979. Oft the nature of expected utility. In: Expected utility hypotheses and the Atiais paradox. Springer, pp 243-257

Fishburn, R C., 3988 Nonlinear preference and utility theory. Vol. 5. Johns Hopkins University Press, Baltimore.

Gams, B., Alvarez, A.. Oliver. G., 2005. Path planning of autonomous amferanter vehicles in current fields with complex spatial variability: an A* approach. In: Proceedings of the 2095 IEEE International Couference on Robotics and

Automati n. 1LL , rr F>4_ Ki ^, ¾

ffwsng. k.. Ahui . N.. I49G < i?rv, mo* ion planning - a ur\ev ACM Computing Surveys (CSt R> 24 (3), 219-291.

Kewlant, G., Ishigami, G , iagnemma. K., Oct 2(S) . Stochastic mobility-based path planning in uncertain environments. In: 2909 fEEE- ^' HSJ Ihtemuu l onference? on intelligent Robots and Systems. pp. 1 183-1189.

Kruger, D., Stoikm, R„ Blum. A.. Briganti. J.. aprii 2007. Optimal AUV path planning to extended missions in complex, fast-flow ing estuarine etn uviuitt m». In: Robotics and Automation. 2007 IEEE International Conference ear. pp. 4265- 4270.

La&trobe, J.-C., 2012. Robot Motion Planning. Vol. 124. Springer Science & Business Media, New York

La Valle. S. M.. 2006. Planning algorithms. Cambridge university press, Cambridge, UK.

Lermustaux, P. F. I, Lo!te, I, Haley, Jr., R J. , Ytgil, K.. Uecfcermatm, M. P., Sotiftergaard. T., Leslie, W. G , 2016. Screnco of autonomy; Time-optimal path planmng and adaptive sampling for swarms of mean vehicles. In; Curtin, T. i Ed.}. Springer Handbook of Ocean Engineering; Autonomous Ocean Vehicles, Subsystems and Conteol. Springer, Ch 21. pp 481 -498

Lermtisiaux, P. F j., Subrsmani, D. R, Lm, L, Kuikttrtu. C S., Gupta. A., Daft, A., Lolla, T., Haley. Jr . P. J. AH. W. R, Mirabito. C.. Jana. S., 2017 A futme for intelligent autonomous ocean observing systems. Journal of Marine .Research 75, special issue: The Science of Ocea Prediction. vol. 17 of The Sea. .In press

LrCalzi, M., Sorato, A..2006. The pearsan system of utility functions. Bat an Journal of Operational Research 172(2), 573.

Loila, T„ Haley. Jr,. P. J.„ Lermusi&tix. P. F J, _« 2014a. Time -optima! path planning in dynamic flows using level set equations: Realistic applications. Ocean Dynamics 6 1 10), 1399-1417.

Lolla. T., Haley. Jr.. R J., Lermusiaua. P. E. J., 2015. Path planning in multiscale ocean flows: coord motion and dynamic h titcles. ( Lean Mocicfoiic 04, 46-66.

Lolla. T.. LcrmcsKUiA. P F J. Ueckertnann. M P.. Haley. Jr.. P J. 2014b Time -optimal path lantting in dynamic Hows using level set equations: Theory and schemes. Ocean Dynamics 64 ( 10), 1373- 1397.,

Mannarini, G.. Pinardi. R, Coppisi. G„ Gdde, P. lafrati. A.. 2016. Visir-i: sma’ 1 vessels-teasl-time nautical routes using wave forecasts. Geoscjentittc Model Development 9 {4t. 1597-1625,

Mirabito, C. Siibramtmi. D. N., Loila. X, Haley. Jr., P. J., Jain, A,, Lermosiaux. P. F. J.. Li, C. Yes, D. K P.. Liu, Y, Hover, F. S., Fulsone, R, Edwards, I , Railey. K. E., Shaw. G., Jun.2017. Autonomy for surface ship interception. In: Oceans 37 MTS/ IEEE Conference. Aberdeen.

Pereira, A. A., Binney, J., Hollinger, G. A , Su hatme, G. S., 2013. Risk-aware path planning for antoaomoas underwater vehicles using predictive ocean, models. Journal of Field Robotics 30 <5j, 741-762.

Pete, C„ Pai!has, Y, Patron, P., Petillot, Y., Evans, J., Lane, i>., April 2007. Path planning for autonomous underwater vehicles. Robotics, IEEE Transactions on 23 (2), 331 -341.

Rao, D., Williams, S B., December 2-4 2009. Large-scale path planning to underwater gliders in ocean ciarents. In: Proceedings of Australasian Conference on .Robotics and Automation, pp. 28-35.

Rosenblatt, J. K., 2000. Optimal selection of uncertain actions by maximizing expected utility. Autonomous Robots 9 (l).

17-25.

Sap&is, T: P., Lennusiaux, P F ^j . Dec. 2009. Dynamically orthos JS itefo equations to continuou stochastic d na ical systems. Physica D: Nonlinear Phenomena 238 (23-24), ^’ 234 ^” —

Schoemaker. P. ^j .. 1982. The expected entity model: Its variants, purposes, evidence andlimilatkais. Journal of economic literature. 529-563.

SetlsisR, J. A., Jan. 1999. Past marching meihods. SLAM Rev. 41 (2), 199-235.

Souhgnac. M . Tuilliix-rt, R, Higher. M. ax 2009. Time- minimal pa& planning in dynamic cutsenl fields. In; IEEE international Conference on Robotics and Automation. pp. 2473 -2479.

Subramani D. N.. Feb. 2018, Probabilistic regional ocean ptetitaRms: Stochastk fields and optimal planning. Pfa.D, thesis, Massachusetts institute of Technology, Department of Mechanical Engineering. Cambridge. Massachusetts.

Suhramtm . D. N.. Haley. Jr,. P, J.. Lemtudaux, R E J., 2017a, Energy-optimal path planning in the coastal «m Journal of Geophysical Research: Oceans 122, 3981-4003.

Subramani. D . I-env.u^au R G. . 2016. Energy-optimal path planning by stochastic dynamically orthogonal level· set optimization. Ocean Modeling 100. 57-77.

Sohramani. D .V , Lermu sux P. FI .1 . Haley. Jr.. P. J., Mirshilo, Jana, S,, Kt kami. C. S., Girard, A., Wiekroait. ., Edwards. J.. Smith. 1. Jun. 2ttl?h. Time-optimal path planning: Real-time sea exercises, in: Oceans T7 MTS/ 1E Ccmfeieace. .Aberdeen.

Subramani. D. N , Leila. T., Haley. Jr., R J.. Lennussaus. P. F. I. 2015. A stochastic optimization mefitod for energy- based path planning. In: Ravels. S., Sandu. A. t&dsJ. DyDESS 2014. Vol 8964 ot LNCS. Springer, pp. 1-12

Subramani. D. N., Wei . J., Lermusisux. P. F. J..2018. Stochastic tixne-optimai path-planning in uncertain, strong, and dynamic Sows. Computer Methods in Applied Mechanics and Engineering 333, 218-237.

Sun. W.. Tsiotras. R. Lolla. T. Subeamam. D. N.. Lermusianx.. P. F. J.. Apr. 2017a. Malt^ te-pursner-txne-evader pursuit evasion game m dynamic flow fields. Journal of Guidance. Control and Dynamics 40 (7 L

Sun. W„ Tsiotras. R. Lolls, T.. Subramani, D, K, Lermasiaux, P. H I„ 2017b, Pursuit-evasion games in d nam c f ow fields via reachability set analysis. 1B: 2017 American Control Conference (ACC), Seattle pp. 4595-4600.

Thornpw!i. D. R.. Chum. S. A . Chao. V Li P . Cahsil. B . Levin. .1.. Sch ttcki. . Bal.ixunya. A. R. PetilSo. S., Arroti, M., Meisinger, M., 2010. Spat otemporal path planning in strong, dynamic.. i certiuu currents. In: Ptxsceedinga of IEEE International Confetence on Robotics and Automation, pp. 4778-4783.

Uecketxnatm. M. P. LennuMaux. P F J.. Sapsts·. T. Ft. Jan, 20 i 3. Numerici'I sx ornet for dynamically orthogonal equa tions of stochastic ftek! and ocean flows, Journal of Computational Physics 233, 272-294.

Vox· Neumann. J.. Moxcen t re, (>., 2!X¾ ⁷ . Theort. < _> rc3tne and economic behaxxor. Pnnce ton university press, Ptinceion, N1

Wang, T. Lr Mailtv. f J. P.. Hotext. L. Ki o. O M., 2016. Path planning in uncertain flow fields using ensemble method. Ocean D namics 66 « 10\ 1231-1251.

Wellman, M. P.. Ford. M., Larson, K.. 1995. Pat ^j t planning undu hhk ,■ s ixieni uncertainly. In: Proceedings of the Eleventh Conference os Uncettanm in ArMicuJ Intelligence t VI M rgan Kaufmaan Publishers lac., an Fras- casco, CA, USA, pp. 532-539.

URL http : //dl„ acm . org/citation . of nCidA20?4iS8. 2074218

Witt, L, Dtmbabin, M., December 2008. Go with the flow: Optimal any path planning in coastal sttviromneots. In: Kim, J., Mahony. R (Eds. *, Proceedings of Australasian Otmfeteace o® Robotics and Automation, pp. 86-94.

V. Additional Considerations

The terms“program” or“software” are used herein in a generic sense to refer to any type of computer code or set of processor-executable instructions that can be employed to program a computer or other processor to implement various aspects of embodiments as discussed above. Additionally, it should be appreciated that according to one aspect, one or more computer programs that when executed perform methods of the disclosure provided herein need not reside on a single computer or processor, but may be distributed in a modular fashion among different computers or processors to implement various aspects of the disclosure provided herein. Processor-executable instructions may be in many forms, such as program modules, executed by one or more computers or other devices. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform particular tasks or implement particular abstract data types. Typically, the functionality of the program modules may be combined or distributed as desired in various embodiments.

Also, data structures may be stored in one or more non-transitory computer-readable storage media in any suitable form. For simplicity of illustration, data structures may be shown to have fields that are related through location in the data structure. Such relationships may likewise be achieved by assigning storage for the fields with locations in a non-transitory computer-readable medium that convey relationship between the fields. However, any suitable mechanism may be used to establish relationships among information in fields of a data structure, including through the use of pointers, tags or other mechanisms that establish relationships among data elements.

Also, various inventive concepts may be embodied as one or more processes, of which examples (e.g., FIG. 4) has been provided. The acts performed as part of each process may be ordered in any suitable way. Accordingly, embodiments may be constructed in which acts are performed in an order different than illustrated, which may include performing some acts simultaneously, even though shown as sequential acts in illustrative embodiments.

All definitions, as defined and used herein, should be understood to control over dictionary definitions, and/or ordinary meanings of the defined terms.

Use of ordinal terms such as“first,”“second,”“third,” etc., in the claims to modify a claim element does not by itself connote any priority, precedence, or order of one claim element over another or the temporal order in which acts of a method are performed. Such terms are used merely as labels to distinguish one claim element having a certain name from another element having a same name (but for use of the ordinal term).

The phraseology and terminology used herein is for the purpose of description and should not be regarded as limiting. The use of "including," "comprising," "having," “containing”,“involving”, and variations thereof, is meant to encompass the items listed thereafter and additional items.

Having described several embodiments of the techniques described herein in detail, various modifications, and improvements will readily occur to those skilled in the art. Such modifications and improvements are intended to be within the spirit and scope of the disclosure. Accordingly, the foregoing description is by way of example only, and is not intended as limiting. The techniques are limited only as defined by the following claims and the equivalents thereto.

What is claimed is: