以最小误差求解非线性方程组

Question

我有一个非线性优化问题（最好用Python来解决）：

在一个二维平面上，有三个圆（它们的中心坐标是x1到x3，y1到y3，半径是d1到d3）。

(x-x1)^2 + (y-y1)^2 - r1^2 = 0
(x-x2)^2 + (y-y2)^2 - r2^2 = 0
(x-x3)^2 + (y-y3)^2 - r3^2 = 0

我们想要找到一个共同的点（x/y），在这种情况下可以用fsolve（scipy.optimize）来计算。但是，如果这三个圆的半径r1到r3都有一定的不确定性u1到u3，该怎么解决这个问题呢？也就是说，圆的真实半径是在r-u到r+u这个范围内。

我该如何找到一个最佳的点（x/y），同时考虑到半径的不确定性呢？

Answer 1

我会这样尝试。对于一个给定的点p，我会计算这个点到三个圆的距离。这个计算可以通过取两个值的绝对差来完成：(1) 点到圆心的距离和 (2) 圆的半径。然后，你的目标就是最小化这三个距离的总和（点p到圆A、圆B和圆C的距离）。我会尝试用简单的加法来计算，但可能还有其他更好的方法来汇总这些距离（比如某种平均计算）。接下来，你可以使用非线性优化工具来最小化这个目标函数（也就是这三个距离的总和）。

现在你有了每个圆的权重。所以你需要修改目标函数，让每个圆的距离根据其不确定性进行调整。这里又需要一些逻辑来决定如何进行加权……一种简单的方法是使用“加权平均”，也就是用1/σ²作为每个距离的权重。这样你的目标函数就变成了：

(weighted average distance) 
= ( distA * sigmaA^-2 + distB * sigmaB^-2 + distC * sigmaC^-2 ) 
   / ( sigmaA^-2 + sigmaB^-2 + sigmaC^-2) 

其中distX是点到圆的距离，σA是圆的位置的标准差（因为圆的位置和大小的不确定性），^-2表示平方后再取倒数。

使用非线性工具包来通过改变点的x和y值来最小化上述目标函数。

Answer 2

其实这并不难，参考这里的描述：

一个建议的算法：

1. 根据链接中的说明找到x。
2. 计算y。

这两个步骤可以用任意两个圆来完成。

1. 将点(x,y)和(x,-y)代入第三个圆。
2. 结果要么是这三个圆在以下点相交：
• x,y
• x,-y
3. 要么根本不相交。

还有另一个建议……我再次阅读了你的问题，意识到你并不是在寻找那个点本身……

不过，如果你在纸上画出所有9个圆（3个相交的圆，加上两个更小和更大的圆，分别代表r+e和r-e，其中e是误差），你会注意到以下几点。你的交点位于一个多边形内。你可以很容易地计算出这个多边形的顶点。然后你的问题就变成了：

• 在多边形内找一个点。
• 或者你写一个目标函数来找到这些顶点，然后最小化这个区域。

想看看我说的圆的情况，可以运行：

# excuse me for the ugly code ...
import pylab
pylab.axes()

cir = pylab.Circle((1,0), radius=1, alpha =.2, fc='b')
cir1 = pylab.Circle((1,0), radius=0.9, alpha =.2, fc='b')
cir2 = pylab.Circle((1,0), radius=1.1, alpha =.2, fc='b')
cir3 = pylab.Circle((-1,0), radius=1, alpha =.2, fc='b')
cir4 = pylab.Circle((-1,0), radius=0.9, alpha =.2, fc='b')
cir5 = pylab.Circle((-1,0), radius=1.1, alpha =.2, fc='b')
cir6 = pylab.Circle((0,-1), radius=0.9, alpha =.2, fc='b')
cir7 = pylab.Circle((0,-1), radius=1.1, alpha =.2, fc='b')
cir8 = pylab.Circle((0,-1), radius=1, alpha =.2, fc='b')
pylab.gca().add_patch(cir)
pylab.gca().add_patch(cir1)
pylab.gca().add_patch(cir2)
pylab.gca().add_patch(cir3)
pylab.gca().add_patch(cir4)
pylab.gca().add_patch(cir5)
pylab.gca().add_patch(cir6)
pylab.gca().add_patch(cir7)
pylab.gca().add_patch(cir8)

pylab.axis('scaled')
pylab.show()