Illustration : Crossing a river#
This example deal with the case of a person who wants to cross a river in a minimum of time.
Problem statement:
the velocity on the ground is 10 km/h
the swimming velocity is 1 km/h
the starting point and the ending point are fixed
the width of the river is fixed (40m)
Based on this problem we will go through the formulation of the optimization method and its resolution. The following points have to be addressed :
Problem formulation
Objective
Variables
Constraints
Objective function
Resolution method
Problem statement : model formulation#
Where is the river : model the river#
The river is an area wher the velocity is lower.
groung_v = 10.0 / 3.6
swim_v = 1.0 / 3.6
y1_river = 40
y2_river = 80
def world_definition(x, y):
v = swim_v * np.ones_like(x)
v[y <= y1_river] = groung_v
v[y > y2_river] = groung_v
return v
x = np.linspace(0, 100, 100)
y = np.linspace(0, 120, 150)
x_world, y_world = np.meshgrid(x, y, indexing="xy")
v_world = world_definition(x_world, y_world)
plt.figure()
plt.contourf(x_world, y_world, v_world, 20, cmap=cm.terrain, vmax=10.0)
plt.title("Velocity map [m/s]")
plt.colorbar()
plt.show()
Defintion of the path#
path = np.array([[5.0, 5.0], [15.0, y1_river], [25, y2_river], [80, 85]])
plt.figure()
plt.contourf(x_world, y_world, v_world, 20, cmap=cm.terrain, vmax=10.0)
plt.plot(path[:, 0], path[:, 1], "go-")
plt.title("Velocity map [m/s]")
plt.colorbar()
<matplotlib.colorbar.Colorbar at 0xffff3a979f10>
Model the river crossing : compute the velocity and the duration to go trougth the path#
first lrt’s work on one segment of the path :
def arc_duration(seg, nb_sub=20):
seg = seg.T
sub = np.linspace(seg[:, 0], seg[:, 1], nb_sub)
sub = sub[:-1] + np.diff(sub, axis=0) / 2
sub_length = np.linalg.norm(np.diff(seg)) / (nb_sub - 1)
v = world_definition(sub[:, 0], sub[:, 1])
duration = (sub_length / v).sum()
return duration, sub, v
# test for one segment of the path
seg = path[:2, :]
duration, sub, v = arc_duration(seg, nb_sub=20)
duration
np.float64(13.104197800704934)
Now we can work on the full path :
total = 0
for i in range(len(path) - 1):
seg = path[i : i + 2, :]
duration, sub, v = arc_duration(seg, nb_sub=20)
total += duration
print("segment {0} : {1:0.2f}s".format(i, duration))
print("Total : {0:0.2f}s".format(total))
segment 0 : 13.10s
segment 1 : 148.43s
segment 2 : 19.88s
Total : 181.42s
def total_duration(path, blabla=False):
total = 0.0
for i in range(len(path) - 1):
seg = path[i : i + 2, :].copy()
duration, _, _ = arc_duration(seg, nb_sub=20)
total += duration
if blabla:
print("segment {0} : {1:0.1f}s".format(i, duration))
if blabla:
print("Total : {0:0.1f}s".format(total))
return total
# test for the path
total_duration(path, blabla=True)
segment 0 : 13.1s
segment 1 : 148.4s
segment 2 : 19.9s
Total : 181.4s
np.float64(181.41765015387776)
Widget plot#
plt.figure()
(l,) = plt.plot(path[:, 0], path[:, 1], "y-o")
plt.contourf(x_world, y_world, v_world, 20, cmap=cm.terrain, vmax=10.0)
title = plt.title("duration")
@ipw.interact(x1=(0.0, 100, 1), x2=(0.0, 100, 1))
def update(x1=30, x2=30):
path[1, 0] = x1
path[2, 0] = x2
t_total = total_duration(path, blabla=False)
l.set_data(path[:, 0], path[:, 1])
title.set_text("Total : {0:0.1f}s".format(t_total))
What is the optimal solution ?#
First idea : test all combination#
Nx, Ny = 30, 30
x1 = np.linspace(5, 100, Nx)
x2 = np.linspace(5, 100, Ny)
X, Y = np.meshgrid(x1, x2)
Z = np.zeros_like(X)
for i in range(len(X)):
for j in range(len(Y)):
path[1, 0] = X[i, j]
path[2, 0] = Y[i, j]
Z[i, j] = total_duration(path, blabla=False)
Find the min#
x1_min, x2_min = np.where(np.min(Z) == Z)
X_min = X[x1_min, x2_min]
Y_min = Y[x1_min, x2_min]
print(X_min, Y_min)
[67.24137931] [70.51724138]
plt.figure()
title = plt.title("")
# plt.contourf(X, Y, Z, 20, cmap=cm.jet)
# plt.colorbar()
# plt.contour(X, Y, Z, 20, cmap=cm.gray)
plt.scatter(X, Y, c=Z, s=50, cmap=cm.jet)
plt.plot(X_min, Y_min, "*y")
plt.colorbar()
plt.xlabel("Diving position [m]")
plt.ylabel("Dout of water [m]")
plt.show()
Is that realistic with a finer grid ? (or in larger dimension)#
%%timeit
total_duration(path, blabla=False)
50.4 μs ± 1.5 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
Usning optimization methode#
def cost(X):
path[1, 0] = X[0]
path[2, 0] = X[1]
return total_duration(path, blabla=False)
from scipy import optimize
sol = optimize.minimize(cost, [10.0, 10.0], method="Nelder-Mead")
plt.figure()
plt.contourf(x_world, y_world, v_world, 20, cmap=cm.terrain, vmax=10.0)
plt.plot(path[:, 0], path[:, 1], "go-", label="optimal path")
plt.title("Velocity map [m/s]")
plt.legend()
plt.colorbar()
<matplotlib.colorbar.Colorbar at 0xffff3a3cf590>
Crossing another river : higer dimension#
The person can swim faster close to the river side.
def world_definition(x, y):
v = swim_v * np.ones_like(x)
v = swim_v + (1 - 0.8 * np.cos(2 * np.pi * (y - 60) / 40))
v[y <= y1_river] = groung_v
v[y > y2_river] = groung_v
return v
x = np.linspace(0, 100, 100)
y = np.linspace(0, 120, 150)
x_world, y_world = np.meshgrid(x, y, indexing="xy")
v_world = world_definition(x_world, y_world)
plt.figure()
plt.contourf(x_world, y_world, v_world, 20, cmap=cm.terrain, vmax=10.0)
plt.title("Velocity map")
plt.colorbar()
plt.show()
Is that still relevant to go straight ?
Add more control points#
path_init = np.linspace([5.0, 5.0], [80.0, 100.0], 10)
path = np.copy(path_init)
def cost_n_points(X):
for i in range(len(X)):
path[i + 1, 0] = X[i]
return total_duration(path, blabla=False)
init_guess = np.copy(path_init[:-2, 0])
sol = optimize.minimize(cost_n_points, init_guess, method="Nelder-Mead")
X = sol.x
for i in range(len(X)):
path[i + 1, 0] = X[i]
sol
message: Optimization terminated successfully.
success: True
status: 0
fun: 71.70370204192615
x: [ 1.684e+01 2.868e+01 4.051e+01 4.781e+01 5.043e+01
5.217e+01 5.719e+01 6.816e+01]
nit: 899
nfev: 1330
final_simplex: (array([[ 1.684e+01, 2.868e+01, ..., 5.719e+01,
6.816e+01],
[ 1.684e+01, 2.868e+01, ..., 5.719e+01,
6.816e+01],
...,
[ 1.684e+01, 2.868e+01, ..., 5.719e+01,
6.816e+01],
[ 1.684e+01, 2.868e+01, ..., 5.719e+01,
6.816e+01]], shape=(9, 8)), array([ 7.170e+01, 7.170e+01, 7.170e+01, 7.170e+01,
7.170e+01, 7.170e+01, 7.170e+01, 7.170e+01,
7.170e+01]))
plt.figure()
plt.contourf(x_world, y_world, v_world, 20, cmap=cm.terrain, vmax=10.0, alpha=0.5)
plt.plot(path_init[:, 0], path_init[:, 1], "o:g")
plt.plot(path[:, 0], path[:, 1], "o-r")
plt.title("Velocity map")
plt.colorbar()
<matplotlib.colorbar.Colorbar at 0xffff3a2ee360>
Add constrainte : Crossing a river with a crocodile#
# Low velocity
def invexp(x, y, tau, cx, cy):
return np.exp(-(((x - cx) / tau) ** 2) - ((y - cy) / tau) ** 2)
def world_definition(x, y):
v = swim_v * np.ones_like(x)
v = swim_v + (1 - 0.8 * np.cos(2 * np.pi * (y - 60) / 40))
tau = 5
v = v * ((1 + 0.1) - invexp(x, y, tau, 50, 60))
v[y <= y1_river] = groung_v
v[y > y2_river] = groung_v
return v
x = np.linspace(0, 100, 100)
y = np.linspace(0, 120, 150)
x_world, y_world = np.meshgrid(x, y, indexing="xy")
v_world = world_definition(x_world, y_world)
plt.figure()
plt.contourf(x_world, y_world, v_world, 20, cmap=cm.terrain, vmax=10.0)
plt.title("Velocity map")
plt.colorbar()
plt.show()
With 2 parametres (2D Problems)#
path = np.array([[5.0, 5.0], [15.0, y1_river], [25, y2_river], [80, 85]])
Nx, Ny = 25, 25
x1 = np.linspace(5, 100, Nx)
x2 = np.linspace(5, 100, Ny)
X, Y = np.meshgrid(x1, x2)
Z = np.zeros_like(X)
for i in range(len(X)):
for j in range(len(Y)):
path[1, 0] = X[i, j]
path[2, 0] = Y[i, j]
Z[i, j] = total_duration(path, blabla=False)
x1_min, x2_min = np.where(np.min(Z) == Z)
X_min = X[x1_min, x2_min]
Y_min = Y[x1_min, x2_min]
print(X_min, Y_min)
[56.45833333] [72.29166667]
plt.figure()
title = plt.title("")
plt.contourf(X, Y, Z, 20, cmap=cm.jet)
plt.colorbar()
plt.contour(X, Y, Z, 20, cmap=cm.gray)
plt.scatter(X, Y, c=Z, s=20, cmap=cm.jet)
plt.plot(X_min, Y_min, "or")
plt.xlabel("Diving position [m]")
plt.ylabel("Dout of water [m]")
plt.show()
With N parametres (ND Problems)#
N = 20
path_init = np.linspace([5.0, 5.0], [80.0, 100.0], N)
path = np.copy(path_init)
def cost_n_points(X):
for i in range(len(X)):
path[i + 1, 0] = X[i]
return total_duration(path, blabla=False)
init_guess
array([ 5. , 13.33333333, 21.66666667, 30. , 38.33333333,
46.66666667, 55. , 63.33333333])
init_guess = 5 + np.random.rand(path_init[:-2, 0].shape[0]) * 80.0
path_init[1:-1, 0] = init_guess
sol = optimize.minimize(cost_n_points, init_guess, method="Nelder-Mead")
X = sol.x
for i in range(len(X)):
path[i + 1, 0] = X[i]
plt.figure()
plt.contourf(x_world, y_world, v_world, 20, cmap=cm.terrain, vmax=10.0, alpha=0.5)
plt.plot(path_init[:, 0], path_init[:, 1], "o:g")
plt.plot(path[:, 0], path[:, 1], "o-r")
plt.title("Velocity map")
plt.colorbar()
<matplotlib.colorbar.Colorbar at 0xffff38947800>