Time-Domain Simulation & Metrics¶

time_domain ¶

Time-domain simulation functions for Linear Time-Invariant (LTI) systems.

step_response ¶

step_response(sys, T=None, X0=None, input_index=0)

Compute the step response of a continuous-time LTI system.

Simulates the system output trajectory for a unit step excitation:

.. math::

u(t) = 1(t) = \begin{cases} 1, & t \ge 0 \\ 0, & t < 0 \end{cases}

Parameters:

Name	Type	Description	Default
`sys`	`LinearTimeInvariant`	The LTI system (TransferFunction or StateSpace).	required
`T`	`Sequence[float] \| NDArray[floating] \| float \| None`	Time vector or duration. If None, automatically determined from dominant poles.	`None`
`X0`	`Sequence[float] \| NDArray[floating] \| None`	Initial state vector (only valid for StateSpace systems).	`None`
`input_index`	`int`	Input channel index for multi-input systems, default is 0.	`0`

Returns:

Type	Description
`TimeResponseData`	Object containing time vector `t`, output trajectory `y`, and state trajectory `x`.

Raises:

Type	Description
`ValueError`	If initial condition X0 is passed to a TransferFunction model or if X0 dimensions do not match state-space model.
`IndexError`	If input_index is out of bounds.
`TypeError`	If sys is not a LinearTimeInvariant instance.

Source code in src/ctrlpy/time_domain.py

def step_response(
    sys: LinearTimeInvariant,
    T: Sequence[float] | NDArray[np.floating] | float | None = None,
    X0: Sequence[float] | NDArray[np.floating] | None = None,
    input_index: int = 0,
) -> TimeResponseData:
    r"""Compute the step response of a continuous-time LTI system.

    Simulates the system output trajectory for a unit step excitation:

    .. math::

        u(t) = 1(t) = \begin{cases} 1, & t \ge 0 \\ 0, & t < 0 \end{cases}

    Parameters
    ----------
    sys : LinearTimeInvariant
        The LTI system (TransferFunction or StateSpace).
    T : Sequence[float] | NDArray[np.floating] | float | None, optional
        Time vector or duration. If None, automatically determined from dominant poles.
    X0 : Sequence[float] | NDArray[np.floating] | None, optional
        Initial state vector (only valid for StateSpace systems).
    input_index : int, optional
        Input channel index for multi-input systems, default is 0.

    Returns
    -------
    TimeResponseData
        Object containing time vector `t`, output trajectory `y`, and state trajectory `x`.

    Raises
    ------
    ValueError
        If initial condition X0 is passed to a TransferFunction model or if
        X0 dimensions do not match state-space model.
    IndexError
        If input_index is out of bounds.
    TypeError
        If sys is not a LinearTimeInvariant instance.
    """
    if not isinstance(sys, LinearTimeInvariant):
        raise TypeError(f"Expected LinearTimeInvariant instance, got {type(sys).__name__}.")

    t_arr = _generate_time_vector(sys, T)

    if isinstance(sys, TransferFunction):
        if X0 is not None:
            raise ValueError(
                "Initial state X0 cannot be specified for TransferFunction models. "
                "Convert to StateSpace first."
            )
        t_out, y_out = signal.step((sys.num, sys.den), T=t_arr)
        return TimeResponseData(t=t_out, y=y_out, x=None, sys=sys, poles=sys.poles())

    if isinstance(sys, StateSpace):
        if not (0 <= input_index < sys.inputs):
            raise IndexError(f"input_index {input_index} out of range for {sys.inputs} inputs.")
        x0_arr = (
            np.zeros(sys.n_states, dtype=np.float64)
            if X0 is None
            else np.asarray(X0, dtype=np.float64).ravel()
        )
        if x0_arr.size != sys.n_states:
            raise ValueError(
                f"Initial state X0 must have length {sys.n_states}, got {x0_arr.size}."
            )

        if sys.inputs == 1:
            u: NDArray[np.float64] = np.ones(len(t_arr), dtype=np.float64)
        else:
            u = np.zeros((len(t_arr), sys.inputs), dtype=np.float64)
            u[:, input_index] = 1.0

        t_out, y_out, x_out = signal.lsim(
            (sys.A, sys.B, sys.C, sys.D),
            u,
            t_arr,
            X0=x0_arr,
        )
        if sys.outputs == 1 and y_out.ndim == 2:
            y_out = y_out.ravel()
        if x_out.ndim == 1:
            x_out = x_out.reshape(-1, sys.n_states)

        return TimeResponseData(t=t_out, y=y_out, x=x_out, sys=sys, poles=sys.poles())

    raise TypeError(f"Unsupported system type: {type(sys).__name__}")

impulse_response ¶

impulse_response(sys, T=None, X0=None, input_index=0)

Compute the impulse response of a continuous-time LTI system.

Simulates the system output trajectory for a Dirac delta excitation:

.. math::

u(t) = \delta(t) \implies y(t) = g(t) = \mathcal{L}^{-1}\{G(s)\}

Parameters:

