archimedes.sysid.pem¶

archimedes.sysid.pem(predictor: KalmanFilterBase, data: Timeseries | tuple[Timeseries, ...], p_guess: T, x0: np.ndarray | tuple[np.ndarray, ...], estimate_x0: bool | tuple[bool] = False, bounds: tuple[T, T] | None = None, P0: np.ndarray | tuple[np.ndarray] = None, method: str = 'lm', options: dict | None = None) → LMResult¶

Estimate parameters using Prediction Error Minimization.

Solves the system identification problem by minimizing the prediction error between model predictions and measured outputs using a Kalman filter framework. This approach provides optimal handling of process and measurement noise while enabling efficient gradient computation through automatic differentiation.

The method implements the discrete-time prediction error objective:

minimize  J = (1/N) Σ[k=1 to N] e[k]ᵀ e[k]

where e[k] = y[k] - ŷ[k|k-1] are the one-step-ahead prediction errors (innovations) from the Kalman filter and N is the number of measurements.

This formulation automatically accounts for:

Noise handling: Process and measurement noise are modeled explicitly
Recursive estimation: Kalman filter provides efficient state propagation
Gradient computation: Automatic differentiation through filter recursions

Parameters:

predictor (KalmanFilterBase) –
Kalman filter implementing the system model. Must provide step(t, x, y, P, args) method. Common choices:
- ExtendedKalmanFilter
- UnscentedKalmanFilter
The predictor encapsulates the system dynamics, observation model, and noise characteristics (Q, R matrices).
data (Timeseries) –
Input-output data containing synchronized time series:
- ts : Time vector of shape (N,)
- us : Input signals of shape (nu, N)
- ys : Output measurements of shape (ny, N)
All arrays must have consistent time dimensions. Multiple Timeseries instances can be provided as a tuple for multi-experiment identification, allowing joint parameter estimation across different datasets.

If multiple datasets are provided, the initial state estimates (x0) should also be a tuple of initial conditions corresponding to each dataset.
p_guess (PyTree) – Initial parameter guess with arbitrary nested structure (e.g., {"mass": 1.0, "damping": {"c1": 0.1, "c2": 0.2}}). The optimization preserves this structure in the result, enabling natural organization of physical parameters.
x0 (array_like, optional) – Initial state estimate of shape (nx,). Used as an initial guess if estimate_x0=True. For multiple datasets, this can be a tuple of initial conditions, allowing different initial states for each dataset.
estimate_x0 (bool, default=False) – Whether to estimate the initial state x0 along with parameters.
bounds (tuple of (PyTree, PyTree), optional) –
Parameter bounds as (lower_bounds, upper_bounds) with the same PyTree structure as p_guess. Enables physical constraints such as:
- Positive masses, stiffnesses, damping coefficients
- Bounded gain parameters, time constants
- Realistic physical parameter ranges
Use -np.inf and np.inf for unbounded parameters.
P0 (array_like, optional) – Initial state covariance matrix of shape (nx, nx). If None, defaults to identity matrix. Represents uncertainty in initial state estimate. Can be a tuple of matrices if multiple datasets are provided, allowing different initial uncertainties for each dataset.
method (str, default="bfgs") – Optimization method. Currently only “lm” (Levenberg-Marquardt), “ipopt”, and “bfgs” (BFGS) are supported. The “lm” method is a custom implementation roughly based on the MINPACK algorithm, and the “bfgs” method dispatches to a SciPy wrapper (“BFGS” or “L-BFGS-B”, depending on whether there are bounds).
options (dict, optional) –
Optimization options passed to the underlying optimization solver.

For the “lm” method, these include:
- ftol : Function tolerance (default: 1e-4)
- xtol : Parameter tolerance (default: 1e-6)
- gtol : Gradient tolerance (default: 1e-6)
- max_nfev : Maximum function evaluations (default: 200)
- nprint : Progress printing interval (default: 0)
For the “bfgs” method, these include:
- gtol : Gradient tolerance (default: 1e-6)
- disp : Whether to print convergence information (default: True)
- maxiter : Maximum iterations (default: 200)
- hess_inv0 : Initial Hessian inverse for BFGS (optional)
If hess_inv0 is not provided, a Gauss-Newton-like Hessian approximation is used to initialize the BFGS inverse-Hessian approximation.

