gds-analysis: step and impulse response computation with time-domain metrics

[rororowyourboat]

Summary

Add response.py to gds_analysis/ providing step response, impulse response computation, and time-domain performance metrics (rise time, settling time, overshoot, steady-state error). Behind the existing [continuous] extra.

Parent issue: #198 (classical control theory stack)

Motivation

Step response metrics are the first thing a controls engineer checks — they are the "unit tests" of control design. GDS can already simulate systems via gds-sim and gds-continuous, producing time-series results. But there is no post-processing to extract the standard performance metrics. This is ~100 lines of code with immediate value for every example in gds-examples/control/.

Proposed API

New file: gds_analysis/response.py (behind [continuous] extra)

Response Generation

def step_response(
    A: list[list[float]],
    B: list[list[float]],
    C: list[list[float]],
    D: list[list[float]],
    t_span: tuple[float, float] = (0.0, 10.0),
    n_points: int = 1000,
    input_index: int = 0,
) -> tuple[list[float], list[list[float]]]:
    """Compute step response from state-space model.
    
    Applies unit step to input_index, simulates via scipy.signal.step.
    Returns (times, outputs) where outputs[i] is the response of output i.
    """

def impulse_response(
    A: list[list[float]],
    B: list[list[float]],
    C: list[list[float]],
    D: list[list[float]],
    t_span: tuple[float, float] = (0.0, 10.0),
    n_points: int = 1000,
    input_index: int = 0,
) -> tuple[list[float], list[list[float]]]:
    """Compute impulse response. Wraps scipy.signal.impulse."""

Performance Metrics

@dataclass(frozen=True)
class StepMetrics:
    """Standard time-domain performance metrics for a step response."""
    rise_time: float          # time from 10% to 90% of final value (seconds)
    settling_time: float      # time to stay within ±2% of final value (seconds)
    overshoot_pct: float      # peak overshoot as percentage of final value
    peak_time: float          # time of first peak (seconds)
    peak_value: float         # value at first peak
    steady_state_value: float # final value (average of last 5% of samples)
    steady_state_error: float # |setpoint - steady_state_value|

def step_response_metrics(
    times: list[float],
    values: list[float],
    setpoint: float = 1.0,
    settling_band: float = 0.02,
    rise_lower: float = 0.1,
    rise_upper: float = 0.9,
) -> StepMetrics:
    """Extract standard performance metrics from a step response.
    
    Works on any time-series data — can be fed from:
    - step_response() output (this module)
    - ODEResults.state_array() (gds-continuous)
    - gds-sim Results column extraction
    
    Args:
        times: time values (monotonically increasing)
        values: response values
        setpoint: desired final value (default 1.0 for unit step)
        settling_band: fractional band for settling time (default ±2%)
        rise_lower/rise_upper: fractional thresholds for rise time (default 10%/90%)
    """

Convenience for GDS Results

def metrics_from_ode_results(
    results: "ODEResults",
    state_name: str,
    setpoint: float = 1.0,
    **kwargs,
) -> StepMetrics:
    """Extract step metrics from an ODEResults object.
    
    Convenience wrapper: calls results.times and results.state_array(state_name),
    then delegates to step_response_metrics().
    """

Implementation Notes

Metric definitions (matching standard control textbooks)

Metric	Definition
Rise time	Time for response to go from rise_lower × final_value to rise_upper × final_value. Linear interpolation between samples.
Settling time	Earliest time after which the response stays within ±settling_band × final_value permanently. Scan backward from end.
Overshoot %	(peak - final_value) / final_value × 100. Zero if no overshoot.
Peak time	Time of maximum value.
Steady-state value	Mean of last 5% of samples (robust to minor oscillation).
Steady-state error	`

Settling time — backward scan algorithm

Rather than scanning forward (which misclassifies early transient crossings), scan backward from the end to find the last sample outside the settling band:

final = steady_state_value
band = settling_band * abs(final)
for i in range(len(values) - 1, -1, -1):
    if abs(values[i] - final) > band:
        settling_time = times[i]
        break

Interface contract

• step_response() / impulse_response() accept list[list[float]] matrices, return list[float]
• step_response_metrics() accepts list[float] — works with any data source
• metrics_from_ode_results() is the only function that imports from gds-continuous; guarded by try/except ImportError

Relation to existing gds_analysis/metrics.py

The existing trajectory_distances() computes user-defined distance functions between successive states. response.py is complementary — it computes standard control-specific metrics on single time-series. These could live in the same file, but a separate file is cleaner given the different domains (generic trajectory metrics vs. control step response).

Key Files

• packages/gds-analysis/gds_analysis/metrics.py — existing metrics infrastructure
• packages/gds-continuous/gds_continuous/results.py — ODEResults type (optional import)
• packages/gds-examples/control/thermostat/ — immediate consumer for testing
• packages/gds-examples/control/double_integrator/ — immediate consumer for testing

Testing

• test_response.py:
  • First-order 1/(s+1): rise time ≈ 2.2s, settling time ≈ 4.0s (±2%), overshoot = 0%
  • Second-order underdamped (ζ=0.5, ω_n=1): overshoot ≈ 16.3%, verify against analytical formula exp(-πζ/√(1-ζ²)) × 100
  • Critically damped (ζ=1): overshoot = 0%, rise time > underdamped
  • Pure integrator step response: no steady state → settling_time = inf or last sample time
  • Pre-computed data: feed known list[float] directly to step_response_metrics() — no scipy needed for this test path

Concepts Addressed (MATLAB Tech Talks)

• Video 3: Step response — rise time, overshoot, settling time, steady-state error, second-order system parameters (ω_n, ζ)