Show code cell source
%matplotlib notebook
#%matplotlib widget
import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize
import ipywidgets as ipw
from matplotlib import cm
import ipywidgets as widgets
from scipy import integrate
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 0x7f197580b340>
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
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
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)
131 µs ± 2.92 µ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 0x7f197559d930>
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]]), 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 0x7f19754c0f40>
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 = 10
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 0x7f197547c820>