For the “ipopt” method, see the archimedes.minimize() documentation for available options. The solver defaults to a limited-memory approximation of the Hessian.

Returns:

result – Optimization result with estimated parameters in result.x preserving the original PyTree structure. Additional fields include:

success : Whether estimation succeeded
fun : Final prediction error objective value
nit : Number of optimization iterations
history : Detailed convergence history

Return type:

scipy.optimize.OptimizeResult

Notes

This implementation provides:

Automatic Gradients:: Efficient gradient computation through automatic differentiation of the Kalman filter recursions. No need for finite differences or manual Jacobian implementation.
Structured Parameters:: Natural handling of nested parameter dictionaries enables intuitive organization of physical parameters and parameter bounds.
Physical Constraints:: Box constraints enable realistic parameter bounds (mass > 0, etc.) without sacrificing convergence properties.
Kalman Filter Integration:: Seamless integration with both Extended and Unscented Kalman Filters enables handling of linear and nonlinear systems with appropriate accuracy-efficiency tradeoffs.
Multi-Experiment Support:: Can handle multiple datasets simultaneously, allowing joint parameter estimation across different experiments or operating conditions.

The method automatically computes gradients with respect to both initial conditions (when estimate_x0=True) and model parameters using efficient automatic differentiation through the Kalman filter recursions. This avoids the computational expense and numerical issues of finite difference approximations.

Examples

>>> import numpy as np
>>> import archimedes as arc
>>> from archimedes.sysid import pem, Timeseries
>>> from archimedes.observers import ExtendedKalmanFilter
>>>
>>> # Define second-order damped oscillator
>>> def dynamics(t, x, u, params):
...     omega_n = params["omega_n"]  # Natural frequency
...     zeta = params["zeta"]        # Damping ratio
...
...     return np.hstack([
...         x[1],  # velocity
...         -omega_n**2 * x[0] - 2*zeta*omega_n*x[1] + omega_n**2 * u[0]
...     ])
>>>
>>> def observation(t, x, u, params):
...     return x[0]  # Measure position only
>>>
>>> # Generate synthetic measurement data
>>> dt = 0.05
>>> ts = np.arange(0, 10, dt)
>>> x0 = np.array([0.0, 0.0])
>>> params_true = {"omega_n": 2.0, "zeta": 0.1}
>>> us = np.ones((1, len(ts)))  # Step input
>>> # Generate step response
>>> xs = arc.odeint(
...     dynamics,
...     (ts[0], ts[-1]),
...     x0,
...     t_eval=ts,
...     args=(np.array([1.0]), params_true),
... )
>>> noise_std = 0.01
>>> ys = xs[:1, :] + np.random.normal(0, noise_std, size=xs.shape[1])
>>>
>>> # Set up identification problem
>>> dyn_discrete = arc.discretize(dynamics, dt, method="rk4")
>>> Q = noise_std ** 2 * np.eye(2)  # Process noise covariance
>>> R = noise_std ** 2 * np.eye(1)  # Measurement noise covariance
>>> ekf = ExtendedKalmanFilter(dyn_discrete, observation, Q, R)
>>>
>>> data = Timeseries(ts=ts, us=us, ys=ys)
>>> p_guess = {"omega_n": 2.5, "zeta": 0.5}
>>>
>>> # Estimate parameters with known initial conditions
>>> result = pem(ekf, data, p_guess, x0=x0)
>>> print(f"Estimated parameters: {result.p}")
Estimated parameters: {'omega_n': array(1.9709515), 'zeta': array(0.11517324)}
>>> print(f"Converged in {result.nit} iterations")
Converged in 26 iterations
>>>
>>> # With physical parameter constraints
>>> bounds = (
...     {"omega_n": 0.0, "zeta": 0.0},     # Lower bounds (positive values)
...     {"omega_n": 10.0, "zeta": 1.0},    # Upper bounds (reasonable ranges)
... )
>>> result = pem(ekf, data, p_guess, x0=x0, bounds=bounds)