State Estimation Analysis

The State-Estimation Analysis performs model-based state estimation on time-series data using Kalman-type filters. Given a dynamical system model and noisy measurements, the analysis computes state estimates and their uncertainties, enabling reconstruction of the internal system state (including individual state variables) that cannot be directly measured.

This analysis leverages LowLevelParticleFiltersMTK.jl, which bridges ModelingToolkit.jl models with state estimation algorithms from LowLevelParticleFilters.jl.

State Estimation Theory

State-Space Model with Noise

The analysis assumes the system follows a discrete-time state-space model with process and measurement noise:

\begin{aligned} x_{k + 1} & = f (x_{k}, u_{k}, p, t_{k}, w_{k}) \\ y_{k} & = g (x_{k}, u_{k}, p, t_{k}) + e_{k} \end{aligned}

where:

Variable $x_{k}$ is the state vector at time $k$
Variable $u_{k}$ is the control input
Variable $y_{k}$ is the measured output
Variable $w_{k} \sim N (0, R_{1})$ is process noise (disturbance)
Variable $e_{k} \sim N (0, R_{2})$ is measurement noise
Functions $f$ and $g$ are the system dynamics and measurement functions

Disturbance Modeling

Unlike most textbook state-estimation procedures where process noise is added directly to state variables, this analysis requires explicit disturbance inputs in the model. Process noise $w_{k}$ enters through designated disturbance input variables, allowing the framework to handle the fact that the state realization is chosen by ModelingToolkit and not by the user.

Filtering vs Smoothing

Filtering computes the state estimate ${\hat{x}}_{k}$ using measurements up to time $k$ :

{\hat{x}}_{k} = E [x_{k} ∣ y_{1}, \dots, y_{k}]

Smoothing computes the state estimate using all available measurements:

{\hat{x}}_{k} = E [x_{k} ∣ y_{1}, \dots, y_{N}]

Smoothing generally provides better estimates when processing data offline, as it incorporates both past and future measurements. The smooth function performs the backward smoothing pass.

Filter Types

Extended Kalman Filter (EKF)

The ExtendedKalmanFilter linearizes the nonlinear dynamics and measurement functions at each time step using first-order Taylor expansion:

\begin{aligned} A_{k} & = \frac{\partial f}{\partial x} |_{{\hat{x}}_{k}} \\ C_{k} & = \frac{\partial g}{\partial x} |_{{\hat{x}}_{k}} \end{aligned}

Advantages:

Computationally efficient
Works well for mildly nonlinear systems

Limitations:

May diverge for highly nonlinear systems
Requires differentiable dynamics

Unscented Kalman Filter (UKF)

The UnscentedKalmanFilter uses a deterministic sampling technique (sigma points) to capture the mean and covariance of the state distribution through nonlinear transformations without requiring Jacobians.

Advantages:

Better accuracy for highly nonlinear systems
No need to compute Jacobians

Limitations:

Slightly more computationally expensive than EKF
May require tuning of UKF-specific parameters (sigma-point weights)

Example Definition

The following example demonstrates state estimation using a simple dynamical system with measured input-output data. We will create a component with explicit disturbance inputs and then perform state estimation using the Extended Kalman Filter.

Adding Required Dependencies

This example requires the following Julia packages:

julia

using Pkg
Pkg.add(["DyadControlSystems", "DyadData", "DyadInterface"])

Dyad Code

First, let's define a simple second-order system component in Dyad. If you already have a component library, add this to a new .dyad file in the dyad folder.

dyad

component SimpleDynamicalSystem
  # State variables (internal dynamics)
  variable x1::Real
  variable x2::Real

  # Control input
  u = RealInput()

  # Measured output
  y = RealOutput()

  # Disturbance inputs (where process noise enters)
  w1 = RealInput()
  w2 = RealInput()

  # Parameters
  parameter damping::Real = 0.5

relations
  initial x1 = 1.0
  initial x2 = 0.0
  der(x1) = -damping*x1 + x2 + w1
  der(x2) = -x1 - damping*x2 + u + w2
  y = x1
end

Before we can estimate the state, we need measurement data. In practice, this would come from real sensors. For this example, we'll generate synthetic data by first simulating the system using a TransientAnalysis.

dyad

using BlockComponents

