Simulation Metrics

Control System Performance Metrics

Overview

This skill provides a small, self-contained Python module (simulation_metrics.py) that computes the four classic time-domain performance metrics for a control system's step response: rise time, overshoot percentage, steady-state error, and settling time. These metrics turn a wiggly response curve into a small set of comparable scalars, enabling objective, repeatable evaluation against design specifications such as "≤10% overshoot, settle inside the 2% band in under 0.5 s."

Two deliberate engineering choices shape the implementation:

• NumPy for the math. Vectorized threshold detection (np.flatnonzero) and stable tail averaging (np.mean) stay fast and numerically well-behaved even on long, fine-grained traces, where a hand-rolled Python loop would be slow and easy to get subtly wrong.
• Explicit type hints and up-front validation. The annotations make the contract obvious to callers and to static checkers (no untyped/Any parameters), and the shared validator catches malformed data before any indexing happens. Every metric assumes times[i] lines up with values[i]; a length mismatch or a non-monotonic time vector would otherwise yield a plausible but incorrect answer instead of an error.

When to Use

Reach for this skill when you need to put hard numbers on a control system's step response — quantitative, time-domain figures you can check against a design spec. Rise time, overshoot, steady-state error, and settling time answer, respectively: how fast does it react, how far does it overshoot, where does it finally land, and when is it effectively done?

Trigger keywords: rise time, overshoot, settling time, steady-state error, step response, PID tuning, control system performance, simulation metrics, time-domain analysis, 2% band, 10-90 rise.

In practice you want these whenever you are:

1. Tuning a PID loop and need quantitative feedback on each iteration.
2. Comparing two controller designs head-to-head on the same step input.
3. Proving that a simulation meets a written requirement (e.g., "≤10% overshoot, settle inside the 2% band in under 0.5 s").

When not to use (and why it matters)

• These are time-domain metrics, not frequency-domain ones. They are computed directly from a response-versus-time curve, so they say nothing about Bode/Nyquist behavior or gain/phase stability margins. A system can post excellent time-domain numbers and still be fragile to phase lag or modeling error — if margins are what you care about, use a frequency-domain tool instead.
• They assume a reasonably clean signal. Every function reads the trajectory you hand it literally, noise and all. On a high-bandwidth, noisy trace a single spike can register as overshoot or kick the settling time far to the right, producing pessimistic-but-wrong results. Low-pass filter or smooth the data first so the metrics describe the underlying response rather than the measurement noise.
• A zero setpoint is rejected on purpose. Overshoot and the tolerance band are normalized by the target, so target == 0 would divide by zero. Rather than return a misleading number, the functions raise ValueError. The same philosophy applies to empty, length-mismatched, or non-finite arrays: a silently-wrong metric buried in a tuning report is far more dangerous than a loud, early failure.

Prerequisites

• Python 3.10+ (uses from __future__ import annotations and X | Y union syntax).
• NumPy installed: pip install numpy
• NumPy type stubs (optional but recommended for static checking): pip install numpy-types
• On Windows (PowerShell), ensure Python is on your PATH:

  python --version
  python -c "import numpy; print(numpy.__version__)"
  

Procedure

Step 1 — Create the module file

Create simulation_metrics.py in your project. The module contains shared validation helpers and four metric functions. All definitions below belong in this single file.

Step 2 — Add shared imports and validation helpers

These definitions are the foundation for every metric that follows; keep them at the top of the same module so the functions below can reuse them. The helpers exist so the why of each guard lives in one place rather than being copy-pasted (and drifting) across four functions.

from __future__ import annotations

from collections.abc import Sequence

import numpy as np
import numpy.typing as npt

# A concrete, dense float array is what every metric operates on internally.
FloatArray = npt.NDArray[np.float64]
# Callers may pass plain Python sequences or numpy arrays; both are accepted.
SignalInput = Sequence[float] | FloatArray


def _as_validated_signal(
    times: SignalInput,
    values: SignalInput,
) -> tuple[FloatArray, FloatArray]:
    """Coerce raw inputs into aligned float64 arrays and reject malformed data.

    Why validate up front: every metric indexes by position and assumes that
    ``times[i]`` corresponds to ``values[i]``. A length mismatch, a stray NaN, or a
    non-monotonic time vector would silently produce a plausible-but-wrong number,
    which is worse in a tuning report than a hard failure. Strictly increasing time
    is required because the interpolation steps below need distinct, ordered x-values.
    """
    t = np.asarray(times, dtype=np.float64)
    v = np.asarray(values, dtype=np.float64)

    if t.ndim != 1 or v.ndim != 1:
        raise ValueError("times and values must be one-dimensional sequences.")
    if t.size == 0 or v.size == 0:
        raise ValueError("times and values must be non-empty.")
    if t.size != v.size:
        raise ValueError(
            f"times and values must be the same length (got {t.size} and {v.size})."
        )
    if not np.all(np.isfinite(t)) or not np.all(np.isfinite(v)):
        raise ValueError("times and values must not contain NaN or infinity.")
    if np.any(np.diff(t) <= 0.0):
        raise ValueError("times must be strictly increasing.")
    return t, v


