Stochastic Optimal Control of a Sailboat

In recent years, advances in numerical and mathematical methods have greatly impacted the sport of sailing. Yacht path planning algorithms are now heavily employed both as a part of match simulation software, and as a tool to provide live guidance on optimal decision-making during a race. Here, we focus on the problem of finding the optimal feedback control to perform live path planning in the presence of a stochastically-evolving weather system.

Weather conditions greatly impact the performance capabilities of a sailboat in a race. While it is possible to produce meteorological forecasts, there is an inevitable stochastic component to how weather conditions evolve during the course of a race which is impossible to predict exactly. It is critical that algorithms used for decision making correctly factor this randomness into their control policy. In [2], this randomness is modeled as a discrete-time Markov chain process which updates weather conditions at a fixed tick rate. In [3] this idea is extended to a continuous-time Markov chain, formulating the system as a free boundary problem whose solution is a switching curve giving the optimal control policy. Here, we will take the approach of [1] and describe optimality through dynamic programming, yielding a system of two coupled quasi-variational inequalities of Hamilton-Jacobi-Bellman (HJB) type with degenerate ellipticity.

For full physical accuracy, a model would have to include both the position and velocity of the sailboat in the state space, and model the forces on the sailboat in different state configurations, leading to at best $5$ coupled first-order equations in $5$ state variables. An approach along these lines can be found in [2]. But in long course matches, inertial effects associated with accelerating up to the equilibrium speed are (in general) secondary, allowing one to reduce the system to $3$ coupled first-order equations by assuming that the velocity vector can be changed instantaneously.

However, this assumption cannot be made in modeling the process of tacking, which is necessary when attempting to travel upwind. To make progress, sailors must travel at angles to the wind, repeatedly switching directions to steer in a zig-zag pattern on ether side of the wind direction. Each time the bow is swung to the other side of the wind (a single tack), the boat is significantly slowed down for the period of time in which the wind pushes against the boat. In previous works ([1],[4]), this is incorporated by modeling the system as a hybrid control problem, with a continuous control corresponding to the angle of the sails and a discrete control corresponding to which tack we should be on. The lost inertia has been modeled as a fixed switching delay, reflecting the time lost due to the reduced speed from performing the tack.

In this work, we extend and improve upon the methods of Ferretti and Festa [1], who solve the problem using a semi-Lagrangian discretization. In Sec. II, we give the original problem statement and recast it to a reduced coordinate system in 2D, which alone greatly reduces the time needed to solve the problem numerically. In Sec. III, we set up the optimal control problem following the structure of Ref. [1], and provide an improvement to their switching operator to incorporate the wind evolution during the process of a tack. In Sec. IV, we describe our semi-Lagrangian numerical methods, including a discretization in which the system is nearly-causal and can be solved in a handful of iterations. We present the results of numerical simulations in Sec. V, producing optimal switching maps and highlighting the effect of our improved switching operator. We conclude by listing directions for future work in Sec. VI.

We consider a problem introduced in [1]: time-optimal sailboat navigation to a closed target set $\mathcal{T}$ when the wind direction evolves stochastically. For simplicity, we focus on their simplified model where the wind speed is fixed, and the angle $\phi$ of the upwind direction (measured counterclockwise from the $y$-axis) evolves under a drift-diffusion process:

where $a$ is the drift rate, $\sigma$ is the strength of the Brownian motion, and $dW$ is the differential of a Wiener process.

We control the angle of the boat’s velocity relative to the upwind direction. But at any given time we are also restricted to an interval of angles determined by the sailboat’s current tack, or configuration of sails. We represent the tack as a controlled discrete variable $q\in\{1,2\}$; $q=1$ represents the starboard tack, which we model as having access to the angles $[0,\pi]$, while the port tack $q=2$ has access to $[-\pi,0]$. For notational convenience, we will use $(-1)^{q}u$ to encode the velocity angle relative to the upwind direction, with the steering control variable $u\in[0,\pi].$ So, the steering angle $u$ is measured counterclockwise for $q=1$ and clockwise for $q=2,$ with $u=0$ corresponding to motion directly against the wind. Due to our inertia-less approximation, whenever $u$ changes, the boat is assumed to instantly start travelling in the new direction with the angle-dependent speed $s(u)$. The speed profile (termed the polar plot) capturing this dependence is determined by the exact geometry of a specific boat. A typical example [5] is provided in Fig. 1(a).

Assuming that the boat’s location is $(x(t),y(t)),$ the system dynamics in absolute coordinates is then

We now simplify by focusing on radially symmetric problems: if we assume that the target $\mathcal{T}$ is circular, we no longer care about an absolute angle of the boat’s position relative to the target. Shifting to a coordinate system centered at $\mathcal{T},$ we can reduce the $(x,y,\phi)$ state space to just two coordinates $(r,\theta)$, where $r=\sqrt{x^{2}+y^{2}}$ is the radial distance of the boat from the center of $\mathcal{T}$ and $\theta=\phi+\operatorname{atan2}(-x,-y)$ is the angle of the wind measured counterclockwise relative to the straight line between the boat and the target; see Fig. 1(b). In this reduced representation, the dynamics are governed by

