Прямо сейчас мы пытаемся создать 5-мерный фильтр Калмана, который получает на вход координаты x и y маленького жука, который бродит по коробке. Мы наблюдаем за этой ошибкой, когда она совершает примерно 2000 движений, а затем прогнозируем ее положение с этого момента. Размеры: x-координата, y-координата, скорость, курс и угловое ускорение. Ниже приведен код, который у нас есть до сих пор.
#x is a list of five variables - x, y, velocity, angularVelocity, angularAcceleration
#currentmeasurement is the x and y that were observed
def kalman_filter(x, P, currentmeasurement, lastMeasurement = None):
prevmeasurement = []
#if there is a lastMeasurement argument, it becomes measurement
if lastMeasurement:
prevmeasurement = [lastMeasurement[0], lastMeasurement[1], x.value[3][0]]
#if there is no lastMeasurement argument, the current measurement becomes measurement.
else:
prevmeasurement = [x.value[0][0], x.value[1][0], x.value[3][0]]
#Prediction Step
a = x.value[3][0]
F = matrix([
[1., 0., cos(a), 0., 0.],
[0., 1., sin(a), 0., 0.],
[0., 0., 1., 0., 0.],
[0., 0., 0., 1., dt],
[0., 0., 0., 0., 1.]])
x = (F * x) + u
P = F * P * F.transpose()
#we can calculate the heading from our observations.
heading = atan2(currentmeasurement[0][1] - prevmeasurement[1], currentmeasurement[0][0] - prevmeasurement[0])
while(abs(heading - prevmeasurement[2]) > pi):
heading = heading + 2*pi*((prevmeasurement[2]-heading)/abs(prevmeasurement[2]-heading))
#perhaps the velocity should also be calculated?
# measurement update
dt = 1.
u = matrix([[0.], [0.], [0.], [0.], [0.]]) # external motion
H = matrix([[1., 0., 0., 0., 0.], #the measurement function
[0., 1., 0., 0., 0.],
[0., 0., 0., 1., 0.]])
R = matrix([[1., 0., .0], #measurement uncertainty
[0., 1., .0],
[0., 0., 1.]])
I = matrix([[]]) #a 5x5 identity matrix
I.identity(5)
prevmeasurement = [currentmeasurement[0][0], currentmeasurement[0][1], heading]
Z = matrix([prevmeasurement])
y = Z.transpose() - (H * x)
S = H * P * H.transpose() + R
K = P * H.transpose() * S.inverse()
x = x + (K * y)
P = (I - (K * H)) * P
return x,P
Это приводит к крайне неверным оценкам. Мы не уверены, что мы здесь делаем неправильно - я думаю, что мы правильно выполняем все шаги, но не уверен, что у нас есть все необходимые матрицы. Любой вклад будет полезен!