def _validate_target(target: float) -> float:
    """Validate the setpoint, guarding the division-by-zero that a 0 target causes.

    Overshoot and the settling tolerance band are both normalized by the target, so a
    zero setpoint is undefined for these metrics rather than merely inconvenient.
    """
    target = float(target)
    if not np.isfinite(target):
        raise ValueError("target must be a finite number.")
    if target == 0.0:
        raise ValueError(
            "target must be non-zero; metrics are normalized by the target and a "
            "zero setpoint would divide by zero."
        )
    return target


def _first_crossing_time(
    times: FloatArray,
    normalized: FloatArray,
    level: float,
) -> float | None:
    """Return the linearly interpolated time at which ``normalized`` first reaches ``level``.

    Interpolating between the two bracketing samples removes the discretization bias
    you get from snapping to the nearest sample time, which matters when the solver
    step is coarse relative to how quickly the response moves. Returns ``None`` when the
    level is never reached within the captured window.
    """
    reached = np.flatnonzero(normalized >= level)
    if reached.size == 0:
        return None

    idx = int(reached[0])
    if idx == 0:
        # Already at/above the level at the first sample; there is no earlier point
        # to interpolate against, so report the first timestamp as-is.
        return float(times[0])

    n_prev, n_curr = float(normalized[idx - 1]), float(normalized[idx])
    t_prev, t_curr = float(times[idx - 1]), float(times[idx])
    if n_curr == n_prev:
        return t_curr

    fraction = (level - n_prev) / (n_curr - n_prev)
    return t_prev + fraction * (t_curr - t_prev)

Step 3 — Add rise_time

Rise time is the interval taken to climb from 10% to 90% of the target. The 10–90 band is the convention rather than 0–100% for a concrete reason: a step response approaches its final value asymptotically, so "time to reach 100%" is often never well-defined, while the first instants near 0% are dominated by sensor dead-zone. The 10–90 window is the standard, reproducible way to capture how fast the system reacts. Returning None when the 90% level is never reached is itself a useful diagnostic — it means the system is too slow or unstable to qualify within the simulated window.

def rise_time(
    times: SignalInput,
    values: SignalInput,
    target: float,
    low: float = 0.10,
    high: float = 0.90,
) -> float | None:
    """Time to climb from ``low`` (10%) to ``high`` (90%) of the target.

    The thresholds are exposed as parameters but default to the standard 10/90 band so
    callers can switch to, say, 5/95 without rewriting the function. The crossings are
    interpolated between samples for sub-step accuracy, and the calculation is sign-aware
    via the normalized signal, so it works for negative setpoints too. Returns ``None``
    if the response never reaches the high threshold in the captured window.
    """
    t, v = _as_validated_signal(times, values)
    target = _validate_target(target)
    if not (0.0 <= low < high <= 1.0):
        raise ValueError(f"Require 0 <= low < high <= 1 (got low={low}, high={high}).")

    # Dividing by the target maps the response onto a 0 -> 1 progress curve regardless
    # of the target's sign, so the same >= comparison detects crossings either way.
    normalized = v / target
    t_low = _first_crossing_time(t, normalized, low)
    t_high = _first_crossing_time(t, normalized, high)
    if t_low is None or t_high is None:
        return None
    return t_high - t_low

Step 4 — Add overshoot_percent

Overshoot is the peak excursion beyond the target, expressed as a percentage of the target. Normalizing by the target is what makes the figure meaningful across setpoints: a 2-unit overshoot is trivial for a setpoint of 500 but enormous for a setpoint of 5. The percentage form is the dimensionless number engineers compare against a spec like "≤10% overshoot." It matters physically because excessive overshoot can saturate or damage actuators, trip safety limits, or signal an under-damped, ringing loop.

def overshoot_percent(
    values: SignalInput,
    target: float,
) -> float:
    """Peak excursion beyond the target, as a percentage of the target.

    Sign-aware: for a negative setpoint the worst-case excursion is the *minimum*, not
    the maximum, so the peak is taken in the direction of the target's sign. Returns
    ``0.0`` for critically- or over-damped responses that never pass the target, which
    is the correct answer rather than a small negative artifact.
    """
    v = np.asarray(values, dtype=np.float64)
    if v.size == 0:
        raise ValueError("values must be non-empty.")
    if not np.all(np.isfinite(v)):
        raise ValueError("values must not contain NaN or infinity.")
    target = _validate_target(target)

    # The most extreme sample in the direction the system was commanded to move.
    peak = float(np.max(v)) if target > 0 else float(np.min(v))
    # Positive only when the peak actually exceeds the target; clamp so a response
    # that never reaches the target reports 0.0 overshoot rather than a negative value.
    excess_fraction = (peak - target) / target
    return max(0.0, excess_fraction * 100.0)

