"""Open-loop step-response identification of FOPDT/SOPDT models. The workflow this module supports: 1. The user performs (or simulates, via :func:`pidtune.simulation.loop.simulate_open_loop`) an open-loop step test: with the loop in manual, the input is stepped from a steady baseline and the output is recorded until it settles. 2. The recording is wrapped in :class:`StepResponseData`, which validates shape, monotonic time, and finiteness *at construction* — garbage data fails loudly here, not deep inside a fit. 3. :func:`fit_fopdt` or :func:`fit_sopdt` estimates model parameters using the method selected from :class:`FitMethod`, returning a :class:`FitResult` that carries the model **and** fit diagnostics (RMSE, R²) so users can judge whether rule-based tuning on that model is justified. Method survey, formulas and applicability limits are specified in ``docs/design/03_tuning_methods.md``, section "Step-response identification" — in brief: * ``TANGENT`` — classic Ziegler–Nichols tangent-at-inflection construction; simple, but noise-sensitive (the inflection slope is a derivative estimate). * ``TWO_POINT`` — Smith's 28.3 % / 63.2 % two-point method; robust to moderate noise, FOPDT only. * ``AREA`` — Åström–Hägglund area method; integrates rather than differentiates, so it is the most noise-tolerant closed-form option. * ``LEAST_SQUARES`` — full nonlinear least squares on the simulated response; most accurate, required for SOPDT, needs an initial guess (seeded internally from ``TWO_POINT``). Fitting functions are **API stubs** in this milestone; data records and validation are final. """ from __future__ import annotations import dataclasses import enum from typing import Optional, TYPE_CHECKING import numpy as np from numpy.typing import NDArray if TYPE_CHECKING: # annotation-only imports from ..models.base import ProcessModel __all__ = ["StepResponseData", "FitMethod", "FitResult", "fit_fopdt", "fit_sopdt"] class FitMethod(enum.Enum): """Parameter-estimation method for step-response fitting. Attributes ---------- TANGENT Tangent at the inflection point (Ziegler–Nichols process reaction curve). FOPDT only. Fast; degrades badly with noise. TWO_POINT Smith's two-point method using the 28.3 % and 63.2 % rise times. FOPDT only. The recommended default for clean-to- moderately-noisy data. AREA Åström–Hägglund area method (average residence time from the integral of the normalised response). FOPDT only. Best closed-form choice for noisy data. LEAST_SQUARES Nonlinear least squares minimising the output-error sum of squares between the recorded response and the model's simulated response. Works for FOPDT and SOPDT; the only method accepted by :func:`fit_sopdt`. """ TANGENT = "tangent" TWO_POINT = "two_point" AREA = "area" LEAST_SQUARES = "least_squares" @dataclasses.dataclass(frozen=True) class StepResponseData: """A validated open-loop step-test recording. Validation happens in ``__post_init__`` so that every fitting routine can assume well-formed input. Arrays are coerced to contiguous ``float64`` copies, making instances safe to share. Attributes ---------- t : numpy.ndarray 1-D sample times. Must be strictly increasing and contain at least 10 samples (fewer cannot constrain even an FOPDT fit). Need *not* be uniformly spaced — field data rarely is; methods that require a uniform grid resample internally. y : numpy.ndarray 1-D recorded output, same length as ``t``. u_step : float Magnitude of the applied input step (signed). Must be non-zero — the process gain is ``Δy / u_step``. step_time : float, optional Time at which the input stepped. Default ``0.0``. Must lie within ``[t[0], t[-1])``. y_baseline : float, optional Steady output level before the step. When ``None`` (default), fitters estimate it as the mean of pre-step samples; pass it explicitly if the recording starts at or after ``step_time``. Raises ------ ValueError If any validation rule above is violated. Messages name the offending field and constraint, per the conventions in ``docs/design/05_errors_extensibility.md``. Examples -------- >>> import numpy as np >>> t = np.linspace(0.0, 50.0, 501) >>> y = 2.0 * (1 - np.exp(-np.clip(t - 3.0, 0.0, None) / 10.0)) >>> data = StepResponseData(t=t, y=y, u_step=1.0) >>> len(data.t) 501 """ t: NDArray[np.float64] y: NDArray[np.float64] u_step: float step_time: float = 0.0 y_baseline: Optional[float] = None def __post_init__(self) -> None: t = np.ascontiguousarray(self.t, dtype=np.float64) y = np.ascontiguousarray(self.y, dtype=np.float64) object.__setattr__(self, "t", t) object.__setattr__(self, "y", y) if t.ndim != 1 or y.ndim != 1: raise ValueError( f"StepResponseData.t and .y must be 1-D arrays, got shapes " f"{t.shape} and {y.shape}." ) if t.shape[0] != y.shape[0]: raise ValueError( f"StepResponseData.t and .y must have equal length, got " f"{t.shape[0]} and {y.shape[0]}." ) if t.shape[0] < 10: raise ValueError( f"StepResponseData requires at least 10 samples to constrain a fit, " f"got {t.shape[0]}." ) if not np.all(np.isfinite(t)) or not np.all(np.isfinite(y)): raise ValueError( "StepResponseData.t and .y must be finite (no NaN/inf); clean the " "recording before constructing StepResponseData." ) if not np.all(np.diff(t) > 0.0): raise ValueError( "StepResponseData.t must be strictly increasing; sort and " "de-duplicate timestamps before constructing StepResponseData." ) if not np.isfinite(self.u_step) or self.u_step == 0.0: raise ValueError( f"StepResponseData.u_step must be a finite non-zero number, got " f"{self.u_step!r} (a zero input step carries no gain information)." ) if not (t[0] <= self.step_time < t[-1]): raise ValueError( f"StepResponseData.step_time ({self.step_time}) must lie within " f"[t[0], t[-1]) = [{t[0]}, {t[-1]})." ) if self.y_baseline is not None and not np.isfinite(self.y_baseline): raise ValueError( f"StepResponseData.y_baseline must be finite or None, got " f"{self.y_baseline!r}." ) @dataclasses.dataclass(frozen=True) class FitResult: """Outcome of a step-response fit: the model plus its diagnostics. Returned (never mutated) by :func:`fit_fopdt` and :func:`fit_sopdt`. Carrying diagnostics alongside the model is a deliberate API decision: rule-based tuning on a poorly fitting model is the most common way to get a bad PID, and the library should make fit quality impossible to overlook. Attributes ---------- model : pidtune.models.base.ProcessModel The identified model (FOPDT or SOPDT instance). method : FitMethod The method that produced it. rmse : float Root-mean-square error between the recorded output and the model's simulated response on the recording's time grid, in output units. r_squared : float Coefficient of determination of the same comparison. Values below ≈ 0.95 indicate the chosen structure is questionable; fitters emit a ``UserWarning`` in that case (they do not raise — judging adequacy is the user's call). residuals : numpy.ndarray or None Pointwise ``y_data - y_model`` when requested via ``return_residuals=True``; otherwise ``None`` to keep results lightweight. """ model: "ProcessModel" method: FitMethod rmse: float r_squared: float residuals: Optional[NDArray[np.float64]] = None def fit_fopdt( data: StepResponseData, method: FitMethod = FitMethod.TWO_POINT, return_residuals: bool = False, ) -> FitResult: """Fit a first-order-plus-dead-time model to a step recording. Estimates gain :math:`K`, time constant :math:`\\tau` and dead time :math:`\\theta` of :class:`pidtune.models.fopdt.FOPDTModel` using the selected :class:`FitMethod`. Regardless of method, the returned :class:`FitResult` always reports diagnostics computed from a full simulated-response comparison, so methods are directly comparable. Parameters ---------- data : StepResponseData Validated step-test recording. method : FitMethod, optional Estimation method. Default :attr:`FitMethod.TWO_POINT`. ``LEAST_SQUARES`` seeds itself from a ``TWO_POINT`` estimate. return_residuals : bool, optional Include pointwise residuals in the result. Default ``False``. Returns ------- FitResult With ``result.model`` an ``FOPDTModel``. Raises ------ pidtune.exceptions.IdentificationError If the response does not permit the construction the method needs — e.g. the output never crosses the 63.2 % level (``TWO_POINT``), no inflection point is found (``TANGENT``), the response has not settled within the recording (``AREA``, detected via a final-segment slope test), or the optimiser fails to converge (``LEAST_SQUARES``). Messages state which construction failed and what to record instead, per ``docs/design/05_errors_extensibility.md``. Warns ----- UserWarning If the final ``r_squared`` is below 0.95, indicating FOPDT may be a poor structure for this process (consider :func:`fit_sopdt`). Examples -------- >>> result = fit_fopdt(data, method=FitMethod.AREA) # doctest: +SKIP >>> result.model # doctest: +SKIP FOPDTModel(gain=2.01, time_constant=9.8, dead_time=3.1) """ raise NotImplementedError( "fit_fopdt is an API stub; implementation is scheduled for the " "identification milestone (docs/design/06_roadmap.md)." ) def fit_sopdt( data: StepResponseData, method: FitMethod = FitMethod.LEAST_SQUARES, return_residuals: bool = False, ) -> FitResult: """Fit a second-order-plus-dead-time model to a step recording. Estimates the parameters of :class:`pidtune.models.sopdt.SOPDTModel` (gain, two time constants — or natural frequency and damping for the underdamped form — and dead time). Only :attr:`FitMethod.LEAST_SQUARES` is accepted: the closed-form constructions in the other methods are FOPDT-specific. The optimiser is seeded from a ``TWO_POINT`` FOPDT fit with the second time constant initialised to one tenth of the first (overdamped branch) or from a log-decrement estimate when the recording overshoots (underdamped branch). Parameters ---------- data : StepResponseData Validated step-test recording. method : FitMethod, optional Must be :attr:`FitMethod.LEAST_SQUARES`; any other value raises ``ValueError`` immediately. return_residuals : bool, optional Include pointwise residuals in the result. Default ``False``. Returns ------- FitResult With ``result.model`` an ``SOPDTModel``. Raises ------ ValueError If ``method`` is not ``LEAST_SQUARES``. pidtune.exceptions.IdentificationError If the optimiser fails to converge or converges to a boundary of the feasible parameter region (which signals structural mismatch rather than a usable model). Warns ----- UserWarning If ``r_squared`` is below 0.95, or if the identified second time constant is small enough (< 5 % of the first) that an FOPDT model would serve equally well — over-parameterised models tune no better and identify less reliably. """ raise NotImplementedError( "fit_sopdt is an API stub; implementation is scheduled for the " "identification milestone (docs/design/06_roadmap.md)." )