which describe the projections of the velocity vector onto the radial and angular components, combined with the drift-diffusion evolution of the wind. We refer to the coefficients on $dt$ in each equation as $\dot{r}$ and $\dot{\theta}$, respectively. New from the original system is the $\dot{\theta}\propto\frac{1}{r}$ relationship, which will force us to utilize an adaptive timestepping scheme that chooses smaller timesteps near the target to counteract the rapidly increasing angular speed. We note that $\dot{\theta}$ still remains bounded since $\mathcal{T}$ has positive radius.

Assuming we start from a position $\boldsymbol{r}=(r,\theta)$ and tack $q$, the total time-to-target $T_{\mu}(\boldsymbol{r},q)=\min\{t\mid\boldsymbol{r}^{\prime}(t)\in\mathcal{T},\boldsymbol{r}^{\prime}(0)=\boldsymbol{r},q^{\prime}(0)=q\}$ is a random variable whose expectation we are trying to minimize¹¹1 Strictly speaking, the objective in [1] is different since they minimize $\mathbb{E}\left[\int_{0}^{T}e^{-\lambda t}dt\right]=\mathbb{E}\left[\left(1-\exp(-\lambda T)\right)/\lambda\right],$ where $\lambda>0$ is a time-discounting coefficient chosen to improve the convergence of value iterations. Our objective is obtained by taking the limit with $\lambda\to 0$ and the numerical tests indicate that value iterations based on our discretization converge even in this “undiscounted” case. by choosing a feedback steering/tack-switching policy $\mu$. This $T_{\mu}(\boldsymbol{r},q)$ can be split into two pieces: the time spent continuously controlling the sails plus the time spent performing tack-switches. Each tack-switch incurs a switching delay $C$ due to the speed loss from the wind pushing directly against the boat’s intended direction during the maneuver²²2For simplicity of modeling, we will treat this delay as constant, though state-dependent $C$ would not present additional computational challenges. . We model the boat’s propelled speed as falling to zero (i.e., $s(u)=0$) during a tack-switch.

We thus define the value function to encode the minimal expected time to $\mathcal{T}$ form each starting configuration:

By Bellman’s Optimality Principle, this $v$ must satisfy

where the first argument to the $\min$ is to be interpreted as the best expected time-to-$\mathcal{T}$ if we stay on the current tack at least for a small time $\tau$, while the second is the best expected time if we switch tacks immediately. More formally, we define the operators $\Sigma_{\mu,\tau}$ and $\mathcal{N}$ as:

where we employ the compact notation $\mathbb{E}_{\boldsymbol{r}^{\prime};\mu,\tau}[\cdot|\boldsymbol{r}]$ to indicate the expectation of a function over evolved positions $\boldsymbol{r}^{\prime}$, starting from position $\boldsymbol{r}$ and using control policy $\mu$ for a time $\tau$. In (7), the symbol $\O$ in place of the control policy indicates that during the tacking time $C$, no steering control is taking place as $s(u)=0$ for the duration. Importantly, the version of operator (7) used in [1] does not treat $\theta$ as evolving during the tack, instead setting $\mathcal{N}v(\boldsymbol{r},q)=v(\boldsymbol{r},3-q)$. As we show in section V, correctly incorporating these switch-duration dynamics is critical in high-drift problems.

A standard argument based on an Ito-Taylor expansion of (5) yields the quasi-variational inequality that must be satisfied by a smooth value function:

with the Hamiltonian

If $v$ is not smooth, one can still interpret it as a weak solution (in particular, the unique viscosity solution [6]) of (8).

As $\theta$ is the only randomly evolving variable and obeys drift-diffusive dynamics, both expectations in (6) and (7) take a similar form as averaging $v$ over Gaussian-distributed final $\theta^{\prime}$. For small $\tau$, an approximation to (6) taking $u=\mu(\boldsymbol{r},q)$ and $\dot{r},\dot{\theta}$ constant for a time $\tau$ yields:

where $v$ is interpreted as $2\pi$-periodic in $\theta^{\prime}$, and $\dot{r},\dot{\theta}$ both implicitly depend on $(r,\theta,q,u)$. Eq. (7) is the same formula under the substitutions $q\to 3-q,\tau\to C,\dot{r}=0,\dot{\theta}=a$.

We now apply our algorithm to solve the optimal control problem (3)-(4) for a variety of system parameters, and sample stochastic optimal trajectories for each scenario³³3 To ensure the computational reproducibility, we provide the full source code, movies illustrating these simulations, and additional convergence data at https://eikonal-equation.github.io/Stochastic_Sailing/ . For all results shown, the target radius is $R=0.1$ and the domain is $(r,\theta)\in[R,2.0]\times[0,2\pi]$. Target boundary conditions are set as $v(r,\theta,q)=0,\;\forall r<R$ and exit boundary conditions are $v(r,\theta,q)=10^{6},\;\forall r>2.0.$