Name	Type	Description	Default
`sys`	`LinearTimeInvariant`	The LTI system (TransferFunction or StateSpace).	required
`T`	`Sequence[float] \| NDArray[floating] \| float \| None`	Time vector or duration. If None, automatically determined from dominant poles.	`None`
`X0`	`Sequence[float] \| NDArray[floating] \| None`	Initial state vector before the impulse (only valid for StateSpace systems).	`None`
`input_index`	`int`	Input channel index for multi-input systems, default is 0.	`0`

Returns:

Type	Description
`TimeResponseData`	Object containing time vector `t`, output trajectory `y`, and state trajectory `x`.

Raises:

Type	Description
`ValueError`	If initial condition X0 is passed to a TransferFunction model or if X0 dimensions do not match state-space model.
`IndexError`	If input_index is out of bounds.
`TypeError`	If sys is not a LinearTimeInvariant instance.

Source code in src/ctrlpy/time_domain.py

def impulse_response(
    sys: LinearTimeInvariant,
    T: Sequence[float] | NDArray[np.floating] | float | None = None,
    X0: Sequence[float] | NDArray[np.floating] | None = None,
    input_index: int = 0,
) -> TimeResponseData:
    r"""Compute the impulse response of a continuous-time LTI system.

    Simulates the system output trajectory for a Dirac delta excitation:

    .. math::

        u(t) = \delta(t) \implies y(t) = g(t) = \mathcal{L}^{-1}\{G(s)\}

    Parameters
    ----------
    sys : LinearTimeInvariant
        The LTI system (TransferFunction or StateSpace).
    T : Sequence[float] | NDArray[np.floating] | float | None, optional
        Time vector or duration. If None, automatically determined from dominant poles.
    X0 : Sequence[float] | NDArray[np.floating] | None, optional
        Initial state vector before the impulse (only valid for StateSpace systems).
    input_index : int, optional
        Input channel index for multi-input systems, default is 0.

    Returns
    -------
    TimeResponseData
        Object containing time vector `t`, output trajectory `y`, and state trajectory `x`.

    Raises
    ------
    ValueError
        If initial condition X0 is passed to a TransferFunction model or if
        X0 dimensions do not match state-space model.
    IndexError
        If input_index is out of bounds.
    TypeError
        If sys is not a LinearTimeInvariant instance.
    """
    if not isinstance(sys, LinearTimeInvariant):
        raise TypeError(f"Expected LinearTimeInvariant instance, got {type(sys).__name__}.")

    t_arr = _generate_time_vector(sys, T)

    if isinstance(sys, TransferFunction):
        if X0 is not None:
            raise ValueError(
                "Initial state X0 cannot be specified for TransferFunction models. "
                "Convert to StateSpace first."
            )
        t_out, y_out = signal.impulse((sys.num, sys.den), T=t_arr)
        return TimeResponseData(t=t_out, y=y_out, x=None, sys=sys, poles=sys.poles())

    if isinstance(sys, StateSpace):
        if not (0 <= input_index < sys.inputs):
            raise IndexError(f"input_index {input_index} out of range for {sys.inputs} inputs.")
        x0_arr = (
            np.zeros(sys.n_states, dtype=np.float64)
            if X0 is None
            else np.asarray(X0, dtype=np.float64).ravel()
        )
        if x0_arr.size != sys.n_states:
            raise ValueError(
                f"Initial state X0 must have length {sys.n_states}, got {x0_arr.size}."
            )

        b_col = sys.B[:, input_index]
        x_init = b_col + x0_arr
        b_zero = np.zeros_like(sys.B)

        t_out, y_out, x_out = signal.lsim(
            (sys.A, b_zero, sys.C, sys.D),
            0.0,
            t_arr,
            X0=x_init,
            interp=False,
        )
        if sys.outputs == 1 and y_out.ndim == 2:
            y_out = y_out.ravel()
        if x_out.ndim == 1:
            x_out = x_out.reshape(-1, sys.n_states)

        return TimeResponseData(t=t_out, y=y_out, x=x_out, sys=sys, poles=sys.poles())

    raise TypeError(f"Unsupported system type: {type(sys).__name__}")

forced_response ¶

forced_response(sys, T, U, X0=None)

Compute the simulation response of an LTI system to arbitrary inputs \(u(t)\).

Parameters:

Name	Type	Description	Default
`sys`	`LinearTimeInvariant`	The LTI system (TransferFunction or StateSpace).	required
`T`	`Sequence[float] \| NDArray[floating]`	Simulation time points (1D array).	required
`U`	`Sequence[float] \| NDArray[floating] \| float \| int`	Input signal \(u(t)\) over time. For SISO systems, a 1D array of length `len(T)` or a scalar. For MIMO systems, a 2D array of shape `(len(T), inputs)`.	required
`X0`	`Sequence[float] \| NDArray[floating] \| None`	Initial state vector \(x(0)\) (only valid for StateSpace systems).	`None`

Returns:

Type	Description
`TimeResponseData`	Object containing time vector `t`, output trajectory `y`, and state trajectory `x`.

Raises:

