предсказание фильтра Калмана t + 2 временного шага

Question

У меня есть вопрос об использовании фильтра Калмана для прогнозирования значений t + 2. Как мы знаем, базовый фильтр Калмана имеет два этапа: прогнозирование и обновление. Часть предсказания может генерировать xt на основе xt-1. Вот несколько примеров кодов, которые я нашел в Интернете.

import numpy as np

class KalmanFilter(object):
    def __init__(self, F = None, B = None, H = None, Q = None, R = None, P = None, x0 = None):

        if(F is None or H is None):
            raise ValueError("Set proper system dynamics.")

        self.n = F.shape[1]
        self.m = H.shape[1]

        self.F = F
        self.H = H
        self.B = 0 if B is None else B
        self.Q = np.eye(self.n) if Q is None else Q
        self.R = np.eye(self.n) if R is None else R
        self.P = np.eye(self.n) if P is None else P
        self.x = np.zeros((self.n, 1)) if x0 is None else x0

    def predict(self, u = 0):
        self.x = np.dot(self.F, self.x) + np.dot(self.B, u)
        self.P = np.dot(np.dot(self.F, self.P), self.F.T) + self.Q
        return self.x

    def update(self, z):
        y = z - np.dot(self.H, self.x)
        S = self.R + np.dot(self.H, np.dot(self.P, self.H.T))
        K = np.dot(np.dot(self.P, self.H.T), np.linalg.inv(S))
        self.x = self.x + np.dot(K, y)
        I = np.eye(self.n)
        self.P = np.dot(np.dot(I - np.dot(K, self.H), self.P), 
            (I - np.dot(K, self.H)).T) + np.dot(np.dot(K, self.R), K.T)

def example():
    dt = 1.0/60
    F = np.array([[1, dt, 0], [0, 1, dt], [0, 0, 1]])
    H = np.array([1, 0, 0]).reshape(1, 3)
    Q = np.array([[0.05, 0.05, 0.0], [0.05, 0.05, 0.0], [0.0, 0.0, 0.0]])
    R = np.array([0.5]).reshape(1, 1)

    x = np.linspace(-10, 10, 100)
    measurements = - (x**2 + 2*x - 2)  + np.random.normal(0, 2, 100)

    kf = KalmanFilter(F = F, H = H, Q = Q, R = R)
    predictions = []

    for z in measurements:
        predictions.append(np.dot(H,  kf.predict())[0])
        kf.update(z)

    import matplotlib.pyplot as plt
    plt.plot(range(len(measurements)), measurements, label = 'Measurements')
    plt.plot(range(len(predictions)), np.array(predictions), label = 'Kalman Filter Prediction')
    plt.legend()
    plt.show()

if __name__ == '__main__':
    example()

В этой задаче мы используем значение t-1 для прогнозирования t и обновляем значение t. Если я хочу предсказать значение t + 1 на основе т. Я что-то изменил соответственно:

import numpy as np

class KalmanFilter(object):
    def __init__(self, F = None, F_1 = None, B = None, H = None, Q = None, R = None, P = None, x0 = None):

        if(F is None or H is None):
            raise ValueError("Set proper system dynamics.")

        self.n = F.shape[1]
        self.m = H.shape[1]

        self.F = F
        self.F_1 = F_1
        self.H = H
        self.B = 0 if B is None else B
        self.Q = np.eye(self.n) if Q is None else Q
        self.R = np.eye(self.n) if R is None else R       
        self.P = np.eye(self.n) if P is None else P
        self.x = np.zeros((self.n, 1)) if x0 is None else x0

    def predict(self, u = 0):
        self.x = np.dot(self.F, self.x) + np.dot(self.B, u)
        self.P = np.dot(np.dot(self.F, self.P), self.F_1) + self.Q
        return self.x

    def update(self, z):
        y = z - np.dot(self.H, self.x)
        S = self.R + np.dot(self.H, np.dot(self.P, self.H.T))
        K = np.dot(np.dot(self.P, self.H.T), np.linalg.inv(S))
        self.x = self.x + np.dot(K, y)
        I = np.eye(self.n)
        self.P = np.dot(np.dot(I - np.dot(K, self.H), self.P), 
                        (I - np.dot(K, self.H)).T) + np.dot(np.dot(K, self.R), K.T)


