Симпи - найти точку, минимизирующую расстояние до трех окружностей

Question

Я экспериментировал с различными методами нахождения точки (x,y), которая минимизирует суммарные расстояния от (x,y) до окружностей трех окружностей.

На рисунке ниже показан пример расположения этих окружностейи позиционирование (x,y) с использованием первых трех из четырех методов.

Сейчас я пробую четвертую технику, описанную в этом стековом посте .Короче говоря, я хотел бы использовать sympy для вычисления производной 1-го порядка функции потерь, чтобы найти ее минимумы / максимумы, а затем вычислить производную 2-го порядка для выделения минимумов.

Этофункция потерь:

$ E (x, y) = \ sum_i \ big ((x-x_i) ^ 2 + (y-y_i) ^ 2 - r_i ^ 2 \ big) ^ 2 $

Это производная 1-го порядка:

$ E '(x, y) = \ sum_i \ frac {y-y_i} {- x + x_i} $

Здесьэто код, который пытается решить производное уравнение 1-го порядка для y в терминах x, а затем подставить его для решения для x:

x, y = sympy.symbols('x y')
x1, y1 = 0, 0
x2, y2 = 3, 0
x3, y3 = 2, 3

def fprime(x,y):
    return (y-y1)/(-x+x1) + (y-y2)/(-x+x2) + (y-y3)/(-x+x3)

sols = sympy.solve(fprime(x,y), y)
y = sols[0]
x_sols = sympy.solve(fprime(x,y), x)
y_sols = []
for x_sol in x_sols:
    y_sols.append(y.evalf(subs={x:x_sol}))

for x,y in zip(x_sols, y_sols):
    plt.scatter(float(x), float(y)) 

Я не верю, что яправильно оцениваю / решаю эти уравнения, потому что сгенерированные точки очень неправильные (см. изображение ниже)

Чтобы продемонстрировать, что градиент не зависитпо радиусу окружностей:

In [2]: x,y = sympy.symbols('x y')

In [3]: xi, yi, ri = sympy.symbols('xi yi ri', constant=True)

In [4]: def f(x,y):
   ...:     return ((x-xi)**2 + (y-yi)**2 - ri**2)**2
   ...:
   ...:

In [5]: sympy.idiff(f(x,y), x, y)
Out[5]: (y - yi)/(-x + xi)

Answer 1

idiff не вычисляет градиент.То, что вы хотите, это sympy.vector.gradient.Я рассмотрю это только для одной пары ( x_i, y_i ):

import sympy
from sympy.vector import CoordSys3D, gradient

xi, yi, ri = sympy.symbols('xi yi ri', constant=True)
R = CoordSys3D('R')
f1 = ((R.x-xi)**2 + (R.y-yi)**2 - ri**2)**2
fprime = gradient(f1)

, что дает вам

In [25]: fprime
Out[25]: ((4*R.x - 4*xi)*(-ri**2 + (R.x - xi)**2 + (R.y - yi)**2))*R.i + ((4*R.y - 4*yi)*(-ri**2 + (R.x - xi)**2 + (R.y - yi)**2))*R.j

Это двухмерное векторное поле.Теперь мы хотим найти R.x и R.y так, чтобы оба компонента были равны нулю.Итак, сначала я решаю компонент R.i для нуля:

sympy.solve(fprime.components[R.i], R.x)

, что дает

[xi,
 xi - sqrt((-R.y + ri + yi)*(R.y + ri - yi)),
 xi + sqrt((-R.y + ri + yi)*(R.y + ri - yi))]

Я просто выбрал одно решение из этого, позже вы можете проверить, действительно ли это минимум.Итак, теперь нам нужно подключить это к компоненту R.j, чтобы получить уравнение для R.y:

eq = fprime.components[R.j].subs(R.x, xi - sympy.sqrt((-R.y + ri + yi)*(R.y + ri - yi)))

Решение этого с помощью

sympy.solve(eq3, R.y)

дает просто y_i, поэтомуу нас есть наше решение.

Я надеюсь, что это обрисовывает в общих чертах то, что вы должны сделать, и что я не допустил никаких ошибок, так как я не являюсь хорошим экспертом.Я открыт для конструктивной критики здесь.