pm21-dragon/lectures/lecture-07/2 - optimization.ipynb
2024-11-29 09:11:00 +01:00

6.7 KiB

None <html> <head> </head>

Optimization

Let's consider the following code:

In [ ]:
import numpy as np
import matplotlib.pyplot as plt

# position of flower in 2D
flower = np.array([29.1, 19.1])

def make_bee_track(t):
    """Find bee position at time t"""
    pos0 = (13.21, 4.56)
    velocity = (3.1, 0.25)
    pos_x = pos0[0] + t*velocity[0]
    pos_y = pos0[1] + t*velocity[1]
    return np.array([pos_x,pos_y])

t = np.linspace(0,15,10)
bee_track = make_bee_track(t)

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4,3))
ax.plot( [flower[0]], [flower[1]], 'or', label='flower' )
ax.plot( bee_track[0], bee_track[1], '.-k', label='bee')
ax.axis('equal')
ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('y')

In the above code, we parameterized the bee trajectory by the variable t in the function make_bee_track(). This means we could get a new point on the track by choosing a new value of t. For example:

In [ ]:
print(make_bee_track(0.0))
print(make_bee_track(0.1))
print(make_bee_track(0.2))
print(make_bee_track(1.0))

Now let's measure the distance between the bee and the flower.

In [ ]:
def compute_distance(a,b):
    a = np.array(a)
    b = np.array(b)
    return np.sqrt(np.sum((a-b)**2))

n_time_steps = bee_track.shape[1]
distance = np.zeros(n_time_steps)
for i in range(n_time_steps):
    bee_pos = bee_track[:,i]
    distance[i] = compute_distance(bee_pos, flower)
display(distance)

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4,3))
ax.plot( t, distance )
ax.set_xlabel('t')
ax.set_ylabel('distance')

Given the plot of distance versus t above, we can see the distance is minimized when t is near 5. What is the bee's position when t is 5?

In [ ]:
print(make_bee_track(5))

We can check back to the xy plot to see, indeed, this point is pretty close to the flower.

What if we want to know, however, exactly the closest point? There are several ways to find this. Here we are going to use a "brute force" approach which will work on many different problems. Specifically, we will use scipy.optimize.minimize_scalar. The overall idea of this kind of numerical optimization is that we find the best fitting parameter to minimize our "error".

In this example, we are not so much concerned with the exact algorithm being used, but with the way we call this algorithm.

Let's create a function which uses the flower location to compute distance:

In [ ]:
def calc_distance_func(t):
    # assume variable 'flower' in our global scope
    x1, y1 = flower
    # calculate bee position (also assuming 'make_bee_track' in scope)
    x2, y2 = make_bee_track(t)
    dist = compute_distance((x1,y1), (x2,y2))
    print(f't: {t} -> dist: {dist}')
    return dist
In [ ]:
calc_distance_func(0)
In [ ]:
calc_distance_func(5)
In [ ]:
import scipy.optimize
result = scipy.optimize.minimize_scalar(calc_distance_func)
In [ ]:
scipy.optimize.minimize_scalar?
In [ ]:
result
In [ ]:
type(result)
In [ ]:
result.x
In [ ]:
print(calc_distance_func(5))
In [ ]:
calc_distance_func(5.468493134264144)
In [ ]:
# Where is the bee for this value of `t`?
make_bee_track(5.468493134264144)
In [ ]:
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4,3))
ax.plot( [flower[0]], [flower[1]], 'or', label='flower' )
ax.plot( bee_track[0], bee_track[1], '.-k', label='bee')

x,y = make_bee_track(5.468493134264144)
ax.plot( [x], [y] ,'x', label='closest')
ax.axis('equal')
ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('y')
</html>