In this notebook, we compute the route that minimizes fuel consumption for a simplified boat model, with time-independent data, and without constraint on the arrival time. The instantaneous fuel consumption of the boat, moving a velocity $v$ at a position $x$, is a quadratic function $$ c(x,v) = \mu(x) + \frac 1 2 \| v - \omega(x) \|_{M(x)}^2. $$ The parameters have the following interpretation:
Discussion of the assumptions
Our assumptions are obviously an strong simplification of reality. Some of them may be relaxed, at the price of more complexity in the implementation and in the modeling phase.
Any cost which is convex and with at least quadratic growth at infinity, can be approximated in this form. However, determining this cost (which depends on the physical properties of the boat), and approximating it (fitting a convex set with an intersection of ellipsoids), may raise difficulties in itself.
Summary of volume Fast Marching Methods, this series of notebooks.
Main summary of the Adaptive Grid Discretizations book of notebooks, including the other volumes.
This Python® notebook is intended as documentation and testing for the HamiltonFastMarching (HFM) library, which also has interfaces to the Matlab® and Mathematica® languages. More information on the HFM library in the manuscript:
Copyright Jean-Marie Mirebeau, University Paris-Sud, CNRS, University Paris-Saclay
import sys; sys.path.insert(0,"..") # Allow import of agd from parent directory (useless if conda package installed)
#from Miscellaneous import TocTools; print(TocTools.displayTOC('BoatRouting','FMM'))
from agd import Eikonal
from agd import LinearParallel as lp
from agd import FiniteDifferences as fd
from agd.Metrics import Rander,Riemann
from agd import AutomaticDifferentiation as ad
from agd.Plotting import savefig,quiver; #savefig.dirName = 'Images/BoatRouting'
import numpy as np; xp = np
from copy import copy
import matplotlib.pyplot as plt
Utility functions
norm_infinity = ad.Optimization.norm_infinity
from agd.ExportedCode.Notebooks_NonDiv.LinearMonotoneSchemes2D import streamplot_ij
Uncomment the following line to use the GPU eikonal solver.
#xp,plt,quiver,Eikonal = map(ad.cupy_friendly,(xp,plt,quiver,Eikonal))
We assume in this section that the cost is independent of the position $c=c(v)$. In other words, the fuel consumption at rest $\mu>0$, the water currents $\omega$, and the problem geometry $M\succ 0$ are independent of time.
Assume that the boat needs to move from the origin to a prescribed position $z$, without time constraint. The optimal path is to go straight toward the objective, and the minimal cost is thus $$ N(z) := \inf_{t>0} \ t \ c(z/t). $$ For the quadratic cost model, this problem reads $$ N(z) = \min_{t>0} \ t \ \Big(\mu + \frac 1 2 \big\|\frac z t -\omega\big\|_M^2\Big). $$ It is easily solvable (exercise), and the optimal time is $$ t_opt = \frac {\|z\|_M}{\sqrt{2 \mu + \|\omega\|_M^2}}, $$ whereas the minimal cost reads $$ N(z) = \|z\|_M \sqrt{2 \mu + \|\omega\|_M^2} - <\omega,z>_M. $$
def route_min(z,params):
z,μ,ω,M = fd.common_field((z,)+params, depths=(1,0,1,2))
z_norm = np.sqrt(lp.dot_VAV(z,M,z))
μω_norm = np.sqrt( 2*μ +lp.dot_VAV(ω,M,ω) )
cost = z_norm*μω_norm - lp.dot_VAV(z,M,ω)
time = z_norm / μω_norm
fuel = cost/time
rvel = z/time - ω
return {
'cost':cost, # minimal cost for this travel
'time':time, # optimal travel time
'fuel':fuel, # instantaneous fuel consumption
'rvel':rvel, # relative velocity, w.r.t current
}
The cost $N(z)$ takes the form of a Rander norm, a.k.a the sum of a distorted euclidean norm and of a linear function. Note that the Rander compatibility condition is satisfied: $N(z)\geq 0$ for all $z$. The more general cost considered in the introduction would yields a supremum of functions of this form.
def metric(params):
μ,ω,M = fd.common_field(params,depths=(0,1,2))
return Rander( M*(2*μ + lp.dot_VAV(ω,M,ω)), -lp.dot_AV(M,ω))
# Generate a cartesian grid
aX = xp.linspace(-1,1)
X = ad.array(np.meshgrid(aX,aX,indexing='ij'))
The following observations can be made.
# Parameters : unit weight for time, unit horizontal drift, euclidean geometry
params = (1.,xp.array((1.,0.)),xp.eye(2))
route = route_min(X,params)
assert np.allclose(metric(params).norm(X),route['cost'])
plt.figure(figsize=[12,10])
plt.subplot(2,2,1); plt.axis('equal')
plt.title('Minimal travel cost')
plt.contour(*X,route['cost']); plt.colorbar()
plt.subplot(2,2,2); plt.axis('equal')
plt.title('Travel time for the optimal path')
plt.contour(*X,route['time']); plt.colorbar()
plt.subplot(2,2,3); plt.axis('equal')
plt.title('Instantaneous fuel consumption')
plt.contour(*X,route['fuel']); plt.colorbar()
plt.subplot(2,2,4); plt.axis('equal')
plt.title('Velocity relative to water currents')
quiver(*X,*route['rvel'],subsampling=(5,5))
<matplotlib.quiver.Quiver at 0x12c715f40>
By setting $\mu=0$, we do not penalize travel time. In that case:
# Parameters : unit weight for time, unit horizontal drift, euclidean geometry
params = (0.,xp.array((1.,0.)),xp.eye(2))
route = route_min(X,params)
plt.figure(figsize=[12,10])
plt.subplot(2,2,1); plt.axis('equal')
plt.title('Minimal travel cost')
plt.contour(*X,route['cost']); plt.colorbar()
plt.subplot(2,2,2); plt.axis('equal')
plt.title('Travel time for the optimal path')
plt.contour(*X,route['time']); plt.colorbar()
plt.subplot(2,2,3); plt.axis('equal')
plt.title('Instantaneous fuel consumption')
plt.contour(*X,route['fuel']); plt.colorbar()
plt.subplot(2,2,4); plt.axis('equal')
plt.title('Velocity relative to water currents')
quiver(*X,*route['rvel'],subsampling=(5,5))
<matplotlib.quiver.Quiver at 0x12d3da750>
By setting $\mu$ to a large value, we strongly penalize travel time. In that case:
# Parameters : unit weight for time, unit horizontal drift, euclidean geometry
params = (1e3,xp.array((1.,0.)),xp.eye(2))
route = route_min(X,params)
plt.figure(figsize=[12,10])
plt.subplot(2,2,1); plt.axis('equal')
plt.title('Minimal travel cost')
plt.contour(*X,route['cost']); plt.colorbar()
plt.subplot(2,2,2); plt.axis('equal')
plt.title('Travel time for the optimal path')
plt.contour(*X,route['time']); plt.colorbar()
plt.subplot(2,2,3); plt.axis('equal')
plt.title('Instantaneous fuel consumption')
plt.contour(*X,route['fuel']); plt.colorbar()
plt.subplot(2,2,4); plt.axis('equal')
plt.title('Velocity relative to water currents')
quiver(*X,*route['rvel'],subsampling=(5,5))
<matplotlib.quiver.Quiver at 0x12d050800>
In this section, we considere a setup with space dependent geometry, water currents, and obstacles. Finding the optimal route requires solving an eikonal equation, which we do using the HFM library.
The first step is to define the location of the obstacles, the starting point, and the target points.
#ReloadPackages()
hfmIn = Eikonal.dictIn({
'model':'Rander2', # Riemannian + drift, what is needed here
'exportValues':1,
'exportGeodesicFlow':1,
'seed':[-1.7,0.6], # Where to start the front propagation
})
hfmIn.SetRect([[-2,2],[-1,1]],dimx=200) # Rectangular domain
X = hfmIn.Grid() # Coordinate system
hfmIn['walls'] = X[1]-np.abs(X[0])>=0 # Obstacles in the domain
hfmIn.SetUniformTips((6,3)) # Points from which to backtrack geodesics
hfmIn['dims']
array([200., 100.])
plt.title('Obstacle inserted in the domain'); plt.axis('equal')
plt.contourf(*X, hfmIn['walls']);
We now define the Riemannian metric $M$. Here, we choose the metric attached to a sphere, a.k.a. the earth, of radius one. Recall that the intrinsic metric on a manifold embedded in Euclidean space is $$ M(x) = df(x)^T df(x) = \nabla f(x) \nabla f(x)^T, $$ where $df$ is the differential of $f$, and $\nabla f(x) := df(x)^T$.
def Spherical(θ,ϕ):
"""Spherical embedding: θ is longitude, ϕ is latitude from equator toward pole"""
return (np.cos(θ)*np.cos(ϕ), np.sin(θ)*np.cos(ϕ), np.sin(ϕ))
def IntrinsicMetric(Embedding,*X):
"""Riemannian metric for a manifold embedded in Euclidean space"""
X_ad = ad.Dense.identity(constant=X,shape_free=(2,)) # First order dense AD variable
Embed_ad = ad.asarray(Embedding(*X_ad)) # Differentiate the embedding
Embed_grad = Embed_ad.gradient()
Embed_M = lp.dot_AA(Embed_grad,lp.transpose(Embed_grad)) # Riemannian metric
return Embed_M
We now choose the water (surface) currents, arbitrarily. They are defined in the local chart $\theta,\phi$, and not on the three dimensional sphere, for simplicity.
def bump(x,y):
"""Gaussian-like bump (not normalized)"""
return np.exp(-(x**2+y**2)/2)
def Currents(θ,ϕ):
"""Some arbitrary vector field (water currents)"""
bump0 = bump(θ+1,(ϕ+0.3)*2); ω0=(0,1) # intensity and direction of the currents
bump1 = 2*bump(2*(θ-0.7),ϕ-0.2); ω1=(1,-1)
bump0,ω0,bump1,ω1 = fd.common_field( (bump0,ω0,bump1,ω1), depths=(0,1,0,1))
return bump0*ω0+bump1*ω1
We have put two strong currents:
Embed_ω = Currents(*X)
Embed_M = IntrinsicMetric(Spherical,*X) # Actually a diagonal matrix, but we don't exploit this fact here
plt.title('Drift vector field')
quiver(*X,*Embed_ω,subsampling=(10,10))
<matplotlib.quiver.Quiver at 0x12d578800>
params = (1.,Embed_ω,Embed_M) # Idle fuel consumption μ = 1
hfmIn['metric'] = metric(params)
It is worth noting that the problem metric is quite strongly anisotropic, due to the presence of strong currents. Recall that the anisotropy, at a given point, is the maximum ratio of the travel cost from that point in two different directions of unit Euclidean norm. A.k.a it is much cheaper to be pushed by the current than to go against it.
# only an upper bound. (Slow on GPU, and may cause weird bugs)
if xp is np: np.max(hfmIn['metric'].anisotropy_bound())
The problem is already anisotropic (but not asymmetric) without the currents, due to the geometry of the sphere. However, that anisotropy is much less pronounced.
#np.max(Riemann(Embed_M).anisotropy())
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field cosAngleMin defaults to 0.5 Field refineStencilAtWallBoundary defaults to 0 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.009177 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45
Some observations :
plt.figure(figsize=[16,4])
plt.subplot(1,2,1); plt.axis('equal')
plt.title('Minimal travel cost')
plt.contour(*X,hfmOut['values'],levels=50)
plt.subplot(1,2,2); plt.axis('equal')
plt.title('Optimal paths, and inserted obstacle')
plt.contourf(*X,hfmIn['walls'],cmap='Greys')
for geo in hfmOut['geodesics']: plt.plot(*geo)
The HFM library exports the geodesic flow, which is the local direction of the minimal geodesics coming from the seed point. Using this quantity we can recover the instantaneous fuel consumption, and the relative velocity of the boat. One observes that:
v = hfmOut['flow']
v[:,ad.Optimization.norm(v,axis=0)==0] = np.nan # Avoid division by zero where there is no velocity (walls, seed)
route = route_min(v,params)
plt.figure(figsize=[16,4])
plt.subplot(1,2,1); plt.axis('equal')
plt.title('Instantaneous fuel consumption')
plt.contour(*X,route['fuel'],levels=50)
plt.subplot(1,2,2); plt.axis('equal')
plt.title('Relative velocity of the boat w.r.t water')
quiver(*X,*route['rvel'],subsampling=(10,10))
<matplotlib.quiver.Quiver at 0x12d8c75f0>
One may also recover the total fuel consumption, and the total travel time, integrated along the geodesics, by solving one additional first order PDE.
We can increase the cost of time. As before, if we increase it a lot, then the drift due to water currents becomes negligible.
params = (1e3,Embed_ω,Embed_M) # Strong idle fuel consumption μ
hfmIn['metric'] = metric(params)
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field cosAngleMin defaults to 0.5 Field refineStencilAtWallBoundary defaults to 0 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.006859 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45
In the limit, as $\alpha \to \infty$, the geodesics coincide with the great circles on the sphere - except when they hit the obstacles or the domain boundary.
plt.figure(figsize=[16,4])
plt.subplot(1,2,1); plt.axis('equal')
plt.title('Minimal travel cost')
plt.contour(*X,hfmOut['values'],levels=50)
plt.subplot(1,2,2); plt.axis('equal')
plt.title('Optimal paths, and inserted obstacle')
plt.contourf(*X,hfmIn['walls'],cmap='Greys')
for geo in hfmOut['geodesics']: plt.plot(*geo)
The fuel consumption varies much less, and the boat heads directly where it needs to go.
v = hfmOut['flow']
v[:,ad.Optimization.norm(v,axis=0)==0] = np.nan # Avoid division by zero where there is no velocity (walls, seed)
route = route_min(v,params)
plt.figure(figsize=[16,4])
plt.subplot(1,2,1); plt.axis('equal')
plt.title('Instantaneous fuel consumption')
plt.contour(*X,route['fuel'],levels=50); plt.colorbar()
plt.subplot(1,2,2); plt.axis('equal')
plt.title('Relative velocity of the boat w.r.t water')
quiver(*X,*route['rvel'],subsampling=(10,10))
<matplotlib.quiver.Quiver at 0x12dc40800>
In one removes the drift (equivalently $\alpha \to \infty$), the instantaneous fuel consumption becomes constant along the journey, which is expected.
params = (1.,(0,0),Embed_M) # Weight α = 1 for time w.r.t fuel
hfmIn['metric'] = metric(params)
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field cosAngleMin defaults to 0.5 Field refineStencilAtWallBoundary defaults to 0 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.006822 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45
v = hfmOut['flow']
v[:,ad.Optimization.norm(v,axis=0)==0] = np.nan # Avoid division by zero where there is no velocity (walls, seed)
route = route_min(v,params)
print(f"Extremes of fuel instantaneous consumption. min={np.nanmin(route['fuel'])}, max={np.nanmax(route['fuel'])}")
Extremes of fuel instantaneous consumption. min=2.0, max=2.0000000000000004
In order to obtain the arrival time at some destination, we must integrate it along the geodesic. We take advantage of the automatic differentiation capabilities of the HFM library, in a somewhat hacky way, to perform this integration and obtain the arrival times.
def ArrivalTime(hfmIn,params):
hfmIn = copy(hfmIn)
# if hfmIn.xp is not np: hfmIn['solver']='AGSI' #TODO : why needed ?
hfmIn['metric'] = metric(params)
hfmIn['exportGeodesicFlow']=1
cache = Eikonal.Cache(needsflow=True)
hfmOut = hfmIn.Run(cache=cache)
flow = hfmOut['flow']
no_flow = np.all(flow==0,axis=0)
flow[:,no_flow]=np.nan # No flow at the seed point, avoid zero divide
route = route_min(flow,params)
costVariation = route['time']
costVariation[no_flow] = 0
hfmIn['costVariation'] = np.expand_dims(costVariation,axis=-1)
hfmOut2 = hfmIn.Run(cache=cache) # cache avoids some recomputations
time = hfmOut2['values'].gradient(0)
return time,hfmOut
params = (1.,Embed_ω,Embed_M) # Idle fuel consumption μ = 1
arrival_time,hfmOut = ArrivalTime(hfmIn,params)
Requesting cacheable data Field verbosity defaults to 1 Field cosAngleMin defaults to 0.5 Field refineStencilAtWallBoundary defaults to 0 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.009161 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45 Filling cache data Providing cached data Field verbosity defaults to 1 Field cosAngleMin defaults to 0.5 Field refineStencilAtWallBoundary defaults to 0 Field order defaults to 1 Field seedRadius defaults to 0 Bypassing fast marching solver based on cached data. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45
plt.title('Arrival time for the optimal path'); plt.axis('equal')
plt.contourf(*X,arrival_time);
We compare the exact analytical solution and the numerical solution in the case of a constant medium, for validation.
hfmIn = Eikonal.dictIn({
'verbosity':0,
'model':'Rander2',
'exportValues':1,
'seed':[0,0],
# 'solver':'AGSI' #TODO : why needed ?
})
hfmIn.SetRect([[-1,1],[-1,1]],dimx=101,sampleBoundary=True)
params = (1.,xp.array((0.,1.)),xp.eye(2))
arrival_time,hfmOut = ArrivalTime(hfmIn,params)
plt.contourf(*hfmIn.Grid(),hfmOut['values']);
X = hfmIn.Grid()
route = route_min(X,params)
assert norm_infinity(arrival_time - route['time']) < 0.02
assert norm_infinity(hfmOut['values'] - route['cost']) < 0.05