# Create a system with a sine wave input for simulation
component SimpleDynamicalSystemWithInput
  sys = SimpleDynamicalSystem()
  input_signal = BlockComponents.Sine(frequency=1.0, amplitude=0.5, phase=0)
relations
  connect(input_signal.y, sys.u)
  # Connect disturbances to zero for nominal simulation
  sys.w1 = 0
  sys.w2 = 0
end

# Simulate to generate "ground truth" data
analysis GenerateData
  extends TransientAnalysis(start=0, stop=10.0, saveat=0.1)
  model = SimpleDynamicalSystemWithInput()
end

Running the Analysis

NOTE

To run the analyses defined above, you need to use the Julia REPL and load the package that corresponds to the component library with using <LibraryName>

First, let's simulate the system to generate measurement data:

julia

using DyadInterface
using Random
Random.seed!(1234)

# Generate ground truth data
sim_result = GenerateData()

Now we extract the data and add measurement noise to simulate real sensor readings:

julia

# Extract the simulation solution
sol_raw = artifacts(sim_result, :RawSolution)
t_data = sol_raw.t
u_data = sol_raw[sim_result.spec.model.input_signal.y]' |> Matrix  # Transpose to get shape [nu, nT]
y_true = sol_raw[sim_result.spec.model.sys.y]

# Add measurement noise to simulate real sensors
measurement_noise_std = 0.1
y_data = (y_true .+ measurement_noise_std * randn(length(y_true)))' |> Matrix  # Shape [ny, nT]

# Sample time (inferred from data)
Ts = t_data[2] - t_data[1]

println("Generated $(length(t_data)) data points with Ts = $Ts")

Generated 101 data points with Ts = 0.1

Now we can create and run a state estimation analysis using the generated noisy measurements. Since we're using runtime-generated data, we'll create the analysis specification directly in Julia:

julia

using DyadControlSystems

# Create state estimation specification with the generated data
spec = StateEstimationAnalysisSpec(
    name = :StateEst,
    model = SimpleDynamicalSystem(name=:simple),
    input_data = u_data,
    output_data = y_data,
    inputs = ["u"],
    outputs = ["y"],
    disturbance_inputs = ["w1", "w2"],
    estimator = "ExtendedKalmanFilter",
    r1_diag = [0.01, 0.01],      # Process noise covariance
    r2_diag = [0.1^2],              # Measurement noise covariance (matches added noise)
    filtering_mode = "filtering",
    discretization = "Rk4",
    Ts = Ts,
    plot_confidence = true,
    confidence_level = 1.96       # 95% confidence bounds
)

# Run state estimation
est_result = run_analysis(spec)

julia

# Plot estimated state trajectories with confidence bounds
artifacts(est_result, :state)

The state plot shows the estimated state variables (x1 and x2) over time with shaded confidence intervals representing the uncertainty in the estimates.

julia

# Plot predicted outputs vs measured outputs
artifacts(est_result, :output)

This plot compares the filter's predicted measurements (one-step-ahead predictions) against the actual measured output data, showing how well the filter tracks the measurements despite the noise.

julia

# Analyze filter innovations (residuals)
artifacts(est_result, :innovation_analysis)

The innovation analysis shows the prediction errors (innovations) and their autocorrelation. For a well-tuned filter, innovations should appear as white noise with no significant autocorrelation.

We can also compare the estimated state with the true state from the simulation:

julia

using Plots

# Get true state trajectories from simulation using sol.spec.model
x1_true = sol_raw[sim_result.spec.model.sys.x1]
x2_true = sol_raw[sim_result.spec.model.sys.x2]

# Get estimated state trajectories using symbolic indexing on the prob.iosys
x1_est = est_result.sol[est_result.prob.iosys.x1]
x2_est = est_result.sol[est_result.prob.iosys.x2]

# Plot comparison
p1 = plot(t_data, x1_true, label="True x1", linewidth=2)
plot!(p1, t_data, x1_est, label="Estimated x1", linestyle=:dash, linewidth=2)
xlabel!(p1, "Time (s)")
ylabel!(p1, "x1")

p2 = plot(t_data, x2_true, label="True x2", linewidth=2)
plot!(p2, t_data, x2_est, label="Estimated x2", linestyle=:dash, linewidth=2)
xlabel!(p2, "Time (s)")
ylabel!(p2, "x2")