Type	Description
`ValueError`	If U or X0 dimensions are incompatible with T or sys.
`TypeError`	If sys is not a LinearTimeInvariant instance.

Source code in src/ctrlpy/time_domain.py

def forced_response(
    sys: LinearTimeInvariant,
    T: Sequence[float] | NDArray[np.floating],
    U: Sequence[float] | NDArray[np.floating] | float | np.number[Any],
    X0: Sequence[float] | NDArray[np.floating] | None = None,
) -> TimeResponseData:
    r"""Compute the simulation response of an LTI system to arbitrary inputs $u(t)$.

    Parameters
    ----------
    sys : LinearTimeInvariant
        The LTI system (TransferFunction or StateSpace).
    T : Sequence[float] | NDArray[np.floating]
        Simulation time points (1D array).
    U : Sequence[float] | NDArray[np.floating] | float | int
        Input signal $u(t)$ over time. For SISO systems, a 1D array of length `len(T)` or a scalar.
        For MIMO systems, a 2D array of shape `(len(T), inputs)`.
    X0 : Sequence[float] | NDArray[np.floating] | None, optional
        Initial state vector $x(0)$ (only valid for StateSpace systems).


    Returns
    -------
    TimeResponseData
        Object containing time vector `t`, output trajectory `y`, and state trajectory `x`.

    Raises
    ------
    ValueError
        If U or X0 dimensions are incompatible with T or sys.
    TypeError
        If sys is not a LinearTimeInvariant instance.
    """
    if not isinstance(sys, LinearTimeInvariant):
        raise TypeError(f"Expected LinearTimeInvariant instance, got {type(sys).__name__}.")

    t_arr = _generate_time_vector(sys, T)
    n_t = len(t_arr)

    # Process and validate input U
    if isinstance(U, (int, float, np.number)):
        if sys.inputs == 1:
            u_arr: NDArray[np.float64] = np.full(n_t, float(U), dtype=np.float64)
        else:
            u_arr = np.full((n_t, sys.inputs), float(U), dtype=np.float64)
    else:
        u_raw = np.asarray(U, dtype=np.float64)
        if sys.inputs == 1:
            u_arr = u_raw.ravel()
            if u_arr.size != n_t:
                raise ValueError(
                    f"Input U length {u_arr.size} must match time vector T length {n_t}."
                )
        else:
            if u_raw.ndim == 1:
                if u_raw.size == sys.inputs:
                    u_arr = np.tile(u_raw, (n_t, 1))
                else:
                    raise ValueError(
                        f"1D input U of size {u_raw.size} does not match {sys.inputs} inputs."
                    )
            elif u_raw.ndim == 2:
                if u_raw.shape == (n_t, sys.inputs):
                    u_arr = u_raw
                elif u_raw.shape == (sys.inputs, n_t) and n_t != sys.inputs:
                    u_arr = u_raw.T
                else:
                    raise ValueError(
                        f"Input U shape {u_raw.shape} is incompatible with T of "
                        f"length {n_t} and {sys.inputs} inputs."
                    )
            else:
                raise ValueError(f"Input U array must be 1D or 2D, got {u_raw.ndim}D.")

    if isinstance(sys, TransferFunction):
        if X0 is not None:
            raise ValueError(
                "Initial state X0 cannot be specified for TransferFunction models. "
                "Convert to StateSpace first."
            )
        t_out, y_out, _ = signal.lsim((sys.num, sys.den), u_arr, t_arr)
        return TimeResponseData(t=t_out, y=y_out, x=None, sys=sys, poles=sys.poles())

    if isinstance(sys, StateSpace):
        x0_arr = (
            np.zeros(sys.n_states, dtype=np.float64)
            if X0 is None
            else np.asarray(X0, dtype=np.float64).ravel()
        )
        if x0_arr.size != sys.n_states:
            raise ValueError(
                f"Initial state X0 must have length {sys.n_states}, got {x0_arr.size}."
            )

        t_out, y_out, x_out = signal.lsim(
            (sys.A, sys.B, sys.C, sys.D),
            u_arr,
            t_arr,
            X0=x0_arr,
        )
        if sys.outputs == 1 and y_out.ndim == 2:
            y_out = y_out.ravel()
        if x_out.ndim == 1:
            x_out = x_out.reshape(-1, sys.n_states)

        return TimeResponseData(t=t_out, y=y_out, x=x_out, sys=sys, poles=sys.poles())

    raise TypeError(f"Unsupported system type: {type(sys).__name__}")

simulation_results ¶

Simulation results data structures and performance metrics extraction.

TimeResponseData `dataclass` ¶

Container for time-domain simulation response data and performance metrics.

Parameters:

Name	Type	Description	Default
`t`	`NDArray[float64]`	1D array containing simulation time points.	required
`y`	`NDArray[float64]`	1D or 2D array containing the system output trajectory.	required
`x`	`NDArray[float64] \| None`	Optional 2D array containing state trajectories (if simulated from StateSpace).	`None`
`sys`	`Any`	Optional reference to the simulated LTI system.	`None`
`poles`	`NDArray[complex128] \| None`	Optional array of system poles used for stability validation.	`None`

