Advanced Filters Comparison

This example compares advanced filtering techniques for challenging nonlinear problems.

Overview

When standard Kalman filters are insufficient, advanced techniques provide better performance:

1. Constrained EKF - Enforces state constraints during estimation

2. Gaussian Sum Filter - Represents multi-modal distributions

3. Rao-Blackwellized Particle Filter - Combines analytic and Monte Carlo methods

Key Concepts

• State constraints: Physical bounds on state variables

• Multi-modality: Distributions with multiple peaks

• Hybrid filters: Combining different estimation techniques

• Marginalization: Analytically integrating out linear states

Code Highlights

The example demonstrates:

• Implementing state constraints in EKF updates

• Gaussian mixture representation and merging

• Rao-Blackwellization for linear substructure

• Performance comparison metrics

Source Code

"""Advanced filters comparison demonstration.

Demonstrates three advanced filtering techniques:
1. Constrained Extended Kalman Filter (CEKF): Enforces state constraints
2. Gaussian Sum Filter (GSF): Models multi-modal posterior distributions
3. Rao-Blackwellized Particle Filter (RBPF): Combines particles with Kalman filters

Scenario: Nonlinear target tracking with constraints on valid state region.
"""

import os
from pathlib import Path

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

from pytcl.dynamic_estimation.gaussian_sum_filter import (
    GaussianComponent,
    GaussianSumFilter,
)
from pytcl.dynamic_estimation.kalman.constrained import (
    ConstrainedEKF,
    ConstraintFunction,
)
from pytcl.dynamic_estimation.rbpf import RBPFFilter

SHOW_PLOTS = os.environ.get("PYTCL_SHOW_PLOTS", "1") != "0"
OUTPUT_DIR = Path("docs/_static/images/examples")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)


class TargetTrackingScenario:
    """Nonlinear target tracking scenario.

    Target moves in 2D with nonlinear dynamics. Measurements are range and
    bearing from a fixed observer.
    """

    def __init__(self, seed: int = 42):
        """Initialize scenario.

        Parameters
        ----------
        seed : int
            Random seed for reproducibility
        """
        np.random.seed(seed)

        # State: [x, y, vx, vy] (position and velocity in Cartesian coords)
        self.state_dim = 4
        self.measurement_dim = 2  # range and bearing

        # System matrices
        self.dt = 0.1
        self.F = np.array(
            [
                [1, 0, self.dt, 0],
                [0, 1, 0, self.dt],
                [0, 0, 1, 0],
                [0, 0, 0, 1],
            ]
        )

        self.Q = np.diag([0, 0, 0.001, 0.001])  # Process noise

        # Measurement observer position
        self.observer = np.array([0.0, 0.0])

        # Measurement noise, as variances: range in units^2, bearing in rad^2.
        # The bearing term used to be 0.01 rad^2, i.e. a 5.7 degree standard
        # deviation. At this target's ~14 unit range that is 1.4 units of
        # cross-range error on a trajectory only 2.5 units long, so the plotted
        # tracks looked like noise even though the filters were behaving. 1e-4
        # is a 0.57 degree sensor, which is still modest but legible.
        self.R = np.diag([0.01, 1e-4])

        # Initial state
        self.x0 = np.array([10.0, 10.0, -1.0, -0.5])
        self.P0 = np.diag([1.0, 1.0, 0.5, 0.5])

    def f(self, x: np.ndarray) -> np.ndarray:
        """Nonlinear state transition with friction.

        Parameters
        ----------
        x : ndarray
            State vector [x, y, vx, vy]

        Returns
        -------
        ndarray
            Next state with velocity friction
        """
        x_next = self.F @ x
        # Add friction to velocity
        x_next[2] *= 0.95
        x_next[3] *= 0.95
        return x_next

    def h(self, x: np.ndarray) -> np.ndarray:
        """Measurement function: range and bearing.

        Parameters
        ----------
        x : ndarray
            State vector [x, y, vx, vy]

        Returns
        -------
        ndarray
            Measurement [range, bearing]
        """
        pos = x[:2]
        delta = pos - self.observer

        # Range
        r = np.linalg.norm(delta)

        # Bearing (angle from East)
        bearing = np.arctan2(delta[1], delta[0])

        return np.array([r, bearing])

    def h_jacobian(self, x: np.ndarray) -> np.ndarray:
        """Jacobian of measurement function.

        Parameters
        ----------
        x : ndarray
            State vector

        Returns
        -------
        ndarray
            Jacobian dh/dx
        """
        pos = x[:2]
        delta = pos - self.observer
        r = np.linalg.norm(delta)

        if r < 0.01:
            # Avoid singularity
            return np.array(
                [
                    [0, 0, 0, 0],
                    [0, 0, 0, 0],
                ]
            )

        H = np.zeros((2, 4))

        # dr/dx = delta[0] / r
        H[0, 0] = delta[0] / r
        H[0, 1] = delta[1] / r

        # dbearing/dx = -delta[1] / r^2, dbearing/dy = delta[0] / r^2
        H[1, 0] = -delta[1] / r**2
        H[1, 1] = delta[0] / r**2

        return H

    def generate_trajectory(self, steps: int = 50):
        """Generate synthetic true trajectory and measurements.

        Parameters
        ----------
        steps : int
            Number of time steps

        Returns
        -------
        x_true : ndarray (steps, 4)
            True state trajectory
        measurements : ndarray (steps, 2)
            Noisy range/bearing measurements
        """
        x_true = np.zeros((steps, 4))
        measurements = np.zeros((steps, 2))

        x_true[0] = self.x0
        measurements[0] = self.h(x_true[0]) + np.random.multivariate_normal(
            np.zeros(2), self.R
        )

        for k in range(1, steps):
            # True dynamics
            x_true[k] = self.f(x_true[k - 1])
            x_true[k] += np.random.multivariate_normal(np.zeros(4), self.Q)

            # Measurement
            z_true = self.h(x_true[k])
            measurements[k] = z_true + np.random.multivariate_normal(
                np.zeros(2), self.R
            )

        return x_true, measurements