plot(p1, p2, layout=(2,1), size=(800, 600))

This comparison shows how well the Extended Kalman Filter recovers the true state from noisy measurements. The estimated state (dashed lines) closely tracks the true state (solid lines) despite only having access to the noisy output measurement y.

Analysis Arguments

The following arguments define a State Estimation Analysis:

Required Arguments

model: The ModelingToolkit system model (must not be structurally simplified).
inputs::Vector{String}: Names of control inputs to the system (can be empty for autonomous systems).
outputs::Vector{String}: Names of measured outputs available from sensors.
disturbance_inputs::Vector{String}: Names of disturbance inputs where process noise enters the dynamics.

Data Input Options

You must provide data using one of these approaches:

Option 1: Direct Arrays

input_data::Real[nu, nT]: Input data matrix where rows are input signals and columns are time points (can be empty for autonomous systems).
output_data::Real[ny, nT]: Output data matrix where rows are output signals and columns are time points.

Option 2: DyadDataset

dataset: A DyadDataset object containing the time series data.
input_cols::String[nu]: Names of input signal columns in the dataset (can be empty).
output_cols::String[ny]: Names of output signal columns in the dataset.

Estimator Configuration

estimator::String = "ExtendedKalmanFilter": Filter type
- "ExtendedKalmanFilter": EKF, linearization-based (faster, good for mildly nonlinear systems)
- "UnscentedKalmanFilter": UKF, sigma-point based (more accurate for highly nonlinear systems)
r1_diag::Vector{Real}: Diagonal elements of process noise covariance matrix R₁. Length must match number of disturbance inputs.
r2_diag::Vector{Real}: Diagonal elements of measurement noise covariance matrix R₂. Length must match number of measured outputs.
filtering_mode::String = "filtering": Estimation mode
- "filtering": Forward pass only, estimates based on past measurements
- "smoothing": Forward-backward pass, estimates using all measurements (better accuracy for offline data)

Discretization Parameters

discretization::String = "Rk4": Discretization method for continuous-time models
- "Rk4": Fourth-order Runge-Kutta (recommended for pure ODEs)
- "Trapezoidal": Trapezoidal rule (supports algebraic equations)
- "SimpleColloc": Simple collocation method
Ts::Real: Sample time of the data. Set to -1 to automatically infer from dataset.

Initial Conditions

sigma0::Real = 1e-4: Initial state standard deviation for uncertainty initialization.

Visualization Parameters

plot_confidence::Boolean = true: Whether to plot confidence bounds on state estimates.
confidence_level::Real = 1.96: Number of standard deviations for confidence bounds (1.96 ≈ 95% confidence).
n_samples::Integer = 0: Number of Monte Carlo samples for uncertainty visualization (0 for deterministic mean).

State Estimation Algorithm

The analysis follows this systematic workflow:

Model Preparation: The ModelingToolkit model is compiled and disturbance inputs are identified
Data Processing: Input-output data is extracted from arrays or DyadDataset
Discretization: Continuous-time dynamics are discretized using the specified method
Filter Construction: EKF or UKF is instantiated with noise covariance matrices
Forward Filtering: The filter processes measurements sequentially using forward_trajectory, computing ${\hat{x}}_{k}$ and covariance $P_{k}$
Smoothing (Optional): If smoothing mode is selected, a backward pass using smooth refines estimates using future measurements
Solution Packaging: Results are wrapped with symbolic indexing capabilities

Design Considerations

Noise Covariance Tuning

Process Noise R₁ (r1_diag):

Higher values: Less trust in model, more reliance on measurements
Lower values: More trust in model, smoother estimates
Should reflect actual disturbance magnitudes and model uncertainty

Measurement Noise R₂ (r2_diag):

Higher values: Less trust in sensors, rely more on model predictions
Lower values: More trust in sensors, track measurements closely
Should match actual sensor noise characteristics

Practical Tuning Guidelines

Start with measurement noise: Estimate R₂ from sensor datasheets or by analyzing measurement variability when the system is at rest
Tune process noise: Adjust R₁ until the filter innovations (residuals) appear as white noise
Check innovation statistics: Use the innovation analysis artifact to verify filter consistency
Balance tracking vs smoothness: Higher R₁ → faster tracking, noisier estimates; Lower R₁ → smoother estimates, slower tracking