Source code in src/ctrlpy/simulation_results.py

@dataclass
class TimeResponseData:
    """Container for time-domain simulation response data and performance metrics.

    Parameters
    ----------
    t : NDArray[np.float64]
        1D array containing simulation time points.
    y : NDArray[np.float64]
        1D or 2D array containing the system output trajectory.
    x : NDArray[np.float64] | None, optional
        Optional 2D array containing state trajectories (if simulated from StateSpace).
    sys : Any, optional
        Optional reference to the simulated LTI system.
    poles : NDArray[np.complex128] | None, optional
        Optional array of system poles used for stability validation.
    """

    t: NDArray[np.float64]
    y: NDArray[np.float64]
    x: NDArray[np.float64] | None = None
    sys: Any = None
    poles: NDArray[np.complex128] | None = None

    def __post_init__(self) -> None:
        """Validate and cast input arrays to float64."""
        self.t = np.asarray(self.t, dtype=np.float64)
        self.y = np.asarray(self.y, dtype=np.float64)
        if self.x is not None:
            x_arr = np.asarray(self.x, dtype=np.float64)
            if x_arr.ndim == 1:
                x_arr = x_arr.reshape(-1, 1)
            self.x = x_arr

    def __iter__(self) -> Iterator[NDArray[np.float64]]:
        """Enable unpacking of (t, y) or (t, y, x).

        Yields
        ------
        NDArray[np.float64]
            Time array `t`, output array `y`, and optionally state array `x` if present.
        """
        yield self.t
        yield self.y
        if self.x is not None:
            yield self.x

    def _get_channel_1d(self, channel: int = 0) -> NDArray[np.float64]:
        """Extract a 1D output slice for a specified output channel.

        Parameters
        ----------
        channel : int, optional
            Output channel index, default is 0.

        Returns
        -------
        NDArray[np.float64]
            1D array of output values over time.

        Raises
        ------
        IndexError
            If channel index is out of bounds.
        ValueError
            If output array has unsupported dimensionality.
        """
        if self.y.ndim == 1:
            if channel != 0:
                raise IndexError(f"Channel index {channel} is invalid for 1D output.")
            return self.y
        if self.y.ndim == 2:
            if 0 <= channel < self.y.shape[1]:
                return self.y[:, channel]
            raise IndexError(
                f"Channel index {channel} out of range for output shape {self.y.shape}."
            )
        raise ValueError(f"Output array has unsupported dimension {self.y.ndim}.")

    def _check_stability(self) -> None:
        """Check if the system is unstable and raise UnstableSystemError."""
        poles = self.poles
        if poles is None and self.sys is not None:
            poles = self.sys.poles()

        if poles is not None and len(poles) > 0:
            if getattr(self.sys, "is_discrete", False):
                unstable = [p for p in poles if np.abs(p) > 1.0 + 1e-6]
                if len(unstable) > 0:
                    raise UnstableSystemError(
                        f"Cannot calculate transient response metrics for an unstable discrete system. "
                        f"System has pole(s) outside the unit circle (|p| > 1): {unstable}"
                    )
            else:
                unstable = [p for p in poles if np.real(p) > 1e-6]
                if len(unstable) > 0:
                    raise UnstableSystemError(
                        f"Cannot calculate transient response metrics for an unstable system. "
                        f"System has pole(s) in the Right-Half Plane (Re(p) > 0): {unstable}"
                    )

    def steady_state_value(self, channel: int = 0) -> float:
        r"""Compute the final steady-state response value $y_{ss}$.

        .. math::

            y_{ss} = \lim_{t \to \infty} y(t) \approx y(t_{\text{final}})

        Parameters
        ----------
        channel : int, optional
            Output channel index, default is 0.

        Returns
        -------
        float
            The final output value at the end of the simulation horizon.

        Raises
        ------
        UnstableSystemError
            If the system is unstable (has poles with $\mathrm{Re}(p) > 0$).
        """
        self._check_stability()
        y_vec = self._get_channel_1d(channel)
        return float(y_vec[-1])

    def rise_time(
        self,
        low: float = 0.1,
        high: float = 0.9,
        channel: int = 0,
    ) -> float:
        r"""Compute the rise time $t_r$ from low% to high% of the steady-state transition.

        By default, computes the 10% to 90% rise time with linear interpolation
        between discrete time points for high numerical accuracy:

        .. math::

            t_r = t_{90\%} - t_{10\%}

        Parameters
        ----------
        low : float, optional
            Lower threshold fraction of the transition (default is 0.1 for 10%).
        high : float, optional
            Upper threshold fraction of the transition (default is 0.9 for 90%).
        channel : int, optional
            Output channel index, default is 0.

        Returns
        -------
        float
            Rise time in seconds, or NaN if the response does not cross the thresholds.

        Raises
        ------
        UnstableSystemError
            If the system is unstable (has poles with $\mathrm{Re}(p) > 0$).
        """
        self._check_stability()
        y_vec = self._get_channel_1d(channel)
        t_vec = self.t
        y0 = float(y_vec[0])
        yss = float(y_vec[-1])
        dy = yss - y0

        if abs(dy) < 1e-12:
            return 0.0

        y_low = y0 + low * dy
        y_high = y0 + high * dy

        if dy > 0:
            idx_low = np.flatnonzero(y_vec >= y_low)
            idx_high = np.flatnonzero(y_vec >= y_high)
        else:
            idx_low = np.flatnonzero(y_vec <= y_low)
            idx_high = np.flatnonzero(y_vec <= y_high)

        if idx_low.size == 0 or idx_high.size == 0:
            return float("nan")

        i_l = int(idx_low[0])
        if i_l == 0:
            t_l = float(t_vec[0])
        else:
            y_prev, y_curr = float(y_vec[i_l - 1]), float(y_vec[i_l])
            t_prev, t_curr = float(t_vec[i_l - 1]), float(t_vec[i_l])
            t_l = t_prev + (t_curr - t_prev) * (y_low - y_prev) / (y_curr - y_prev)

        i_h = int(idx_high[0])
        if i_h == 0:
            t_h = float(t_vec[0])
        else:
            y_prev, y_curr = float(y_vec[i_h - 1]), float(y_vec[i_h])
            t_prev, t_curr = float(t_vec[i_h - 1]), float(t_vec[i_h])
            t_h = t_prev + (t_curr - t_prev) * (y_high - y_prev) / (y_curr - y_prev)

        return float(t_h - t_l)

    def settling_time(
        self,
        tolerance: float = 0.02,
        channel: int = 0,
    ) -> float:
        r"""Compute the settling time $t_s$ within a specified tolerance band around steady-state.

        Settling time is the time after which the response remains permanently within
        the error band:

        .. math::

            |y(t) - y_{ss}| \le \text{tolerance} \cdot |y_{ss} - y_0|, \quad \forall t \ge t_s

        The boundary crossing is linearly interpolated for sub-sample precision.

        Parameters
        ----------
        tolerance : float, optional
            Settling tolerance fraction, default is 0.02 (2% error band).
        channel : int, optional
            Output channel index, default is 0.

        Returns
        -------
        float
            Settling time in seconds, or NaN if the system has not settled by the
            end of the simulation.

        Raises
        ------
        UnstableSystemError
            If the system is unstable (has poles with $\mathrm{Re}(p) > 0$).
        """
        self._check_stability()
        y_vec = self._get_channel_1d(channel)
        t_vec = self.t
        y0 = float(y_vec[0])
        yss = float(y_vec[-1])
        dy = abs(yss - y0)
        delta = tolerance * dy if dy >= 1e-12 else tolerance * max(abs(yss), 1.0)

        # Check if the signal is actively sweeping or divergent at the end
        k_tail = max(2, int(0.02 * len(t_vec)))
        if np.max(np.abs(y_vec[-k_tail:] - yss)) > delta:
            return float("nan")

        outside = np.flatnonzero(np.abs(y_vec - yss) > delta)
        if outside.size == 0:
            return float(t_vec[0])

        last_idx = int(outside[-1])
        if last_idx == len(t_vec) - 1:
            return float("nan")

        k = last_idx
        y_k = float(y_vec[k])
        y_next = float(y_vec[k + 1])
        t_k = float(t_vec[k])
        t_next = float(t_vec[k + 1])
        target_y = yss + np.sign(y_k - yss) * delta

        if abs(y_next - y_k) > 1e-15:
            return float(t_k + (t_next - t_k) * (target_y - y_k) / (y_next - y_k))
        return float(t_next)

    def overshoot(self, channel: int = 0) -> float:
        r"""Compute the maximum percent overshoot $\%OS$ above steady-state.

        .. math::

            \%OS = \frac{y_{\text{peak}} - y_{ss}}{|y_{ss} - y_0|} \times 100\%

        Parameters
        ----------
        channel : int, optional
            Output channel index, default is 0.

        Returns
        -------
        float
            Percent overshoot (e.g. 16.3 for 16.3% overshoot). Returns 0.0 if no
            overshoot occurs.

        Raises
        ------
        UnstableSystemError
            If the system is unstable (has poles with $\mathrm{Re}(p) > 0$).
        """
        self._check_stability()
        y_vec = self._get_channel_1d(channel)
        y0 = float(y_vec[0])
        yss = float(y_vec[-1])
        dy = yss - y0

        if abs(dy) < 1e-12:
            return 0.0

        if dy > 0:
            y_peak = float(np.max(y_vec))
            if y_peak <= yss:
                return 0.0
            return float(((y_peak - yss) / dy) * 100.0)
        else:
            y_peak = float(np.min(y_vec))
            if y_peak >= yss:
                return 0.0
            return float(((yss - y_peak) / abs(dy)) * 100.0)

    def peak_time(self, channel: int = 0) -> float:
        r"""Compute the time $t_p$ at which the maximum peak response occurs.

        .. math::

            t_p = \arg\max_t |y(t) - y_0|

        Parameters
        ----------
        channel : int, optional
            Output channel index, default is 0.

        Returns
        -------
        float
            Peak time in seconds.

        Raises
        ------
        UnstableSystemError
            If the system is unstable (has poles with $\mathrm{Re}(p) > 0$).
        """
        self._check_stability()
        y_vec = self._get_channel_1d(channel)
        t_vec = self.t
        y0 = float(y_vec[0])
        yss = float(y_vec[-1])
        dy = yss - y0

        if dy >= 0:
            k = int(np.argmax(y_vec))
        else:
            k = int(np.argmin(y_vec))

        return float(t_vec[k])

    def _repr_latex_(self) -> str:
        """Return a LaTeX table of transient response performance metrics for Jupyter."""
        try:
            yss_val = self.steady_state_value()
            yss_str = f"{yss_val:.4f}"
            tr_val = self.rise_time()
            tr_str = f"{tr_val:.4f}\\text{{ s}}" if not np.isnan(tr_val) else r"\text{N/A}"
            ts_val = self.settling_time(tolerance=0.02)
            ts_str = f"{ts_val:.4f}\\text{{ s}}" if not np.isnan(ts_val) else r"\text{N/A}"
            os_val = self.overshoot()
            os_str = f"{os_val:.2f}\\%"
            tp_val = self.peak_time()
            tp_str = f"{tp_val:.4f}\\text{{ s}}"
        except UnstableSystemError:
            return r"$$\textbf{\textcolor{red}{Unstable System: Transient metrics not applicable (divergent response)}}$$"
        except (ValueError, IndexError, RuntimeError, TypeError):
            final_y = self.y[-1] if self.y.ndim == 1 else self.y[-1, 0]
            return (
                rf"$$\text{{TimeResponseData: {len(self.t)} points, final value: {final_y:.4f}}}$$"
            )

        lines = [
            r"$$\begin{array}{|l|c|}",
            r"\hline",
            r"\textbf{Transient Response Metric} & \textbf{Value} \\",
            r"\hline",
            rf"\text{{Steady-State Value }} (y_{{ss}}) & {yss_str} \\",
            rf"\text{{Rise Time }} (t_r, 10\%-90\%) & {tr_str} \\",
            rf"\text{{Settling Time }} (t_s, 2\%) & {ts_str} \\",
            rf"\text{{Percent Overshoot }} (\%OS) & {os_str} \\",
            rf"\text{{Peak Time }} (t_p) & {tp_str} \\",
            r"\hline",
            r"\end{array}$$",
        ]
        return "\n".join(lines)

    def _repr_markdown_(self) -> str:
        """Return a Markdown representation for Jupyter environments."""
        return self._repr_latex_()