We first focus on problems with zero drift $(a=0)$, but varying the diffusivity coefficient $\sigma$. In Fig. 3, we plot switchgrid solutions for a collection of $\sigma$ values, noting that these grids are defined and plotted in $(r,\theta)$ state space, not the absolute position space. Points are marked red if $\texttt{switchgrid}(r,\theta,q=1)=1$, blue if $\texttt{switchgrid}(r,\theta,q=2)=1$, and white if both are $0$. The darker shaded regions correspond to the switchgrids obtained if we use the zero-evolution switching operator of [1] as compared to our (7). The black circle at the center marks the target set $\mathcal{T}$. In general, stochastic optimal trajectories “bounce” between these red and blue switching fronts as they approach $\mathcal{T}$.

With $\sigma=0,$ the switching fronts lie at constant $\theta$, corresponding to the angles that locally maximize the projection of the polar speed plot (Fig. 1(a)) onto the axis of the wind. Once a sailboat hits this front, it can switch tack and stay at this constant angle that maximizes its speed towards the target (see Fig. 4(a)). As $\sigma$ increases, these switching fronts contract more tightly around $\theta=0^{\circ},180^{\circ}$. As the future wind state grows more uncertain with larger $\sigma$, the optimal policy becomes more conservative and constrains the boat to smaller excursions in $\theta$ before switching to the opposite tack. We observe that in these zero-drift problems, our improved switching operator makes the switching grids modestly more conservative, switching earlier by a few degrees in Fig. 3(b,c).

In Fig. 4, we plot samples of stochastic optimal trajectories corresponding to these control policies, with blue markers indicating where tack switches occur. Similarly to [1], we initialize six boats at positions $(r,\theta)\in\{(1.8,0.0),(1.93,\pm 0.39)\},$ and tacks $q\in\{1,2\}$. We then evolve and control each boat’s state until it reaches its target as described in Sec. IV-C. We can see that, as foreshadowed by the switchgrids, at higher values of $\sigma$ the trajectories tend to stay within tighter cones from the target.

Next, in Figs. 5 and 6 we show switchgrids and sampled trajectories for problems in which the drift parameter $a$ is nonzero, with a fixed $\sigma=0.05$. As the drift increases, we observe the switchgrids and the trajectories becoming skewed and asymmetric towards negative $\theta$. This reflects that it is optimal to “lead” the wind in anticipation of future wind evolution, allowing the boat to stay in the portion of state space that allows it to steer at the optimal $~{}40^{\circ}\mathchar 45\relax 45^{\circ}$ angle while the wind evolves.

Additionally, some “wave”-like features start to emerge in the switchgrids for large $a$. Empirically, we observe that the number of “waves” in the switchgrid corresponds to the number of $180^{\circ}$ turns the wind tends to make by the time an average boat at $r=R$ reaches the target. Each “wave” can thus be interpreted as corresponding to the optimal control during one (average) half-cycle of the wind. We notice that accurately modeling the switching operator is critical in these high-drift scenarios, as the two options produce dramatically different wave-like structures in the switchgrids.

We have introduced an accelerated semi-Lagrangian approach to solve a hybrid stochastic optimal control problem arising in match race sailing. By working in a reduced coordinate system and designing a nearly-causal discretization, we were able to produce control policies on high-resolution grids which provide insight into tack-switching strategies under unpredictable weather.

While the conversion to reduced coordinates greatly accelerates our solution, it also introduces several limitations since we had to assume that the problem is radially symmetric. This restricts the target shape to a circle, and does not allow the inclusion of state constraints such as obstacles or nearby coastline. However, it does still allow for modeling match races between multiple agents. Our approach will require keeping track of $(r,\theta)$ coordinates for each agent, yielding an overall state space with one less dimension than in [4]. For this reason, we hope that our approach will allow to rapidly test the implications of various weather models or match strategies in “ideal” conditions far away from any barriers.

This paper is focused on a risk-neutral task of finding a feedback policy $\mu$ that minimizes the expected time to target $\mathbb{E}\left[T_{\mu}\right]$. However, practitioners often prefer robust control policies more suitable when bad outcomes are costly. Two such popular approaches are the “risk-sensitive” control [12, 13] and $H^{\infty}$ control [14]. But neither of these provides any direct bounds on the likelihood of bad scenarios; i.e., $w_{\mu}(c)=\mathbb{P}\left(T_{\mu}\geq c\right).$ An efficient method for minimizing $w_{\mu}(c)$ for all threshold values $c$ simultaneously was recently introduced for piecewise-deterministic Markov processes in [15]. We hope that it can be similarly extended to the current context. It would be also interesting to extend the methods of multiobjective optimal control (e.g., [16, 17]) to find Pareto-optimal policies reflecting the rational tradeoffs between conflicting optimization criteria (e.g., $\mathbb{E}\left[T_{\mu}\right]$ vs $w_{\mu}(c)$ or even $w_{\mu}(c_{1})$ vs $w_{\mu}(c_{2})$).