def example():
    dt = 1.0/60
    F_0 = np.array([[1, dt, 0], [0, 1, dt], [0, 0, 1]])
    F = np.dot(F_0, F_0)
    F_1 = np.dot(F_0.T, F_0.T)
    H = np.array([1, 0, 0]).reshape(1, 3)
    Q = np.array([[0.05, 0.05, 0.0], [0.05, 0.05, 0.0], [0.0, 0.0, 0.0]])
    R = np.array([0.5]).reshape(1, 1)

    x = np.linspace(-10, 10, 100)
    measurements = - (x**2 + 2*x - 2)  + np.random.normal(0, 2, 100)

    kf = KalmanFilter(F = F, F_1 = F_1, H = H, Q = Q, R = R)
    predictions = []

    for i in range(1, len(measurements), 2):
            predictions.append(np.dot(H,  kf.predict())[0])
            kf.update(measurements[i])

    import matplotlib.pyplot as plt
    plt.plot(range(len(measurements)), measurements, label = 'Measurements')
    plt.plot(range(len(predictions)), np.array(predictions), label = 'Kalman Filter Prediction')
    plt.legend()
    plt.show()

if __name__ == '__main__':
    example()

Основные изменения - это два:

1. Я изменил матрицу F.

2. Я использовал значение t + 1 временный шаг, чтобы обновить мой результат.

Однако длина результатов, которые я получил, составляет только половину исходных измерений. Из-за того, что я как бы прыгаю, чтобы обновить их.

Я немного растерялся. У кого-нибудь есть предложения или решения? Большое вам спасибо

Answer 1

Я вижу следующие проблемы:

1. Прогноз x(t+1) на основе x(t) фактически совпадает с прогнозом x(t) на основе x(t-1).Все зависит от определения временного шага (dt)
2. Вам не нужно менять F-матрицу, это было правильно в первом коде.Но это зависит от вашего dt.
3. В вашем втором коде вы заменили self.F.T на self.F_1.Но T означает transpose из F.Это не хорошо.Вот почему, вероятно, ваш результирующий вектор имеет другие измерения.

Так что, если вы хотите прогнозировать в течение некоторого промежутка времени, вам нужно знать только это количество времени.Это будет ваш новый dt, который, конечно, изменит ваши F и Q.

Если ваши измерения будут выполнены с шагом по времени dt, и вы хотите посмотреть, что произойдет, если каждыйВторое измерение падает, у вас есть две опции:

Опция 1

Изменение dt шага прогнозирования, чтобы оно равнялось разнице во времени между двумя последними измерениями ипересчитайте матрицу F с новым значением dt (учтите: в этом случае вам нужно будет изменить матрицу Q, потому что система приобретает большую неопределенность при большем временном шаге).

Вариант 2

Выполнить прогнозирование с исходным шагом по времени dt несколько раз, чтобы заполнить промежуток времени между измерениями.Поэтому, если вы хотите обрабатывать только каждое второе измерение, ваш код будет выглядеть так:

for i in range(1, len(measurements), 2):
    kf.predict()
    kf.predict()
    kf.update(measurements[i])

ОБНОВЛЕНИЕ

К вашему вопросу из комментария:

i       t      todo

1    0*dt      init
2    1*dt      predict
3    2*dt      predict, update
4    3*dt      predict
5    4*dt      predict, update

Вы имели в виду этот случай?

Посмотрите на этот пост .Он показывает, что произойдет, если вы многое прогнозируете без новых обновлений.