__post_init__ ¶

__post_init__()

Validate and cast input arrays to float64.

Source code in src/ctrlpy/simulation_results.py

def __post_init__(self) -> None:
    """Validate and cast input arrays to float64."""
    self.t = np.asarray(self.t, dtype=np.float64)
    self.y = np.asarray(self.y, dtype=np.float64)
    if self.x is not None:
        x_arr = np.asarray(self.x, dtype=np.float64)
        if x_arr.ndim == 1:
            x_arr = x_arr.reshape(-1, 1)
        self.x = x_arr

iter ¶

__iter__()

Enable unpacking of (t, y) or (t, y, x).

Yields:

Type	Description
`NDArray[float64]`	Time array `t`, output array `y`, and optionally state array `x` if present.

Source code in src/ctrlpy/simulation_results.py

def __iter__(self) -> Iterator[NDArray[np.float64]]:
    """Enable unpacking of (t, y) or (t, y, x).

    Yields
    ------
    NDArray[np.float64]
        Time array `t`, output array `y`, and optionally state array `x` if present.
    """
    yield self.t
    yield self.y
    if self.x is not None:
        yield self.x

steady_state_value ¶

steady_state_value(channel=0)

Compute the final steady-state response value \(y_{ss}\).

.. math::

y_{ss} = \lim_{t \to \infty} y(t) \approx y(t_{\text{final}})