def run_cekf_filter(
    scenario: TargetTrackingScenario,
    measurements: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
    """Run Constrained EKF with position constraint.

    Parameters
    ----------
    scenario : TargetTrackingScenario
        Tracking scenario
    measurements : ndarray
        Measurements

    Returns
    -------
    x_est : ndarray
        State estimates
    P_est : ndarray
        Covariance estimates
    """
    cekf = ConstrainedEKF()

    # Operating-area constraint: the target is known to stay inside a circle
    # centered at (5, 5).
    #
    # The radius has to be chosen against the actual trajectory. It was 10.0,
    # but the true track never gets further than 7.07 from the center, so the
    # constraint was inactive at every step and the "constrained" EKF returned
    # exactly what a plain EKF would -- which is why its curve sat invisibly
    # underneath the GSF's.
    #
    # 7.10 sits just above that 7.07 maximum: the true track stays feasible,
    # so the constraint is honest, while estimates that wander outward get
    # pulled back. Tightening it further (6.5, say) would exclude the real
    # target and make the estimate worse, which is the correct behavior for
    # a constraint that is simply wrong.
    CENTER = np.array([5.0, 5.0])
    RADIUS = 7.1

    def g_circle(x):
        # Negative means inside the region.
        return (x[0] - CENTER[0]) ** 2 + (x[1] - CENTER[1]) ** 2 - RADIUS**2

    # Jacobian
    def G_circle(x):
        jac = np.zeros((1, 4))
        jac[0, 0] = 2 * (x[0] - CENTER[0])
        jac[0, 1] = 2 * (x[1] - CENTER[1])
        return jac

    cekf.add_constraint(ConstraintFunction(g_circle, G=G_circle))

    # Initialize
    x = scenario.x0.copy()
    P = scenario.P0.copy()

    x_est = np.zeros((len(measurements), 4))
    P_est = np.zeros((len(measurements), 4, 4))

    for k, z in enumerate(measurements):
        # Predict
        def f_wrapper(x_):
            return scenario.f(x_)

        pred = cekf.predict(x, P, f_wrapper, scenario.F, scenario.Q)
        x = pred.x
        P = pred.P

        # Update
        def h_wrapper(x_):
            return scenario.h(x_)

        upd = cekf.update(x, P, z, h_wrapper, scenario.h_jacobian(x), scenario.R)
        x = upd.x
        P = upd.P

        x_est[k] = x
        P_est[k] = P

    # Report how often the constraint actually did something. A "constrained"
    # filter whose constraint never activates is just an EKF, and saying so is
    # more useful than quietly plotting a third identical curve.
    dist = np.linalg.norm(x_est[:, :2] - CENTER, axis=1)
    n_active = int(np.sum(dist >= RADIUS - 1e-9))
    print(f"  Constraint |x-(5,5)| <= {RADIUS} active on {n_active}/{len(x_est)} steps")

    return x_est, P_est


def run_gsf_filter(
    scenario: TargetTrackingScenario,
    measurements: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
    """Run Gaussian Sum Filter.

    Parameters
    ----------
    scenario : TargetTrackingScenario
        Tracking scenario
    measurements : ndarray
        Measurements

    Returns
    -------
    x_est : ndarray
        State estimates
    P_est : ndarray
        Covariance estimates
    """
    gsf = GaussianSumFilter(max_components=5)

    # Initialize with multiple modes
    gsf.initialize(scenario.x0, scenario.P0, num_components=3)

    x_est = np.zeros((len(measurements), 4))
    P_est = np.zeros((len(measurements), 4, 4))

    for k, z in enumerate(measurements):
        # Predict
        def f_wrapper(x_):
            return scenario.f(x_)

        gsf.predict(f_wrapper, scenario.F, scenario.Q)

        # Update
        def h_wrapper(x_):
            return scenario.h(x_)

        # Get current estimate for Jacobian
        x_pred, _ = gsf.estimate()
        gsf.update(z, h_wrapper, scenario.h_jacobian(x_pred), scenario.R)

        # Estimate
        x, P = gsf.estimate()
        x_est[k] = x
        P_est[k] = P

    return x_est, P_est


def run_rbpf_filter(
    scenario: TargetTrackingScenario,
    measurements: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
    """Run Rao-Blackwellized Particle Filter.

    Parameters
    ----------
    scenario : TargetTrackingScenario
        Tracking scenario
    measurements : ndarray
        Measurements

    Returns
    -------
    x_est : ndarray
        State estimates
    P_est : ndarray
        Covariance estimates
    """
    rbpf = RBPFFilter(max_particles=50)

    # Partition: nonlinear (position), linear (velocity)
    y0 = scenario.x0[:2]  # position
    x0 = scenario.x0[2:]  # velocity
    P0 = scenario.P0[2:, 2:]

    rbpf.initialize(y0, x0, P0, num_particles=30)

    x_est = np.zeros((len(measurements), 4))
    P_est = np.zeros((len(measurements), 4, 4))

    # The RBPF API propagates the nonlinear state as y[k+1] = g(y[k]) + noise,
    # with no velocity-to-position coupling, so position uses a random-walk
    # proposal whose Qy must cover the per-step target motion (~|v0|*dt).
    def g(y):
        return y

    Qy = np.eye(2) * 0.02

    # Linear (Kalman) block: velocity with friction
    F_v = np.eye(2) * 0.95
    Qx = np.eye(2) * 0.001

    def f_linear(v, y):
        return F_v @ v

    def h_rbpf(v, y):
        # Full state from position and velocity
        x_full = np.concatenate([y, v])
        return scenario.h(x_full)

    # Range/bearing depend only on position (the particle state), so the
    # measurement has zero sensitivity to the linear velocity block; the
    # particle weights carry all the position information.
    H_v = np.zeros((2, 2))

    for k, z in enumerate(measurements):
        rbpf.predict(g, Qy, f_linear, F_v, Qx)
        rbpf.update(z, h_rbpf, H_v, scenario.R)

        # Estimate
        y_est, v_est, P_v = rbpf.estimate()
        x_est[k] = np.concatenate([y_est, v_est])

        # Full covariance (approximate)
        P_est[k, :2, :2] = np.eye(2) * 0.1
        P_est[k, 2:, 2:] = P_v
        P_est[k, :2, 2:] = 0
        P_est[k, 2:, :2] = 0

    return x_est, P_est


def plot_filter_comparison(
    x_true: np.ndarray,
    x_cekf: np.ndarray,
    x_gsf: np.ndarray,
    x_rbpf: np.ndarray,
    P_cekf: np.ndarray,
    P_gsf: np.ndarray,
    P_rbpf: np.ndarray,
) -> None:
    """Create interactive Plotly visualizations for filter comparison."""
    # Compute errors
    err_cekf = np.linalg.norm(x_cekf - x_true, axis=1)
    err_gsf = np.linalg.norm(x_gsf - x_true, axis=1)
    err_rbpf = np.linalg.norm(x_rbpf - x_true, axis=1)

    # Compute uncertainties
    unc_cekf = np.array([np.trace(P_cekf[k]) for k in range(len(x_true))])
    unc_gsf = np.array([np.trace(P_gsf[k]) for k in range(len(x_true))])
    unc_rbpf = np.array([np.trace(P_rbpf[k]) for k in range(len(x_true))])

    time = np.arange(len(x_true))

    # Create subplot figure
    fig = make_subplots(
        rows=2,
        cols=2,
        subplot_titles=(
            "Estimated Trajectories",
            "State Estimation Error",
            "Estimated Uncertainty",
            "Error Distribution",
        ),
        specs=[
            [{"type": "scatter"}, {"type": "scatter"}],
            [{"type": "scatter"}, {"type": "box"}],
        ],
    )

    # Plot 1: Trajectories
    fig.add_trace(
        go.Scatter(
            x=x_true[:, 0],
            y=x_true[:, 1],
            mode="lines+markers",
            name="True Trajectory",
            line=dict(color="black", width=3, dash="dash"),
            marker=dict(size=5),
            hovertemplate="<b>True Path</b><br>X: %{x:.2f}<br>Y: %{y:.2f}<extra></extra>",
        ),
        row=1,
        col=1,
    )

    fig.add_trace(
        go.Scatter(
            x=x_cekf[:, 0],
            y=x_cekf[:, 1],
            mode="lines",
            name="CEKF Estimate",
            line=dict(color="blue", width=2),
            hovertemplate="<b>CEKF</b><br>X: %{x:.2f}<br>Y: %{y:.2f}<extra></extra>",
        ),
        row=1,
        col=1,
    )

    fig.add_trace(
        go.Scatter(
            x=x_gsf[:, 0],
            y=x_gsf[:, 1],
            mode="lines",
            name="GSF Estimate",
            line=dict(color="green", width=2),
            hovertemplate="<b>GSF</b><br>X: %{x:.2f}<br>Y: %{y:.2f}<extra></extra>",
        ),
        row=1,
        col=1,
    )

    fig.add_trace(
        go.Scatter(
            x=x_rbpf[:, 0],
            y=x_rbpf[:, 1],
            mode="lines",
            name="RBPF Estimate",
            line=dict(color="red", width=2),
            hovertemplate="<b>RBPF</b><br>X: %{x:.2f}<br>Y: %{y:.2f}<extra></extra>",
        ),
        row=1,
        col=1,
    )

    # Plot 2: Position errors
    fig.add_trace(
        go.Scatter(
            x=time,
            y=err_cekf,
            mode="lines",
            name="CEKF Error",
            line=dict(color="blue", width=2),
            hovertemplate="<b>Time:</b> %{x}<br><b>CEKF Error:</b> %{y:.4f}<extra></extra>",
        ),
        row=1,
        col=2,
    )

    fig.add_trace(
        go.Scatter(
            x=time,
            y=err_gsf,
            mode="lines",
            name="GSF Error",
            line=dict(color="green", width=2),
            hovertemplate="<b>Time:</b> %{x}<br><b>GSF Error:</b> %{y:.4f}<extra></extra>",
        ),
        row=1,
        col=2,
    )

    fig.add_trace(
        go.Scatter(
            x=time,
            y=err_rbpf,
            mode="lines",
            name="RBPF Error",
            line=dict(color="red", width=2),
            hovertemplate="<b>Time:</b> %{x}<br><b>RBPF Error:</b> %{y:.4f}<extra></extra>",
        ),
        row=1,
        col=2,
    )

    # Plot 3: Uncertainty estimates
    fig.add_trace(
        go.Scatter(
            x=time,
            y=unc_cekf,
            mode="lines",
            name="CEKF Uncertainty",
            line=dict(color="blue", width=2),
            hovertemplate="<b>Time:</b> %{x}<br><b>CEKF Covariance Trace:</b> %{y:.4f}<extra></extra>",
        ),
        row=2,
        col=1,
    )

    fig.add_trace(
        go.Scatter(
            x=time,
            y=unc_gsf,
            mode="lines",
            name="GSF Uncertainty",
            line=dict(color="green", width=2),
            hovertemplate="<b>Time:</b> %{x}<br><b>GSF Covariance Trace:</b> %{y:.4f}<extra></extra>",
        ),
        row=2,
        col=1,
    )

    fig.add_trace(
        go.Scatter(
            x=time,
            y=unc_rbpf,
            mode="lines",
            name="RBPF Uncertainty",
            line=dict(color="red", width=2),
            hovertemplate="<b>Time:</b> %{x}<br><b>RBPF Covariance Trace:</b> %{y:.4f}<extra></extra>",
        ),
        row=2,
        col=1,
    )

    # Plot 4: Error distribution (box plot)
    fig.add_trace(
        go.Box(
            y=err_cekf,
            name="CEKF",
            marker_color="blue",
            hovertemplate="<b>CEKF</b><br>Error: %{y:.4f}<extra></extra>",
        ),
        row=2,
        col=2,
    )

    fig.add_trace(
        go.Box(
            y=err_gsf,
            name="GSF",
            marker_color="green",
            hovertemplate="<b>GSF</b><br>Error: %{y:.4f}<extra></extra>",
        ),
        row=2,
        col=2,
    )

    fig.add_trace(
        go.Box(
            y=err_rbpf,
            name="RBPF",
            marker_color="red",
            hovertemplate="<b>RBPF</b><br>Error: %{y:.4f}<extra></extra>",
        ),
        row=2,
        col=2,
    )

    # Update layout
    fig.update_xaxes(title_text="X Position", row=1, col=1)
    fig.update_yaxes(title_text="Y Position", row=1, col=1)

    fig.update_xaxes(title_text="Time Step", row=1, col=2)
    fig.update_yaxes(title_text="Position Error (Norm)", row=1, col=2)

    fig.update_xaxes(title_text="Time Step", row=2, col=1)
    fig.update_yaxes(title_text="Covariance Trace", row=2, col=1)

    fig.update_xaxes(title_text="Filter Algorithm", row=2, col=2)
    fig.update_yaxes(title_text="Position Error", row=2, col=2)

    fig.update_layout(
        title_text="Advanced Filter Comparison: CEKF vs GSF vs RBPF",
        height=900,
        showlegend=True,
        hovermode="closest",
        plot_bgcolor="rgba(240,240,240,0.5)",
    )

    if SHOW_PLOTS:
        fig.show()
    else:
        fig.write_html(
            str(OUTPUT_DIR / "advanced_filters_comparison.html"),
            include_plotlyjs="cdn",
            div_id="advanced_filters_comparison",
        )


def main():
    """Run comparison and generate plots."""
    # Create scenario
    scenario = TargetTrackingScenario()

    # Generate data
    print("Generating synthetic trajectory...")
    x_true, measurements = scenario.generate_trajectory(steps=50)

    # Run filters
    print("Running CEKF...")
    x_cekf, P_cekf = run_cekf_filter(scenario, measurements)

    print("Running GSF...")
    x_gsf, P_gsf = run_gsf_filter(scenario, measurements)

    print("Running RBPF...")
    x_rbpf, P_rbpf = run_rbpf_filter(scenario, measurements)

    # Print statistics
    err_cekf = np.linalg.norm(x_cekf - x_true, axis=1)
    err_gsf = np.linalg.norm(x_gsf - x_true, axis=1)
    err_rbpf = np.linalg.norm(x_rbpf - x_true, axis=1)

    unc_cekf = np.array([np.trace(P_cekf[k]) for k in range(len(x_true))])
    unc_gsf = np.array([np.trace(P_gsf[k]) for k in range(len(x_true))])
    unc_rbpf = np.array([np.trace(P_rbpf[k]) for k in range(len(x_true))])

    print("\n" + "=" * 60)
    print("FILTER COMPARISON RESULTS")
    print("=" * 60)
    print(
        f"CEKF - Mean Error: {np.mean(err_cekf):.4f}, Mean Uncertainty: {np.mean(unc_cekf):.4f}"
    )
    print(
        f"GSF  - Mean Error: {np.mean(err_gsf):.4f}, Mean Uncertainty: {np.mean(unc_gsf):.4f}"
    )
    print(
        f"RBPF - Mean Error: {np.mean(err_rbpf):.4f}, Mean Uncertainty: {np.mean(unc_rbpf):.4f}"
    )
    print("=" * 60)

    # Generate interactive Plotly visualizations
    plot_filter_comparison(x_true, x_cekf, x_gsf, x_rbpf, P_cekf, P_gsf, P_rbpf)


OUTPUT_DIR = Path("docs/_static/images/examples")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

if __name__ == "__main__":
    main()

Running the Example

python examples/advanced_filters_comparison.py

See Also

• Kalman Filter Comparison - Basic Kalman filter variants

• Particle Filters - Standard particle filters

• Gaussian Mixtures and Clustering - Gaussian mixture operations