Step 5 — Add steady_state_error

Steady-state error is the absolute residual gap between the target and the value the system finally settles at. The function averages the last final_fraction of the record rather than reading the single last sample, and the reason is robustness: the final sample can still carry ripple or measurement noise, whereas averaging the tail gives a stable estimate of where the response truly settled. A non-zero result on a loop with integral action usually points to actuator saturation or an un-modeled disturbance — i.e. it is diagnostic, not just a pass/fail number.

def steady_state_error(
    values: SignalInput,
    target: float,
    final_fraction: float = 0.10,
) -> float:
    """Absolute residual error between the target and the settled output.

    Averaging the final ``final_fraction`` of the record smooths out late ripple so the
    estimate reflects the steady value rather than a single noisy endpoint. The window is
    clamped to contain at least one sample, so short runs do not produce an empty slice.
    """
    v = np.asarray(values, dtype=np.float64)
    if v.size == 0:
        raise ValueError("values must be non-empty.")
    if not np.all(np.isfinite(v)):
        raise ValueError("values must not contain NaN or infinity.")
    target = _validate_target(target)
    if not (0.0 < final_fraction <= 1.0):
        raise ValueError(f"final_fraction must be in (0, 1] (got {final_fraction}).")

    n = v.size
    # Guarantee at least one sample in the averaging window even when final_fraction
    # rounds the window size below one on a short run.
    start = min(n - 1, int(n * (1.0 - final_fraction)))
    final_avg = float(np.mean(v[start:]))
    return abs(target - final_avg)

Step 6 — Add settling_time

Settling time is the earliest time after which the response stays inside a ±tolerance band around the target for the rest of the record. The subtle part is why we key off the last exit from the band rather than the first entry: an under-damped system can dip into the band, overshoot back out, and re-enter. The metric is defined as the instant after which it never leaves again, so the code locates the last out-of-band sample and reports the crossing into the band that follows it. The tolerance is a fraction of the target magnitude (0.02 → the classic 2% band, 0.05 → 5%), because "close enough" is naturally expressed relative to the commanded value.

def settling_time(
    times: SignalInput,
    values: SignalInput,
    target: float,
    tolerance: float = 0.02,
) -> float | None:
    """Earliest time after which the response stays within ±tolerance of the target.

    Keying off the last out-of-band sample correctly handles a response that re-enters
    the band after an overshoot. The exact band-crossing instant is interpolated between
    the last out-of-band sample and the first in-band sample for sub-step accuracy.
    Returns ``None`` when the response is still outside the band at the end of the window
    (it has not settled within the captured time), and ``times[0]`` when it is already
    inside the band from the very first sample.
    """
    t, v = _as_validated_signal(times, values)
    target = _validate_target(target)
    if not (0.0 < tolerance < 1.0):
        raise ValueError(f"tolerance must be in (0, 1) (got {tolerance}).")

    band = abs(target) * tolerance
    outside = np.flatnonzero(np.abs(v - target) > band)

    if outside.size == 0:
        return float(t[0])  # Within tolerance from the very first sample.

    last_outside = int(outside[-1])
    settled_idx = last_outside + 1
    if settled_idx >= t.size:
        return None  # Never re-entered the band before the run ended.

    # Interpolate the exact moment the trajectory re-enters the band, crossing whichever
    # boundary (upper or lower) it last violated.
    v_out = float(v[last_outside])
    v_in = float(v[settled_idx])
    boundary = (target + band) if v_out > target else (target - band)
    denom = v_in - v_out
    if denom == 0.0:
        return float(t[settled_idx])

    fraction = (boundary - v_out) / denom
    fraction = min(1.0, max(0.0, fraction))  # Stay within the bracketing interval.
    return float(t[last_outside]) + fraction * (float(t[settled_idx]) - float(t[last_outside]))

Step 7 — Call the functions on your simulation output

Import the module and call each function on paired times and values arrays plus the commanded target setpoint. Always handle None return values from rise_time and settling_time as first-class outcomes, not crashes.

Examples

Worked example: second-order-like step response

from __future__ import annotations