Parameters:

Name	Type	Description	Default
`channel`	`int`	Output channel index, default is 0.	`0`

Returns:

Type	Description
`float`	The final output value at the end of the simulation horizon.

Raises:

Type	Description
`UnstableSystemError`	If the system is unstable (has poles with \(\mathrm{Re}(p) > 0\)).

Source code in src/ctrlpy/simulation_results.py

def steady_state_value(self, channel: int = 0) -> float:
    r"""Compute the final steady-state response value $y_{ss}$.

    .. math::

        y_{ss} = \lim_{t \to \infty} y(t) \approx y(t_{\text{final}})

    Parameters
    ----------
    channel : int, optional
        Output channel index, default is 0.

    Returns
    -------
    float
        The final output value at the end of the simulation horizon.

    Raises
    ------
    UnstableSystemError
        If the system is unstable (has poles with $\mathrm{Re}(p) > 0$).
    """
    self._check_stability()
    y_vec = self._get_channel_1d(channel)
    return float(y_vec[-1])

rise_time ¶

rise_time(low=0.1, high=0.9, channel=0)

Compute the rise time \(t_r\) from low% to high% of the steady-state transition.

By default, computes the 10% to 90% rise time with linear interpolation between discrete time points for high numerical accuracy:

.. math::

t_r = t_{90\%} - t_{10\%}

