时间戳

首页 > 解决方案 > Sympy - 找到将距离最小化为三个圆周的点

python - Sympy - 找到将距离最小化为三个圆周的点

问题描述

我一直在尝试各种技术来找到一个点，以最小化从三个圆到圆周的总(x,y)距离。(x,y)

下图显示了这些圆圈的示例排列以及(x,y)使用四种技术中的前三种的定位。

演示

我现在正在尝试此 stackexchange 帖子中描述的第四种技术。简而言之，我想用sympy计算损失函数的一阶导数来找到它的最小值/最大值，然后计算二阶导数来隔离最小值。

这是损失函数：

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

这是一阶导数：

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

这是尝试求解一阶导数方程的代码，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)

标签： pythonsympy

解决方案

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的，所以我们有我们的解决方案。

我希望这概述了您必须做的事情，并且我没有犯任何错误，因为我不是 sympy 专家。我愿意在这里接受建设性的批评。

推荐阅读