时间戳

首页 > 解决方案 > 使用 Numpy linalg.lstsq 求解线性系统时获得（显着）不准确的值

python - 使用 Numpy linalg.lstsq 求解线性系统时获得（显着）不准确的值

问题描述

我正在尝试使用Numpy'slinalg.lstsq()来遵循这篇文章中的逻辑。我会在这里重复一遍。

假设我在 3 维空间中有四个麦克风（或光接收器等）。我想使用每个麦克风记录枪声之间的时间差，或者到达每个麦克风的时间，来确定射手的坐标。因此，我们知道了四个麦克风的坐标x_i, y_i，以及声音的速度v。我们不知道射击者的坐标，也不知道射击的时间。

对于每个麦克风，我们写：

(X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 = (v(t_i - T))**2

X,Y,Z射手的坐标在哪里，T是射门的时间。

i/j如果我们考虑所有可能性，i从方程中减去方程，这个问题就可以简化j。我们最终得到以下形式的6方程式（通常，给定n麦克风，n(n-1)/2方程式）：

2*(x_j - x_i)*X + 2*(y_j - y_i)*Y + 2*(z_j - z_i)*Z + 2 * v**2 * (t_i - t_j) = 2 * v**2 ( t_i**2 - t_j**2) + (x_j**2 + y_j**2 + z_j**2) - (x_i**2 + y_i**2 + z_i**2)

这些方程的形式是Xv_1 + Y_v2 + Z_v3 + T_v4 = b，其中v_i是向量。然后，我们（应该）能够使用线性最小二乘拟合来“解决”（或至少粗略估计）X,Y,Z，并且作为副产品，T.

我已经尝试了很多不同的方法来解决这个问题，并且所有迹象都指向使用链接帖子中描述的方法（即，通过线性最小二乘法解决这个问题，然后使用这个“解决方案”作为更多“精确”算法）。因此，我编写了以下（诚然，有些粗糙）Python代码来开始。它解决了问题的线性最小二乘部分：

from dataclasses import dataclass
from numpy import *
from random import randrange

@dataclass
class Vertexer:

    @staticmethod
    def find():

        # Pick microphones to be at random locations
        x_1 = randrange(10); y_1 = randrange(10); z_1 = randrange(10)
        x_2 = randrange(10); y_2 = randrange(10); z_2 = randrange(10)
        x_3 = randrange(10); y_3 = randrange(10); z_3 = randrange(10)
        x_4 = randrange(10); y_4 = randrange(10); z_4 = randrange(10)

        # Pick shooter to be at random location
        x = randrange(100); y = randrange(100); z = randrange(100)

        # Set velocity (ok, it's a ray gun...)
        c = 299792 # km/s

        # Generate simulated source
        t_1 = math.sqrt( (x - x_1)**2 + (y - y_1)**2 + (z - z_1)**2 ) / c
        t_2 = math.sqrt( (x - x_2)**2 + (y - y_2)**2 + (z - z_2)**2 ) / c
        t_3 = math.sqrt( (x - x_3)**2 + (y - y_3)**2 + (z - z_3)**2 ) / c
        t_4 = math.sqrt( (x - x_4)**2 + (y - y_4)**2 + (z - z_4)**2 ) / c


        A = array([[2 * (x_1 - x_2), 2 * (y_1 - y_2), 2 * (z_1 - z_2), 2 * c * c * (t_2 - t_1)],
                   [2 * (x_1 - x_3), 2 * (y_1 - y_3), 2 * (z_1 - z_3), 2 * c * c * (t_3 - t_1)],
                   [2 * (x_1 - x_4), 2 * (y_1 - y_4), 2 * (z_1 - z_4), 2 * c * c * (t_4 - t_1)],
                   [2 * (x_2 - x_3), 2 * (y_2 - y_3), 2 * (z_2 - z_3), 2 * c * c * (t_3 - t_2)],
                   [2 * (x_2 - x_4), 2 * (y_2 - y_4), 2 * (z_2 - z_4), 2 * c * c * (t_4 - t_2)],
                   [2 * (x_3 - x_4), 2 * (y_3 - y_4), 2 * (z_3 - z_4), 2 * c * c * (t_4 - t_3)]
                   ])
        b = array([2 * c * c * (t_2 * t_2 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_2 * x_2 + y_2 * y_2 + z_2 * z_2),
                   2 * c * c * (t_3 * t_3 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
                   2 * c * c * (t_4 * t_4 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
                   2 * c * c * (t_3 * t_3 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
                   2 * c * c * (t_4 * t_4 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
                   2 * c * c * (t_4 * t_4 - t_3 * t_3) + (x_3 * x_3 + y_3 * y_3 + z_3 * z_3) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4)])

        solved, residuals, rank, s = linalg.lstsq(A, b)

        print(solved)


myVertexer = Vertexer()

myVertexer.find()

然而不幸的是，它似乎并没有产生非常合理的估计。如果愿意，您可以自己运行代码，或者这里是运行的集合，将找到的值与算法产生的“模拟”源的坐标进行比较。

Predicted: [ 2.26739519e+00  4.80191502e+00 -5.07181020e+00 ]
Actual:    78 35 57

Predicted: [ 7.72173135e-01 -2.61250803e+00  4.10306750e+00 ]
Actual:    30 75 48

Predicted: [-1.65368110e+00  8.82919123e-01  1.43648336e+00 ]
Actual:    67 56 44

Predicted: [ 4.08698715e+00 -3.19377148e-01  8.81706364e-01 ]
Actual:    9 82 13

诚然，我不确定从最初的最小二乘解决方案中期望得到什么样的精度。但是，这些值对我来说似乎不是很好。我哪里错了？

请让我知道我是否可以提供任何进一步的澄清。

请注意，我知道还有其他方法可以解决此问题（值得注意的是，只需插入带有非线性求解器的多个软件包之一即可获得非常好的解决方案）。我有兴趣确定为什么这个解决方案（我知道链接的海报已经成功使用，还有其他几个）不起作用，既是出于好奇，也是因为我想更进一步，其中“自定义”代码将很有用。

标签： pythonnumpylinear-algebraphysics

解决方案

有两个单独的问题：

  1. 您在原始帖子中引用的减法方程有一个错误：在 RHS 上它应该是c * c * (t_i * t_i - t_j * t_j)而不是2 * c * c * (t_i * t_i - t_j * t_j).
  2. 第二个问题是，在矩阵中A你有一列tau（第四列），而在你的计算中t_1你t_4隐含地假设tau=0，所以你不能把它作为矩阵中的变量A；因此需要删除第四列。

修改后的代码如下：

from random import randrange, seed
from numpy import array, linalg, set_printoptions
from dataclasses import dataclass
set_printoptions(suppress=True, linewidth=1000, precision=6)
seed(5)

@dataclass
class Vertexer:

    @staticmethod
    def find():

        # Pick microphones to be at random locations
        x_1 = randrange(10); y_1 = randrange(10); z_1 = randrange(10)
        x_2 = randrange(10); y_2 = randrange(10); z_2 = randrange(10)
        x_3 = randrange(10); y_3 = randrange(10); z_3 = randrange(10)
        x_4 = randrange(10); y_4 = randrange(10); z_4 = randrange(10)

        # Pick shooter to be at random location
        x = randrange(100); y = randrange(100); z = randrange(100)

        # Set velocity (ok, it's a ray gun...)
        c = 299792 # km/ns

        # Generate simulated source
        t_1 = math.sqrt( (x - x_1)**2 + (y - y_1)**2 + (z - z_1)**2 ) / c
        t_2 = math.sqrt( (x - x_2)**2 + (y - y_2)**2 + (z - z_2)**2 ) / c
        t_3 = math.sqrt( (x - x_3)**2 + (y - y_3)**2 + (z - z_3)**2 ) / c
        t_4 = math.sqrt( (x - x_4)**2 + (y - y_4)**2 + (z - z_4)**2 ) / c

        A = array([[2 * (x_1 - x_2), 2 * (y_1 - y_2), 2 * (z_1 - z_2)],
                   [2 * (x_1 - x_3), 2 * (y_1 - y_3), 2 * (z_1 - z_3)],
                   [2 * (x_1 - x_4), 2 * (y_1 - y_4), 2 * (z_1 - z_4)],
                   [2 * (x_2 - x_3), 2 * (y_2 - y_3), 2 * (z_2 - z_3)],
                   [2 * (x_2 - x_4), 2 * (y_2 - y_4), 2 * (z_2 - z_4)],
                   [2 * (x_3 - x_4), 2 * (y_3 - y_4), 2 * (z_3 - z_4)]
                   ])
        b = array([1 * c * c * (t_2 * t_2 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_2 * x_2 + y_2 * y_2 + z_2 * z_2),
                   1 * c * c * (t_3 * t_3 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
                   1 * c * c * (t_4 * t_4 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
                   1 * c * c * (t_3 * t_3 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
                   1 * c * c * (t_4 * t_4 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
                   1 * c * c * (t_4 * t_4 - t_3 * t_3) + (x_3 * x_3 + y_3 * y_3 + z_3 * z_3) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4)])

        solved, residuals, rank, s = linalg.lstsq(A, b, rcond=None)

        print('Actual:', [x, y, z])
        print('Predicted:', solved)

myVertexer = Vertexer()

myVertexer.find()

输出是：

Actual: [31, 48, 69]
Predicted: [31. 48. 69.]

推荐阅读