Parameters:

Name	Type	Description	Default
`low`	`float`	Lower threshold fraction of the transition (default is 0.1 for 10%).	`0.1`
`high`	`float`	Upper threshold fraction of the transition (default is 0.9 for 90%).	`0.9`
`channel`	`int`	Output channel index, default is 0.	`0`

Returns:

Type	Description
`float`	Rise time in seconds, or NaN if the response does not cross the thresholds.

Raises:

Type	Description
`UnstableSystemError`	If the system is unstable (has poles with \(\mathrm{Re}(p) > 0\)).

Source code in src/ctrlpy/simulation_results.py

def rise_time(
    self,
    low: float = 0.1,
    high: float = 0.9,
    channel: int = 0,
) -> float:
    r"""Compute the rise time $t_r$ from low% to high% of the steady-state transition.

    By default, computes the 10% to 90% rise time with linear interpolation
    between discrete time points for high numerical accuracy:

    .. math::

        t_r = t_{90\%} - t_{10\%}

    Parameters
    ----------
    low : float, optional
        Lower threshold fraction of the transition (default is 0.1 for 10%).
    high : float, optional
        Upper threshold fraction of the transition (default is 0.9 for 90%).
    channel : int, optional
        Output channel index, default is 0.

    Returns
    -------
    float
        Rise time in seconds, or NaN if the response does not cross the thresholds.

    Raises
    ------
    UnstableSystemError
        If the system is unstable (has poles with $\mathrm{Re}(p) > 0$).
    """
    self._check_stability()
    y_vec = self._get_channel_1d(channel)
    t_vec = self.t
    y0 = float(y_vec[0])
    yss = float(y_vec[-1])
    dy = yss - y0

    if abs(dy) < 1e-12:
        return 0.0

    y_low = y0 + low * dy
    y_high = y0 + high * dy

    if dy > 0:
        idx_low = np.flatnonzero(y_vec >= y_low)
        idx_high = np.flatnonzero(y_vec >= y_high)
    else:
        idx_low = np.flatnonzero(y_vec <= y_low)
        idx_high = np.flatnonzero(y_vec <= y_high)

    if idx_low.size == 0 or idx_high.size == 0:
        return float("nan")

    i_l = int(idx_low[0])
    if i_l == 0:
        t_l = float(t_vec[0])
    else:
        y_prev, y_curr = float(y_vec[i_l - 1]), float(y_vec[i_l])
        t_prev, t_curr = float(t_vec[i_l - 1]), float(t_vec[i_l])
        t_l = t_prev + (t_curr - t_prev) * (y_low - y_prev) / (y_curr - y_prev)

    i_h = int(idx_high[0])
    if i_h == 0:
        t_h = float(t_vec[0])
    else:
        y_prev, y_curr = float(y_vec[i_h - 1]), float(y_vec[i_h])
        t_prev, t_curr = float(t_vec[i_h - 1]), float(t_vec[i_h])
        t_h = t_prev + (t_curr - t_prev) * (y_high - y_prev) / (y_curr - y_prev)

    return float(t_h - t_l)

settling_time ¶

settling_time(tolerance=0.02, channel=0)

Compute the settling time \(t_s\) within a specified tolerance band around steady-state.

Settling time is the time after which the response remains permanently within the error band:

.. math::

|y(t) - y_{ss}| \le \text{tolerance} \cdot |y_{ss} - y_0|, \quad \forall t \ge t_s

The boundary crossing is linearly interpolated for sub-sample precision.

Parameters:

Name	Type	Description	Default
`tolerance`	`float`	Settling tolerance fraction, default is 0.02 (2% error band).	`0.02`
`channel`	`int`	Output channel index, default is 0.	`0`

Returns:

Type	Description
`float`	Settling time in seconds, or NaN if the system has not settled by the end of the simulation.

Raises:

Type	Description
`UnstableSystemError`	If the system is unstable (has poles with \(\mathrm{Re}(p) > 0\)).

Source code in src/ctrlpy/simulation_results.py

