Python - Kalman Filter (Unscented) not working corretly

I am trying to implement a Kalman Filter which is used to filter calculated values of Speed, Heading, Acceleration and Yaw Rate. I also use the Kalman Filter to predict future values and use those to calculate future geographical positions of an object.

It seems that the Filter is resetting the predicted values to 0. Currently, my next state update is always maintained (nv = v, na = a, etc..), but if I change to nv = v + a * dt I get similar behaviour, values are very close to zero (~ e-14).

Main code:

# Initialize or retrieve the UKF instance
ukf = get_ukf(station_id)
# Update UKF with current raw measurements
ukf.update([speed, heading, acceleration, yaw_rate])

new_predictions = predict_future_positions(
    station_id,
    data[station_id].at[curr_index, 'latitude'], 
    data[station_id].at[curr_index, 'longitude']
)

Prediction code:

def fx(state, dt):
    """ Predict next state based on simple motion model assumptions """
    v, theta, a, psi_dot = state
    nv = v
    ntheta = theta          # Assuming constant heading for simplicity
    na = a                  # Assuming constant acceleration for simplicity
    npsi_dot = psi_dot      # Assuming constant yaw rate for simplicity
    return np.array([nv, ntheta, na, npsi_dot])


def hx(state):
    """ Measurement function returns the state itself, assumes you measure all states directly """
    return state


def initialize_ukf():
    """ Initialize and return a UKF for tracking speed, heading, acceleration, and yaw rate """
    sigmas = MerweScaledSigmaPoints(n=4, alpha=0.1, beta=2., kappa=1.)
    ukf = UKF(dim_x=4, dim_z=4, fx=fx, hx=hx, dt=1/30 , points=sigmas)
    ukf.x = np.array([0., 0., 0., 0.])  # initial state [speed, heading, acceleration, yaw_rate]
    ukf.P *= 10  # initial covariance
    ukf.R = np.eye(4) * 0.1  # measurement noise
    ukf.Q = np.eye(4) * 0.1  # process noise
    return ukf


def get_ukf(station_id):
    """ Retrieve or initialize a UKF for a specific station_id """
    if station_id not in ukfs:
        ukfs[station_id] = initialize_ukf()
    return ukfs[station_id]


def calculate_next_position(lat: float, lon: float, speed: float, heading: float, acceleration: float, interval) -> tuple:
    # Convert latitude, longitude and heading from degrees to radians
    lat_rad, lon_rad, heading_rad = map(math.radians, [lat, lon, heading])

    # Distance covered in the interval with the updated speed
    distance = speed * interval + 0.5 * acceleration * interval ** 2

    # Calculate new latitude and longitude in radians
    new_lat_rad = math.asin(math.sin(lat_rad) * math.cos(distance / EARTH_RADIUS) +
                            math.cos(lat_rad) * math.sin(distance / EARTH_RADIUS) * math.cos(heading_rad))
    new_lon_rad = lon_rad + math.atan2(math.sin(heading_rad) * math.sin(distance / EARTH_RADIUS) * math.cos(lat_rad),
                                       math.cos(distance / EARTH_RADIUS) - math.sin(lat_rad) * math.sin(new_lat_rad))

    # Convert new latitude and longitude from radians to degrees
    new_lat, new_lon = map(math.degrees, [new_lat_rad, new_lon_rad])

    return new_lat, new_lon


# def predict_future_positions(lat, lon, speed, heading, acceleration, curve_info, frequency=30, duration=5) -> list:
def predict_future_positions(station_id, lat, lon, frequency=30, duration=5) -> list:
    predicted_positions = []

    # Calculate the number of points and the duration of each interval
    points = frequency * duration
    interval_duration = 1 / frequency

    ukf = get_ukf(station_id)
    ukf.dt = interval_duration          # Set UKF's dt to the interval duration
    print(f"interval_duration: {ukf.dt}")

    for i in range(1, points + 1):        
        # Predict the next state using the UKF for the current interval
        ukf.predict()
        new_speed, new_heading, new_acceleration, new_yaw_rate = ukf.x
        print(f"Predicted: {new_speed}, {new_heading}, {new_acceleration}, {new_yaw_rate}")

        new_lat, new_lon = calculate_next_position(lat, lon, new_speed, new_heading, new_acceleration, interval_duration)
        predicted_positions.append({'lat': new_lat, 'lon': new_lon})
        lat, lon = new_lat, new_lon

    return predicted_positions

My output:

Speed: 10.254806246670183, Heading: 319.27693339350554, Acceleration: -0.23540827984615223, Yaw rate: -0.005145571837527121
interval_duration: 0.03333333333333333
Predicted: 10.166326468873194, 316.5221712417051, -0.2333771471468742, -0.005145058097190569
Predicted: 10.166326468873166, 316.52217124170784, -0.2333771471468764, -0.005145058097190472
...

Speed: 10.467646761934825, Heading: 319.23780992232264, Acceleration: 0.21254958566790286, Yaw rate: -0.03906999369681518
interval_duration: 0.03333333333333333
Predicted: 10.449913038293573, 319.0779853536578, 0.18630528143125535, -0.03707439752961145
Predicted: 10.449913038293431, 319.0779853536569, 0.1863052814312549, -0.037074397529612
...

Prediction accuracy: 0.10094114807863384 meters
Speed: 10.189528933971618, Heading: 319.25287353184626, Acceleration: -0.2778124154209819, Yaw rate: 0.015047067558968702
Predicted: 10.204898856119058, 319.2425502531232, -0.2504165308446913, 0.011970453237527212
Predicted: 10.204898856119286, 319.2425502531214, -0.25041653084469173, 0.011970453237527268
...

Prediction accuracy: 0.808794667925268 meters
Speed: 10.413539616314386, Heading: 322.54565576546065, Acceleration: 0.22384171917679088, Yaw rate: 3.2902986069174447
Predicted: 10.401223867329037, 322.3506784353949, 0.19584697089942216, 3.0967838480884247
Predicted: 10.401223867329037, 322.3506784353949, 0.19584697089942038, 3.0967838480884033
Predicted: 10.401223867329037, 322.3506784353949, 0.19584697089941927, 3.0967838480884176
...

Prediction accuracy: 0.5443148460430843 meters
Speed: 10.539412327559626, Heading: 323.6073064276114, Acceleration: 0.12570394929129958, Yaw rate: 1.0602272699128528
Predicted: 0.0, 0.0, 0.0, 0.0
Predicted: 0.0, 0.0, 0.0, 0.0
...

Prediction accuracy: 10.36880732693119 meters
Speed: 10.361936129429777, Heading: 335.33964303577073, Acceleration: -0.17724217428198916, Yaw rate: 11.71686610233176
...

Upvotes: 0

Views: 81

Answers (1)

Dawson Beatty
Dawson Beatty

Reputation: 688

(some caveats) I only have limited information about what you are trying to do. The code that you posted also does not reproduce your exact problem, so there may be other things going on in your code that I don't know about.

The kalman filter is useful when you have noisy dynamics and noisy measurements, and are trying to get your best estimate of the current state. In your filter you initialize the estimated state (which is different than the true state) as all zeros, then propagate it (which keeps all of the states constant in your case). I'm not surprised that the propagated state remains zero; that's the expected behavior. The biggest thing that you are missing is a call to ukf.update with a measurement of the true state. If you never update the state, then the filter will just keep propagating the initial state and never do any better than that.

Why don't you have latitude and longitude as part of your state?

Upvotes: 0

Related Questions