geometrycurvespython

How do I compute a python function which calculates/plot minimum distance between a point and a curve?

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fmin_cobyla
P = (2487,1858)  
def f(x):  
    return (1.77425666496e-05 * x**3 + -0.125555135823 * x**2 + 296.656015233 * x + -232082.382985)
def objective(X):
    x,y = X 
    return np.sqrt((x - P[0])**2 + (y - P[1])**2)
def c1(X):
    x,y = X
    return f(x) - y
X = fmin_cobyla(objective, x0=[0.1, 0.01], cons=[c1])
print ("The minimum distance is {0:1.2f}".format(objective(X)))
v1 = np.array(P) - np.array(X)
v2 = np.array([1, 2.0 * X[0]])
print ("dot(v1, v2) = ",np.dot(v1, v2))
x = np.linspace(1000, 3000, 500000)
plt.plot(x, f(x), 'r-', label='f(x)')
plt.plot(P[0], P[1], 'bo', label='point')
plt.plot([P[0], X[0]], [P[1], X[1]], 'b-', label='shortest distance')
plt.plot([X[0], X[0] + 1], [X[1], X[1] + 2.0 * X[0]], 'g-', label='tangent')
plt.axis('equal')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

The figure that shows computing distance which is obviously wrong. Can someone tell me what's wrong with my code ?

Solution

  • I changed the problem from a two-dimensional optimisation problem to a one-dimensional one:

    import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fmin_cobyla


P = (2487,1858)

def f(x):
    return (1.77425666496e-05 * x**3 + -0.125555135823 * x**2 + 296.656015233 * x + -232082.382985)

def objective(x):
    return np.sqrt((x[0] - P[0])**2 + (f(x[0]) - P[1])**2)

def c1(x):
    return 0.0


x_c = fmin_cobyla(objective, x0=2300, cons=c1)
X = [x_c , f(x_c)]

print ("The minimum distance is {0:1.2f}".format(objective(X)))


v1 = np.array(P) - np.array(X)
v2 = np.array([1, 2.0 * X[0]])

print ("dot(v1, v2) = ",np.dot(v1, v2))

x = np.linspace(1000, 3000, 500000)


plt.plot(x, f(x), 'r-', label='f(x)')

plt.plot(P[0], P[1], 'bo', label='point')

plt.plot([P[0], X[0]], [P[1], X[1]], 'b-', label='shortest distance')

plt.plot([X[0], X[0] + 1], [X[1], X[1] + 2.0 * X[0]], 'g-', label='tangent')
plt.axis('equal')

plt.xlabel('x')

plt.ylabel('y')

plt.show()