def settling_time(
    self,
    tolerance: float = 0.02,
    channel: int = 0,
) -> float:
    r"""Compute the settling time $t_s$ within a specified tolerance band around steady-state.

    Settling time is the time after which the response remains permanently within
    the error band:

    .. math::

        |y(t) - y_{ss}| \le \text{tolerance} \cdot |y_{ss} - y_0|, \quad \forall t \ge t_s

    The boundary crossing is linearly interpolated for sub-sample precision.

    Parameters
    ----------
    tolerance : float, optional
        Settling tolerance fraction, default is 0.02 (2% error band).
    channel : int, optional
        Output channel index, default is 0.

    Returns
    -------
    float
        Settling time in seconds, or NaN if the system has not settled by the
        end of the simulation.

    Raises
    ------
    UnstableSystemError
        If the system is unstable (has poles with $\mathrm{Re}(p) > 0$).
    """
    self._check_stability()
    y_vec = self._get_channel_1d(channel)
    t_vec = self.t
    y0 = float(y_vec[0])
    yss = float(y_vec[-1])
    dy = abs(yss - y0)
    delta = tolerance * dy if dy >= 1e-12 else tolerance * max(abs(yss), 1.0)

    # Check if the signal is actively sweeping or divergent at the end
    k_tail = max(2, int(0.02 * len(t_vec)))
    if np.max(np.abs(y_vec[-k_tail:] - yss)) > delta:
        return float("nan")

    outside = np.flatnonzero(np.abs(y_vec - yss) > delta)
    if outside.size == 0:
        return float(t_vec[0])

    last_idx = int(outside[-1])
    if last_idx == len(t_vec) - 1:
        return float("nan")

    k = last_idx
    y_k = float(y_vec[k])
    y_next = float(y_vec[k + 1])
    t_k = float(t_vec[k])
    t_next = float(t_vec[k + 1])
    target_y = yss + np.sign(y_k - yss) * delta

    if abs(y_next - y_k) > 1e-15:
        return float(t_k + (t_next - t_k) * (target_y - y_k) / (y_next - y_k))
    return float(t_next)

overshoot ¶

overshoot(channel=0)

Compute the maximum percent overshoot \(\%OS\) above steady-state.

.. math::

\%OS = \frac{y_{\text{peak}} - y_{ss}}{|y_{ss} - y_0|} \times 100\%

Parameters:

Name	Type	Description	Default
`channel`	`int`	Output channel index, default is 0.	`0`

Returns:

Type	Description
`float`	Percent overshoot (e.g. 16.3 for 16.3% overshoot). Returns 0.0 if no overshoot occurs.

Raises:

Type	Description
`UnstableSystemError`	If the system is unstable (has poles with \(\mathrm{Re}(p) > 0\)).

Source code in src/ctrlpy/simulation_results.py

def overshoot(self, channel: int = 0) -> float:
    r"""Compute the maximum percent overshoot $\%OS$ above steady-state.

    .. math::

        \%OS = \frac{y_{\text{peak}} - y_{ss}}{|y_{ss} - y_0|} \times 100\%

    Parameters
    ----------
    channel : int, optional
        Output channel index, default is 0.

    Returns
    -------
    float
        Percent overshoot (e.g. 16.3 for 16.3% overshoot). Returns 0.0 if no
        overshoot occurs.

    Raises
    ------
    UnstableSystemError
        If the system is unstable (has poles with $\mathrm{Re}(p) > 0$).
    """
    self._check_stability()
    y_vec = self._get_channel_1d(channel)
    y0 = float(y_vec[0])
    yss = float(y_vec[-1])
    dy = yss - y0

    if abs(dy) < 1e-12:
        return 0.0

    if dy > 0:
        y_peak = float(np.max(y_vec))
        if y_peak <= yss:
            return 0.0
        return float(((y_peak - yss) / dy) * 100.0)
    else:
        y_peak = float(np.min(y_vec))
        if y_peak >= yss:
            return 0.0
        return float(((yss - y_peak) / abs(dy)) * 100.0)

peak_time ¶

peak_time(channel=0)

Compute the time \(t_p\) at which the maximum peak response occurs.

.. math::

t_p = \arg\max_t |y(t) - y_0|

Parameters:

Name	Type	Description	Default
`channel`	`int`	Output channel index, default is 0.	`0`

Returns:

Type	Description
`float`	Peak time in seconds.

Raises:

Type	Description
`UnstableSystemError`	If the system is unstable (has poles with \(\mathrm{Re}(p) > 0\)).

Source code in src/ctrlpy/simulation_results.py

def peak_time(self, channel: int = 0) -> float:
    r"""Compute the time $t_p$ at which the maximum peak response occurs.

    .. math::

        t_p = \arg\max_t |y(t) - y_0|

    Parameters
    ----------
    channel : int, optional
        Output channel index, default is 0.

    Returns
    -------
    float
        Peak time in seconds.

    Raises
    ------
    UnstableSystemError
        If the system is unstable (has poles with $\mathrm{Re}(p) > 0$).
    """
    self._check_stability()
    y_vec = self._get_channel_1d(channel)
    t_vec = self.t
    y0 = float(y_vec[0])
    yss = float(y_vec[-1])
    dy = yss - y0

    if dy >= 0:
        k = int(np.argmax(y_vec))
    else:
        k = int(np.argmin(y_vec))

    return float(t_vec[k])

Time-Domain Simulation & Metrics¶

time_domain ¶

step_response ¶

impulse_response ¶

forced_response ¶

simulation_results ¶

TimeResponseData dataclass ¶

__post_init__ ¶

__iter__ ¶

steady_state_value ¶

rise_time ¶

settling_time ¶

overshoot ¶

peak_time ¶

TimeResponseData `dataclass` ¶

iter ¶