For detailed troubleshooting advice and additional tuning strategies, see the LowLevelParticleFilters.jl Troubleshooting and Tuning Guide.

Choosing Between EKF and UKF

Use EKF when:

System is mildly nonlinear
Computational efficiency is critical
Dynamics are smooth and well-behaved

Use UKF when:

System is highly nonlinear
Higher accuracy is needed
You want to avoid computing Jacobians

Filtering vs Smoothing

Use filtering when:

Processing data online (real-time)
Implementing in an embedded system
Causal estimates are required

Use smoothing when:

Processing data offline (batch)
Maximum accuracy is needed
Computational cost is not a concern

Artifacts

A StateEstimationAnalysis returns the following artifacts:

Plots

:state: Filtered or smoothed state trajectories with optional confidence bounds
:output: Comparison of predicted outputs vs measured outputs
:innovation_analysis: Innovation (residual) sequence and autocorrelation for filter consistency checking
:native_filter_plot: Comprehensive visualization including predictions, filtered estimates, covariances, inputs, and outputs

Tables

:filter_metrics: DataFrame containing:
- RMS innovation (prediction error) for each output
- Log-likelihood of the data given the model
- Number of state variables, inputs, and outputs
- Sample time and number of time points

Custom Visualization

The analysis supports symbolic indexing for custom visualization:

julia

using DyadInterface

# Create custom visualization specification
vizdef = PlotlyVisualizationSpec(
    [:x1, :x2, :y],  # Variables to plot
    (; plot_confidence = true, confidence_level = 1.96),
    Attribute[]
)

# Generate custom plot
customizable_visualization(sol, vizdef)

Notes

The model must include disturbance input variables marked with [disturbance = true, input = true] to specify where process noise enters
For autonomous systems (no control inputs), set inputs = String[] and include all uncertain inputs as disturbance inputs
The state realization is determined automatically by ModelingToolkit during compilation
Sample time must be uniform; non-uniformly sampled data is not currently supported
For systems with algebraic equations (DAEs), use discretization = "Trapezoidal" instead of "Rk4". See the DAE State Estimation Example for a complete demonstration

Scalar Types

Partial Types

Scalar Types

Function Types

Partial Types

Partial Types

Partial Types

Partial Types

Scalar Types

Partial Types

Partial Types

State Estimation Analysis ​

State Estimation Theory ​

State-Space Model with Noise ​

Disturbance Modeling ​

Filtering vs Smoothing ​

Filter Types ​

Extended Kalman Filter (EKF) ​

Unscented Kalman Filter (UKF) ​

Example Definition ​

Adding Required Dependencies ​

Dyad Code ​

Running the Analysis ​

Analysis Arguments ​

Required Arguments ​

Data Input Options ​

Option 1: Direct Arrays ​

Option 2: DyadDataset ​

Estimator Configuration ​

Discretization Parameters ​

Initial Conditions ​

Visualization Parameters ​

State Estimation Algorithm ​

Design Considerations ​

Noise Covariance Tuning ​

Practical Tuning Guidelines ​

Choosing Between EKF and UKF ​

Filtering vs Smoothing ​

Artifacts ​

Plots ​

Tables ​

Custom Visualization ​

Notes ​

Further Reading ​

Documentation ​

Advanced Examples ​

Background Theory ​

State Estimation Analysis

State Estimation Theory

State-Space Model with Noise

Disturbance Modeling

Filtering vs Smoothing

Filter Types

Extended Kalman Filter (EKF)

Unscented Kalman Filter (UKF)

Example Definition

Adding Required Dependencies

Dyad Code

Running the Analysis

Analysis Arguments

Required Arguments

Data Input Options

Option 1: Direct Arrays

Option 2: DyadDataset

Estimator Configuration

Discretization Parameters

Initial Conditions

Visualization Parameters

State Estimation Algorithm

Design Considerations

Noise Covariance Tuning

Practical Tuning Guidelines

Choosing Between EKF and UKF

Filtering vs Smoothing

Artifacts

Plots

Tables

Custom Visualization

Notes

Further Reading

Documentation

Advanced Examples

Background Theory