# Each row is one simulation timestep captured from the solver.
results: list[dict[str, float]] = [
    {"time": 0.0, "value": 0.0},
    {"time": 0.1, "value": 5.0},
    {"time": 0.2, "value": 15.0},
    {"time": 0.3, "value": 32.0},   # transient peak -> overshoot
    {"time": 0.4, "value": 30.5},
    {"time": 0.5, "value": 30.1},
    {"time": 0.6, "value": 30.0},
]

times: list[float] = [row["time"] for row in results]
values: list[float] = [row["value"] for row in results]
target: float = 30.0

rt = rise_time(times, values, target)
st = settling_time(times, values, target)

print(f"Rise time:     {rt:.4f} s" if rt is not None else "Rise time:     never reached 90%")
print(f"Overshoot:     {overshoot_percent(values, target):.2f} %")
print(f"SS error:      {steady_state_error(values, target):.4f}")
print(f"Settling time: {st:.4f} s" if st is not None else "Settling time: never settled")

# Expected output:
# Rise time:     0.2106 s
# Overshoot:     6.67 %
# SS error:      0.0000
# Settling time: 0.3933 s

Running on Windows (PowerShell)

python simulation_metrics.py

Or run a quick inline check:

python -c "from simulation_metrics import rise_time, overshoot_percent, steady_state_error, settling_time; print('all imports OK')"

Pitfalls

1. Noisy signals produce wrong metrics. A single noise spike can register as overshoot or push the settling time far right. Always low-pass filter or smooth the data before computing metrics so they describe the underlying response, not measurement noise.

2. target == 0 raises ValueError by design. Overshoot and the tolerance band are normalized by the target. Rather than returning a misleading number, the functions fail loudly. If your setpoint is zero, you need a different metric definition (e.g., absolute tolerance band instead of relative).

3. Non-monotonic time vectors are rejected. The interpolation steps require strictly increasing, distinct x-values. Duplicate or backwards timestamps raise ValueError immediately — this is intentional, as interpolation on non-monotonic data would produce silently wrong results.

4. settling_time keys off the last out-of-band sample, not the first entry. An under-damped system can dip into the band, overshoot back out, and re-enter. The metric is defined as the instant after which it never leaves again. If you mistakenly look for the first entry, you will report a settling time that is too optimistic.

5. steady_state_error on a short run still works. The averaging window is clamped to contain at least one sample even when final_fraction rounds the window below one. However, on very short runs the "steady-state" estimate is necessarily coarse — use a longer simulation window when possible.

6. overshoot_percent returns 0.0 for over-damped responses. This is the correct answer, not a bug. A response that never reaches the target has zero overshoot, not a negative overshoot.

7. Negative setpoints are handled, but the peak direction flips. For target < 0, the worst-case excursion is the minimum of the signal, not the maximum. The code handles this automatically via the sign check, but be aware when interpreting results.

8. Plain Python lists and numpy arrays produce identical results. The _as_validated_signal helper coerces both to float64 arrays. Do not worry about which type you pass — but be aware that very large lists will consume memory during conversion.

Verification

Run these checks after implementing the module to confirm correctness:

• Import check — all four functions import without error:

  python -c "from simulation_metrics import rise_time, overshoot_percent, steady_state_error, settling_time; print('OK')"
  

• Worked example — run the worked example above and confirm output matches:

  Rise time:     0.2106 s
  Overshoot:     6.67 %
  SS error:      0.0000
  Settling time: 0.3933 s
  

• rise_time returns None when the 90% level is never reached, and interpolates the 10%/90% crossings between samples rather than snapping to the nearest sample time.

• overshoot_percent returns 0.0 for over- or critically-damped responses and is sign-correct for negative setpoints (it measures the minimum excursion when target < 0).

• settling_time keys off the last out-of-band sample (so re-entry after an overshoot is handled) and interpolates the band-crossing instant; it returns None when the response is still outside the band at the end of the window.

• steady_state_error averages at least one sample even when final_fraction rounds the window below one sample on a short run.

• Input validation rejects empty, length-mismatched, non-finite, or non-monotonic-time data with a clear ValueError instead of returning a plausible-but-wrong number:

  python -c "from simulation_metrics import rise_time; rise_time([0,1,2], [1,2], 10)"
  # Should raise: ValueError: times and values must be the same length
  

• target == 0 raises (the division-by-zero guard):

  python -c "from simulation_metrics import overshoot_percent; overshoot_percent([1,2,3], 0)"
  # Should raise: ValueError: target must be non-zero...
  

• Inputs supplied as plain Python lists and as numpy arrays produce identical results.

• Validate against a known second-order step response: compare the computed metrics to the closed-form values for a chosen damping ratio (ζ) and natural frequency (ωn).

Related skills

• system-stability-analysis
• pid-tuning-optimization
• signal-processing-filters
• time-series-analysis