使用NumPy最小化此误差函数

本教程将介绍使用NumPy最小化此误差函数的处理方法,这篇教程是从别的地方看到的,然后加了一些国外程序员的疑问与解答,希望能对你有所帮助,好了,下面开始学习吧。

使用NumPy最小化此误差函数 教程 第1张

问题描述

背景

我在3-dimensions和使用4节点中尝试解决(出了名的痛苦的)到达时间差(TDOA)多边问题已经有一段时间了。如果您不熟悉这个问题,它是在给定n节点的坐标、信号到达每个节点的时间以及信号v的速度的情况下,确定某个信号源(X,Y,Z)的坐标。

我的解决方案如下:

对于每个节点,我们编写(X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 = (v(t_i - T)**2

其中(x_i, y_i, z_i)ith节点的坐标,T是发射时间。

我们现在有4个未知数中的4个方程。四个节点显然不够。我们尝试直接解决这个系统,然而,考虑到问题的高度非线性性质(事实上,我已经尝试了直接技术),这似乎几乎是不可能的。并且失败了)。,我们通过考虑所有i/j可能性,从等式j中减去i,将其简化为一个线性问题。我们得到(n(n-1))/2=6形式的方程:

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) = 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。现在我们尝试应用标准线性最小二乘法,其中解是A^T Ax = A^T b中的矩阵向量x。不幸的是,如果您尝试将其添加到任何标准的线性最小二乘算法中,它将无法正常工作。那么,我们现在该怎么办?

...

信号到达节点i的时间(当然)由

给出:

sqrt( (X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 ) / v

这个方程意味着到达时间T0。如果我们有T = 0,我们就可以删除矩阵中T列,问题就大大简化了。事实上,NumPy'slinalg.lstsq()给出的结果非常准确。

...

所以,我所做的是输入时间,方法是从每个方程式中减去最早的时间。然后,我所要做的就是确定我可以每次相加的dt,以使通过线性最小二乘找到的点的平方和误差的残差最小化。

我将某些dt的误差定义为将input times+dt提供给最小二乘算法预测的点的到达时间之间的平方差减去所有4节点的输入时间(归一化)。

for node, time in nodes, times:
 error += ( (sqrt( (X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 ) / v) - time) ** 2

我的问题:

通过使用暴力,我能够令人满意地做到这一点。我从dt = 0开始,然后向上移动到某个最大迭代次数,直到达到某个最小的RSS误差,也就是我添加到归一化时间以获得解的dt。得到的解决方案非常准确和精确,但速度很慢。

实际上,我希望能够在解决此问题,因此需要一个速度快得多的解决方案。我首先假设误差函数(即上述定义的dtvserror)将是高度非线性的--即刻,这对我来说是有意义的。

由于我没有实际的数学函数,我可以自动排除需要区分的方法(例如Newton-Raphson)。误差函数将始终为正,因此我可以排除bisection等。相反,我尝试一个简单的近似搜索。不幸的是,这一努力惨遭失败。然后,我尝试了禁忌搜索,然后是遗传算法,以及其他几种搜索算法。他们都惨败了。

所以,我决定做一些调查。事实证明,误差函数与DT的关系图看起来有点像平方根,只是根据信号源到节点的距离而向右移动:

其中dt位于水平轴,错误位于垂直轴

而且,事后看来,。我将误差函数定义为涉及平方根,因此,至少在我看来,这似乎是合理的。

要做什么?

所以,我现在的问题是,我怎么确定误差函数的最小值对应的dt

我的第一次尝试(非常粗略)是在错误图上获取一些点(如上所述),使用numpy.polyfit进行拟合,然后将结果提供给numpy.root。该根对应于dt。不幸的是,这也失败了。我尝试了各种degrees,也尝试了各种点,甚至达到了荒谬的数量,以至于我可能会使用暴力。

怎么确定此误差函数的最小值对应的dt

由于我们处理的是高速(无线电信号),因此和很重要,因为dt中的微小变化可能会偏离结果点。

我相信我在这里所做的事情中隐藏着一些无限简单的方法,忽略其他一切,我怎么找到dt

我的要求:

    速度至上

    我只能访问将在其中运行此程序的环境中的纯PythonNumPy

编辑:

以下是我的代码。诚然,有点凌乱。这里,我使用的是polyfit技术。它将为您模拟&Q;源,并比较结果:

from numpy import poly1d, linspace, set_printoptions, array, linalg, triu_indices, roots, polyfit
from dataclasses import dataclass
from random import randrange
import math

@dataclass
class Vertexer:

 receivers: list

 # Defaults
 c = 299792

 # Receivers:
 # [x_1, y_1, z_1]
 # [x_2, y_2, z_2]
 # [x_3, y_3, z_3]

 # Solved:
 # [x, y, z]

 def error(self, dt, times):
  solved = self.linear([time + dt for time in times])

  error = 0
  for time, receiver in zip(times, self.receivers):
error += ((math.sqrt( (solved[0] - receiver[0])**2 + 
 (solved[1] - receiver[1])**2 +
 (solved[2] - receiver[2])**2 ) / c ) - time)**2

  return error

 def linear(self, times):
  X = array(self.receivers)
  t = array(times)

  x, y, z = X.T
  i, j = triu_indices(len(x), 1)

  A = 2 * (X[i] - X[j])
  b = self.c**2 * (t[j]**2 - t[i]**2) + (X[i]**2).sum(1) - (X[j]**2).sum(1)

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

  return(solved)

 def find(self, times):
  # Normalize times
  times = [time - min(times) for time in times]
  # Fit the error function

  y = []
  x = []
  dt = 1E-10
  for i in range(50000):
x.append(self.error(dt * i, times))
y.append(dt * i) 

  p = polyfit(array(x), array(y), 2)
  r = roots(p)
  return(self.linear([time + r for time in times]))




# SIMPLE CODE FOR SIMULATING A SIGNAL

# Pick nodes 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 source to be at random location
x = randrange(1000); y = randrange(1000); z = randrange(1000)

# Set velocity
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

print('Actual:', x, y, z)

myVertexer = Vertexer([[x_1, y_1, z_1],[x_2, y_2, z_2],[x_3, y_3, z_3],[x_4, y_4, z_4]])
solution = myVertexer.find([t_1, t_2, t_3, t_4])
print(solution)

推荐答案

似乎班克罗夫特方法适用于此问题?下面是一个纯NumPy实现。

# Implementation of the Bancroft method, following
# https://gssc.esa.int/navipedia/index.php/Bancroft_Method
M = np.diag([1, 1, 1, -1])


def lorentz_inner(v, w):
 return np.sum(v * (w @ M), axis=-1)


B = np.array(
 [
  [x_1, y_1, z_1, c * t_1],
  [x_2, y_2, z_2, c * t_2],
  [x_3, y_3, z_3, c * t_3],
  [x_4, y_4, z_4, c * t_4],
 ]
)
one = np.ones(4)
a = 0.5 * lorentz_inner(B, B)
B_inv_one = np.linalg.solve(B, one)
B_inv_a = np.linalg.solve(B, a)
for Lambda in np.roots(
 [
  lorentz_inner(B_inv_one, B_inv_one),
  2 * (lorentz_inner(B_inv_one, B_inv_a) - 1),
  lorentz_inner(B_inv_a, B_inv_a),
 ]
):
 x, y, z, c_t = M @ np.linalg.solve(B, Lambda * one + a)
 print("Candidate:", x, y, z, c_t / c)

好了关于使用NumPy最小化此误差函数的教程就到这里就结束了,希望趣模板源码网找到的这篇技术文章能帮助到大家,更多技术教程可以在站内搜索。