In control engineering, we often summarize a time response by a small set of scalars that capture “speed” and “quality.” For a step input, the most common performance metrics are:

• Rise time \( t_r \): how quickly the response moves from near its initial value to near its steady value.
• Settling time \( t_s \): how long it takes to enter (and remain in) a specified error band around steady state.
• Overshoot \( M_p \): how much the transient exceeds the steady value.
• Steady-state error \( e_{ss} \): the asymptotic difference between input and output (tracking error) or between desired and achieved final value.

flowchart LR
  A["Input: step of amplitude A"] --> B["Output y(t)"]
  B --> C["y_inf = final value"]
  B --> D["t_r: between 10% and 90% of y_inf (common)"]
  B --> E["t_p: time of first peak"]
  B --> F["M_p: (y_peak - y_inf)/y_inf"]
  B --> G["t_s: first time response stays \nwithin +/-2% (or 5%) band"]
        

Throughout this lesson we focus on step responses, because they are a canonical excitation and provide a consistent basis for comparing systems. The definitions and derivations extend naturally to other inputs, but step-based metrics are the standard foundation.

Let \( r(t) \) be a step input of amplitude \( A \), i.e. \( r(t)=A\,u(t) \), where \( u(t) \) is the Heaviside step. Let \( y(t) \) be the output. Define the (tracking) error: \( e(t)=r(t)-y(t) \).

A key quantity is the steady (final) value \( y_\infty \), defined as the limit (when it exists):

\[ y_\infty \triangleq \lim_{t \to \infty} y(t). \]

If the system is LTI with transfer function \( G(s) \) from \( r \) to \( y \), then for a step input \( R(s)=A/s \), we have \( Y(s)=G(s)R(s) \). Under standard conditions for the final value theorem (namely that all poles of \( sY(s) \) satisfy \( \Re(p) < 0 \)),

\[ y_\infty = \lim_{s \to 0} sY(s) = \lim_{s \to 0} s\,G(s)\frac{A}{s} = A \lim_{s \to 0} G(s) = A\,G(0). \]

Therefore, for a stable system with finite DC gain \( G(0) \), the steady output to a step is determined purely by the DC gain.

The steady-state tracking error (for the same input-output channel) is:

\[ e_{ss} \triangleq \lim_{t \to \infty} e(t) = \lim_{t \to \infty} \big(r(t) - y(t)\big) = A - y_\infty. \]

In many control configurations, the step is applied as a reference through a feedback loop. For unity feedback with forward-path \( G(s) \), the error transfer is \( E(s)=R(s)/(1+G(s)) \), yielding:

\[ e_{ss} = \lim_{s \to 0} sE(s) = \lim_{s \to 0} s \frac{A/s}{1+G(s)} = \frac{A}{1+ \lim_{s \to 0} G(s)}. \]

This motivates the DC loop-gain constant \( K_p \triangleq \lim_{s \to 0}G(s) \) (when finite), giving the compact expression \( e_{ss} = A/(1+K_p) \).

Assume a step response with finite steady value \( y_\infty \) and initial value \( y(0^+) \). Define the normalized response: \( \tilde{y}(t) = \frac{y(t)-y(0^+)}{y_\infty - y(0^+)} \) when \( y_\infty \neq y(0^+) \).

Rise time. A common convention is “10% to 90% rise time”: define \( t_{10} \) as the first time \( \tilde{y}(t) \) reaches \(0.1\) and \( t_{90} \) as the first time it reaches \(0.9\). Then:

\[ t_r \triangleq t_{90} - t_{10}, \quad t_{10} \triangleq \inf\{t \ge 0 : \tilde{y}(t) \ge 0.1\},\quad t_{90} \triangleq \inf\{t \ge 0 : \tilde{y}(t) \ge 0.9\}. \]

Overshoot. Let \( y_{\max} \triangleq \sup_{t \ge 0} y(t) \). The peak overshoot ratio is:

\[ M_p \triangleq \frac{y_{\max}-y_\infty}{|y_\infty - y(0^+)|}, \quad \text{and percent overshoot } \%OS = 100\,M_p. \]

Settling time. Fix a tolerance band \( \delta \in (0,1) \) (e.g., \(0.02\) for 2%). Define the settling time as:

\[ t_s(\delta) \triangleq \inf\left\{T \ge 0 : \forall t \ge T,\; \left|\frac{y(t)-y_\infty}{y_\infty - y(0^+)}\right| \le \delta \right\}. \]

These definitions are robust and apply to any response trajectory. Closed-form formulas emerge for standard first- and second-order models, which we derive next.

Consider the standard first-order transfer function (introduced in Lesson 1):

\[ G(s) = \frac{K}{\tau s + 1}, \quad K \in \mathbb{R},\; \tau > 0. \]

For a step input of amplitude \(A\), the time response (for zero initial conditions) is:

\[ y(t) = KA\left(1 - e^{-t/\tau}\right), \quad t \ge 0, \quad y_\infty = KA. \]

10%–90% rise time. Solve \( y(t_{p}) = p\,y_\infty \) for \(p \in (0,1)\):

\[ KA\left(1 - e^{-t_p/\tau}\right)=p\,KA \;\Rightarrow\; e^{-t_p/\tau}=1-p \;\Rightarrow\; t_p = -\tau \ln(1-p). \]

Therefore, \( t_{10} = -\tau\ln(0.9) \) and \( t_{90} = -\tau\ln(0.1) \), giving:

\[ t_r = t_{90}-t_{10} = -\tau\ln(0.1) + \tau\ln(0.9) = \tau \ln\!\left(\frac{0.9}{0.1}\right) = \tau \ln(9). \]

Settling time (2% band). Require \( |y(t)-y_\infty| \le 0.02|y_\infty|\):

\[ |KAe^{-t/\tau}| \le 0.02|KA| \;\Rightarrow\; e^{-t/\tau} \le 0.02 \;\Rightarrow\; t \ge -\tau \ln(0.02). \]

Hence \( t_s(0.02) = -\tau\ln(0.02) \approx 3.912\,\tau \).

Overshoot. Since \( y(t) \) is monotone increasing for \(K A > 0\), we have \( y_{\max}=y_\infty \) and \( M_p=0 \).

Steady-state error. For this input-output mapping, \( e_{ss}=A-y_\infty=A(1-K) \). In unity feedback, the earlier result gives \( e_{ss} = A/(1+K_p) \), with \( K_p = \lim_{s \to 0} \frac{K}{\tau s+1} = K \).

Consider the standard second-order closed-loop form (Lesson 2):

\[ T(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}, \quad \omega_n > 0,\; 0 < \zeta < 1. \]

For a unit step input, the response is \( y(t)=\mathcal{L}^{-1}\{T(s)/s\} \). A standard derivation proceeds via the equivalent ODE:

\[ \ddot{y}(t) + 2\zeta\omega_n \dot{y}(t) + \omega_n^2 y(t) = \omega_n^2 u(t),\quad y(0)=0,\;\dot{y}(0)=0. \]

Step-response formula (sketch of proof). The particular solution for a constant RHS is \( y_p(t)=1 \). The homogeneous solution solves \( r^2 + 2\zeta\omega_n r + \omega_n^2 = 0 \), with complex roots:

\[ r_{1,2} = -\zeta\omega_n \pm j\omega_d,\quad \omega_d \triangleq \omega_n\sqrt{1-\zeta^2}. \]

Hence \( y_h(t)=e^{-\zeta\omega_n t}\big(C_1\cos(\omega_d t)+C_2\sin(\omega_d t)\big) \). Applying \( y(0)=0 \) and \( \dot{y}(0)=0 \) yields constants that can be expressed compactly using the angle \( \phi \) defined by:

\[ \phi \triangleq \arctan\!\left(\frac{\sqrt{1-\zeta^2}}{\zeta}\right), \quad \text{so that}\quad \sin(\phi)=\sqrt{1-\zeta^2},\;\cos(\phi)=\zeta. \]

One obtains the standard closed form:

\[ y(t) = 1 - \frac{1}{\sqrt{1-\zeta^2}} e^{-\zeta\omega_n t}\, \sin(\omega_d t + \phi),\quad t \ge 0. \]

Peak time \(t_p\) and overshoot \(M_p\). Differentiate:

\[ \dot{y}(t)=\frac{1}{\sqrt{1-\zeta^2}}e^{-\zeta\omega_n t} \left(\zeta\omega_n\sin(\omega_d t+\phi) - \omega_d\cos(\omega_d t+\phi)\right). \]

A (local) extremum occurs when the bracket equals zero: \( \zeta\omega_n\sin(\theta) = \omega_d\cos(\theta) \) with \( \theta=\omega_d t + \phi \), i.e. \( \tan(\theta) = \omega_d/(\zeta\omega_n)=\sqrt{1-\zeta^2}/\zeta \). Using the definition of \( \phi \), this implies \( \theta = \pi \) at the first positive peak, hence:

\[ \omega_d t_p + \phi = \pi \;\Rightarrow\; t_p = \frac{\pi - \phi}{\omega_d}. \]

Many texts present the commonly used approximation \( t_p \approx \pi/\omega_d \) because \( \phi \) is often moderate; however, the exact first peak time for the first overshoot relative to unity final value is:

\[ t_p = \frac{\pi}{\omega_d}. \]

(This follows equivalently by analyzing the first maximum of \( \sin(\omega_d t + \phi) \) relative to the decaying envelope; the classical result is \( t_p=\pi/\omega_d \).)

Evaluate \( y(t) \) at \( t_p=\pi/\omega_d \). Since \( \sin(\omega_d t_p + \phi)=\sin(\pi+\phi)=-\sin(\phi) \), we get:

\[ y(t_p) = 1 - \frac{1}{\sqrt{1-\zeta^2}} e^{-\zeta\omega_n (\pi/\omega_d)}\,(-\sin\phi) = 1 + e^{-\frac{\zeta\pi}{\sqrt{1-\zeta^2}}}. \]

Therefore, the overshoot ratio is:

\[ M_p = y(t_p)-1 = \exp\!\left(-\frac{\zeta\pi}{\sqrt{1-\zeta^2}}\right), \quad \%OS = 100\,\exp\!\left(-\frac{\zeta\pi}{\sqrt{1-\zeta^2}}\right). \]

Settling time \(t_s(\delta)\) via an envelope bound. From the response expression:

\[ |y(t)-1| = \left|\frac{1}{\sqrt{1-\zeta^2}}e^{-\zeta\omega_n t} \sin(\omega_d t + \phi)\right| \le \frac{1}{\sqrt{1-\zeta^2}}e^{-\zeta\omega_n t}. \]

To guarantee \( |y(t)-1| \le \delta \) for all later times, it is sufficient to enforce:

\[ \frac{1}{\sqrt{1-\zeta^2}}e^{-\zeta\omega_n t} \le \delta \;\Rightarrow\; t \ge \frac{1}{\zeta\omega_n}\ln\!\left(\frac{1}{\delta\sqrt{1-\zeta^2}}\right). \]

Thus an analytically justified conservative settling-time bound is:

\[ t_s(\delta) \le \frac{1}{\zeta\omega_n}\ln\!\left(\frac{1}{\delta\sqrt{1-\zeta^2}}\right). \]

For the common 2% band and moderate damping, the widely used design approximation is:

\[ t_s(0.02) \approx \frac{4}{\zeta\omega_n}. \]

Rise time. For underdamped responses, “0%–100% rise time” is often defined as the first time the response reaches the final value, i.e. the smallest \( t \) with \( y(t)=1 \). From the response formula:

\[ y(t)=1 \;\Rightarrow\; \sin(\omega_d t + \phi)=0 \;\Rightarrow\; \omega_d t + \phi = \pi \;\Rightarrow\; t_r^{(0\text{–}100)}=\frac{\pi-\phi}{\omega_d}. \]

The 10%–90% rise time generally requires solving \( y(t)=0.1 \) and \( y(t)=0.9 \), for which numerical root-finding is typically used.

From Lesson 4, many stable higher-order LTI systems can be approximated (for transient behavior) by their dominant poles. A generic stable step response can often be written (conceptually) as:

\[ y(t) = y_\infty + \sum_{k=1}^{n} c_k e^{p_k t}, \quad \Re(p_k) < 0. \]

Define the slowest decay rate \( \sigma \triangleq \min_k (-\Re(p_k)) \), so that \( \sigma > 0 \). Then \( \Re(p_k) \le -\sigma \) for all \(k\). Using the triangle inequality:

\[ |y(t)-y_\infty| = \left|\sum_{k=1}^n c_k e^{p_k t}\right| \le \sum_{k=1}^n |c_k|\,|e^{p_k t}| = \sum_{k=1}^n |c_k|\,e^{\Re(p_k)t} \le \left(\sum_{k=1}^n |c_k|\right)e^{-\sigma t}. \]

Let \( C \triangleq \sum_{k=1}^n |c_k| \). A sufficient condition to be within a relative band \( |y(t)-y_\infty| \le \delta |y_\infty| \) is:

\[ C e^{-\sigma t} \le \delta |y_\infty| \;\Rightarrow\; t \ge \frac{1}{\sigma}\ln\!\left(\frac{C}{\delta |y_\infty|}\right). \]

This inequality makes the dominant-pole intuition precise: the smallest magnitude real part (closest pole to the imaginary axis) governs the asymptotic transient decay and therefore dominates settling time, up to a logarithmic constant.

In practice, when a dominant complex conjugate pair exists, engineers often map the pair to an equivalent second-order model and use the closed-form approximations for \(t_s\) and \(M_p\).

In simulation and experiments, \( y(t) \) is typically available as discrete samples \( (t_i, y_i) \). A robust approach estimates \( y_\infty \) from the tail of the response and then computes thresholds and extrema.

flowchart TD
  S["Obtain sampled response: (t[i], y[i])"] --> F["Estimate y_inf from last N samples"]
  F --> N1["Normalize: y_norm = (y - y0)/(y_inf - y0)"]
  N1 --> R1["Find first crossing of 0.1 \nand 0.9 -> t10, t90"]
  R1 --> TR["Rise time: t_r = t90 - t10"]
  N1 --> OS["Overshoot: y_peak = max(y); \nM_p = (y_peak - y_inf)/abs(y_inf - y0)"]
  N1 --> TS["Settling: find smallest T such that \n|y(t)-y_inf| <= delta*|y_inf-y0| \nfor all t>=T"]
  F --> ESS["Steady-state error: e_ss = A - y_inf \n(or reference minus final output)"]
        

Two practical details matter:

• Interpolation: threshold crossings are between samples; linear interpolation reduces discretization bias.
• Settling logic: scan backward from the end to find the last time the response exits the tolerance band; the settling time is the time immediately after this last exit.

We use two approaches: (i) library-based simulation via python-control and (ii) metric computation from scratch.

import numpy as np

def step_metrics(t, y, A=1.0, delta=0.02, tail_frac=0.1):
    """
    Compute rise time (10-90), settling time (delta band), overshoot, steady-state error.
    Assumes t is increasing and y is step response to a step of amplitude A.
    """
    t = np.asarray(t).astype(float)
    y = np.asarray(y).astype(float)
    if t.ndim != 1 or y.ndim != 1 or t.size != y.size:
        raise ValueError("t and y must be 1D arrays of the same length")

    # Estimate y_inf from the last tail_frac portion
    n = t.size
    tail_n = max(5, int(np.ceil(tail_frac * n)))
    y_inf = float(np.mean(y[-tail_n:]))

    y0 = float(y[0])
    denom = (y_inf - y0)
    if abs(denom) < 1e-12:
        raise ValueError("y_inf is too close to y0; cannot normalize")

    y_norm = (y - y0) / denom

    def first_cross(level):
        idx = np.where(y_norm >= level)[0]
        if idx.size == 0:
            return np.nan
        k = int(idx[0])
        if k == 0:
            return float(t[0])
        # linear interpolation between k-1 and k
        t1, t2 = t[k-1], t[k]
        y1, y2 = y_norm[k-1], y_norm[k]
        if abs(y2 - y1) < 1e-15:
            return float(t2)
        return float(t1 + (level - y1) * (t2 - t1) / (y2 - y1))

    t10 = first_cross(0.1)
    t90 = first_cross(0.9)
    tr = t90 - t10 if np.isfinite(t10) and np.isfinite(t90) else np.nan

    # Overshoot
    y_peak = float(np.max(y))
    Mp = (y_peak - y_inf) / abs(denom)

    # Settling time: find last index outside band, then settling is next time sample
    band = delta * abs(denom)
    outside = np.where(np.abs(y - y_inf) > band)[0]
    if outside.size == 0:
        ts = float(t[0])
    else:
        last_out = int(outside[-1])
        ts = float(t[min(last_out + 1, n - 1)])

    ess = float(A - y_inf)

    return {
        "y_inf": y_inf,
        "t10": t10, "t90": t90, "t_r_10_90": tr,
        "t_s": ts, "M_p": Mp,
        "e_ss": ess,
        "y_peak": y_peak
    }

# --- Example: second-order standard form using python-control (if installed) ---
# pip install control
try:
    import control as ct

    zeta = 0.3
    wn = 10.0
    num = [wn**2]
    den = [1.0, 2*zeta*wn, wn**2]
    sys = ct.TransferFunction(num, den)

    t = np.linspace(0, 3, 3000)
    t, y = ct.step_response(sys, T=t)

    m = step_metrics(t, y, A=1.0, delta=0.02)
    for k, v in m.items():
        print(f"{k}: {v}")

except Exception as ex:
    print("python-control not available or failed to run:", ex)
      

Notes: For underdamped systems, rise time depends strongly on the chosen convention (0–100 vs. 10–90) and is therefore best computed from data. Settling time depends on the tolerance band \( \delta \).

This C++ example generates the underdamped second-order step response analytically and computes the same metrics. Compile with a modern C++ compiler (e.g., g++ -O2 main.cpp -o metrics).

#include <iostream>
#include <vector>
#include <cmath>
#include <limits>

struct Metrics {
    double y_inf;
    double t10, t90, tr_10_90;
    double ts;
    double Mp;
    double ess;
    double y_peak;
};

static double first_cross(const std::vector<double>& t,
                          const std::vector<double>& y_norm,
                          double level) {
    for (size_t k = 0; k < y_norm.size(); ++k) {
        if (y_norm[k] >= level) {
            if (k == 0) return t[0];
            double t1 = t[k-1], t2 = t[k];
            double y1 = y_norm[k-1], y2 = y_norm[k];
            double denom = (y2 - y1);
            if (std::abs(denom) < 1e-15) return t2;
            return t1 + (level - y1) * (t2 - t1) / denom;
        }
    }
    return std::numeric_limits<double>::quiet_NaN();
}

static Metrics step_metrics(const std::vector<double>& t,
                            const std::vector<double>& y,
                            double A = 1.0,
                            double delta = 0.02,
                            double tail_frac = 0.1) {
    if (t.size() != y.size() || t.empty()) {
        throw std::runtime_error("t and y must have same nonzero length");
    }

    size_t n = t.size();
    size_t tail_n = std::max<size_t>(5, static_cast<size_t>(std::ceil(tail_frac * n)));

    double y_inf = 0.0;
    for (size_t i = n - tail_n; i < n; ++i) y_inf += y[i];
    y_inf /= static_cast<double>(tail_n);

    double y0 = y[0];
    double denom = (y_inf - y0);
    if (std::abs(denom) < 1e-12) {
        throw std::runtime_error("y_inf too close to y0; cannot normalize");
    }

    std::vector<double> y_norm(n);
    for (size_t i = 0; i < n; ++i) y_norm[i] = (y[i] - y0) / denom;

    double t10 = first_cross(t, y_norm, 0.1);
    double t90 = first_cross(t, y_norm, 0.9);
    double tr = (std::isfinite(t10) && std::isfinite(t90)) ? (t90 - t10)
                                                           : std::numeric_limits<double>::quiet_NaN();

    double y_peak = y[0];
    for (double yi : y) y_peak = std::max(y_peak, yi);

    double Mp = (y_peak - y_inf) / std::abs(denom);

    double band = delta * std::abs(denom);
    long long last_out = -1;
    for (size_t i = 0; i < n; ++i) {
        if (std::abs(y[i] - y_inf) > band) last_out = static_cast<long long>(i);
    }
    double ts = (last_out < 0) ? t[0] : t[std::min<size_t>(static_cast<size_t>(last_out + 1), n - 1)];

    double ess = A - y_inf;

    Metrics m;
    m.y_inf = y_inf;
    m.t10 = t10; m.t90 = t90; m.tr_10_90 = tr;
    m.ts = ts; m.Mp = Mp; m.ess = ess; m.y_peak = y_peak;
    return m;
}

int main() {
    // Underdamped second-order: T(s)=wn^2/(s^2 + 2*zeta*wn*s + wn^2)
    double zeta = 0.3;
    double wn = 10.0;
    double wd = wn * std::sqrt(1.0 - zeta*zeta);
    double phi = std::atan2(std::sqrt(1.0 - zeta*zeta), zeta);

    // y(t) = 1 - (1/sqrt(1-zeta^2)) * exp(-zeta*wn*t) * sin(wd*t + phi)
    auto y_of_t = [&](double t) {
        double a = 1.0 / std::sqrt(1.0 - zeta*zeta);
        return 1.0 - a * std::exp(-zeta*wn*t) * std::sin(wd*t + phi);
    };

    // Sample response
    double T_end = 3.0;
    int N = 3000;
    std::vector<double> t(N), y(N);
    for (int i = 0; i < N; ++i) {
        t[i] = (T_end * i) / (N - 1);
        y[i] = y_of_t(t[i]);
    }

    Metrics m = step_metrics(t, y, 1.0, 0.02);

    std::cout << "y_inf: " << m.y_inf << "\n";
    std::cout << "t10: " << m.t10 << ", t90: " << m.t90 << ", tr(10-90): " << m.tr_10_90 << "\n";
    std::cout << "ts(2%): " << m.ts << "\n";
    std::cout << "Mp: " << m.Mp << " (percent OS: " << 100.0*m.Mp << "%)\n";
    std::cout << "e_ss: " << m.ess << "\n";
    std::cout << "y_peak: " << m.y_peak << "\n";

    return 0;
}
      

This Java example mirrors the C++ approach: generate an analytic second-order step response, then compute rise time, overshoot, settling time, and steady-state error. (In larger projects, numerical ODE solvers such as Apache Commons Math can replace analytic expressions.)

import java.util.Arrays;

public class StepMetrics {

    static class Metrics {
        double yInf;
        double t10, t90, tr;
        double ts;
        double Mp;
        double ess;
        double yPeak;
    }

    static double firstCross(double[] t, double[] yNorm, double level) {
        for (int k = 0; k < yNorm.length; k++) {
            if (yNorm[k] >= level) {
                if (k == 0) return t[0];
                double t1 = t[k - 1], t2 = t[k];
                double y1 = yNorm[k - 1], y2 = yNorm[k];
                double denom = (y2 - y1);
                if (Math.abs(denom) < 1e-15) return t2;
                return t1 + (level - y1) * (t2 - t1) / denom;
            }
        }
        return Double.NaN;
    }

    static Metrics stepMetrics(double[] t, double[] y, double A, double delta, double tailFrac) {
        if (t.length != y.length || t.length == 0) throw new IllegalArgumentException("Bad input arrays");

        int n = t.length;
        int tailN = Math.max(5, (int) Math.ceil(tailFrac * n));

        double yInf = 0.0;
        for (int i = n - tailN; i < n; i++) yInf += y[i];
        yInf /= tailN;

        double y0 = y[0];
        double denom = (yInf - y0);
        if (Math.abs(denom) < 1e-12) throw new IllegalArgumentException("yInf too close to y0");

        double[] yNorm = new double[n];
        for (int i = 0; i < n; i++) yNorm[i] = (y[i] - y0) / denom;

        double t10 = firstCross(t, yNorm, 0.1);
        double t90 = firstCross(t, yNorm, 0.9);
        double tr = (Double.isFinite(t10) && Double.isFinite(t90)) ? (t90 - t10) : Double.NaN;

        double yPeak = y[0];
        for (double yi : y) yPeak = Math.max(yPeak, yi);

        double Mp = (yPeak - yInf) / Math.abs(denom);

        double band = delta * Math.abs(denom);
        int lastOut = -1;
        for (int i = 0; i < n; i++) {
            if (Math.abs(y[i] - yInf) > band) lastOut = i;
        }
        double ts = (lastOut < 0) ? t[0] : t[Math.min(lastOut + 1, n - 1)];

        double ess = A - yInf;

        Metrics m = new Metrics();
        m.yInf = yInf; m.t10 = t10; m.t90 = t90; m.tr = tr;
        m.ts = ts; m.Mp = Mp; m.ess = ess; m.yPeak = yPeak;
        return m;
    }

    public static void main(String[] args) {
        double zeta = 0.3;
        double wn = 10.0;
        double wd = wn * Math.sqrt(1.0 - zeta * zeta);
        double phi = Math.atan2(Math.sqrt(1.0 - zeta * zeta), zeta);

        // y(t) = 1 - (1/sqrt(1-zeta^2)) exp(-zeta wn t) sin(wd t + phi)
        java.util.function.DoubleUnaryOperator yOfT = (double tt) -> {
            double a = 1.0 / Math.sqrt(1.0 - zeta * zeta);
            return 1.0 - a * Math.exp(-zeta * wn * tt) * Math.sin(wd * tt + phi);
        };

        int N = 3000;
        double TEnd = 3.0;
        double[] t = new double[N];
        double[] y = new double[N];
        for (int i = 0; i < N; i++) {
            t[i] = (TEnd * i) / (N - 1);
            y[i] = yOfT.applyAsDouble(t[i]);
        }

        Metrics m = stepMetrics(t, y, 1.0, 0.02, 0.1);

        System.out.println("y_inf: " + m.yInf);
        System.out.println("t10: " + m.t10 + ", t90: " + m.t90 + ", tr(10-90): " + m.tr);
        System.out.println("ts(2%): " + m.ts);
        System.out.println("Mp: " + m.Mp + " (percent OS: " + (100.0 * m.Mp) + "%)");
        System.out.println("e_ss: " + m.ess);
        System.out.println("y_peak: " + m.yPeak);
    }
}
      

MATLAB’s Control System Toolbox provides step and stepinfo, which compute rise time, settling time, overshoot, and more using standardized conventions (configurable via name-value parameters).

% Second-order standard form
zeta = 0.3;
wn = 10;
T = tf([wn^2], [1 2*zeta*wn wn^2]);

% Step response
t = linspace(0, 3, 3000);
y = step(T, t);

% Library metrics (default definitions use 10-90 rise time and 2% settling band)
info = stepinfo(y, t, 1);   % third argument is step amplitude (1)
disp(info);

% Manual steady-state error for a unity step at the output of T(s)
y_inf = mean(y(end-200:end));
e_ss  = 1 - y_inf;
fprintf('Estimated y_inf = %.6f, e_ss = %.6f\n', y_inf, e_ss);

% Unity feedback steady-state error (reference tracking) for plant G(s)
% If using plant G(s) and unity feedback: e_ss = 1/(1+dcgain(G))
G = tf(5, [1 1]);  % Example plant
Kp = dcgain(G);
e_ss_unity_feedback = 1/(1 + Kp);
fprintf('Unity-feedback step e_ss (for G): %.6f\n', e_ss_unity_feedback);
      

Simulink workflow (recommended for later labs):

1. Create a model with blocks: Step → Transfer Fcn → To Workspace (or logging).
2. Set Transfer Fcn numerator to [wn^2] and denominator to [1 2*zeta*wn wn^2].
3. Run the simulation and retrieve t, y from the workspace.
4. Call stepinfo(y,t,1) and/or implement a custom metric function (as in the Python section) to match your course conventions.

Mathematica can generate step responses symbolically/numerically and extract metrics with root finding and maximization.

(* Second-order standard form *)
zeta = 0.3;
wn = 10;
sys = TransferFunctionModel[wn^2/(s^2 + 2 zeta wn s + wn^2), s];

(* Step response *)
tmax = 3;
resp = OutputResponse[sys, UnitStep[t], {t, 0, tmax}];

(* Sample for numerical metric extraction *)
tsamp = Subdivide[0, tmax, 3000];
ysamp = resp /@ tsamp;

y0 = ysamp[[1]];
yInf = Mean[Take[ysamp, -200]];
den = yInf - y0;
yNorm = (ysamp - y0)/den;

firstCross[level_] := Module[{k},
  k = FirstPosition[yNorm, _?(# >= level &), Missing["NotFound"]];
  If[k === Missing["NotFound"], Indeterminate,
    If[k[[1]] == 1, tsamp[[1]],
      (* linear interpolation *)
      With[{i = k[[1]]},
        t1 = tsamp[[i - 1]]; t2 = tsamp[[i]];
        y1 = yNorm[[i - 1]]; y2 = yNorm[[i]];
        t1 + (level - y1) (t2 - t1)/(y2 - y1)
      ]
    ]
  ]
];

t10 = firstCross[0.1];
t90 = firstCross[0.9];
tr = t90 - t10;

yPeak = Max[ysamp];
Mp = (yPeak - yInf)/Abs[den];

delta = 0.02;
band = delta Abs[den];
lastOut = Last@Flatten@Position[Abs[ysamp - yInf], _?(# > band &), {1}, Heads -> False];
ts = If[MissingQ[lastOut], 0, tsamp[[Min[lastOut + 1, Length[tsamp]]]]];

ess = 1 - yInf;

{yInf, t10, t90, tr, ts, Mp, 100 Mp, ess}
      

Problem 1 (First-order rise and settling time): A first-order system has \( G(s)=\frac{1}{\tau s + 1} \) with \( \tau = 0.5 \) seconds. For a unit step input, compute \( t_r \) (10%–90%) and \( t_s \) (2% band).

Solution:

\[ t_r = \tau \ln(9) = 0.5\ln(9) \approx 1.0986\;\text{s}. \]

\[ t_s(0.02) = -\tau\ln(0.02) = -0.5\ln(0.02) \approx 1.956\;\text{s}. \]

Problem 2 (Overshoot and peak time from \( \zeta, \omega_n \)): Consider \( T(s)=\frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2} \) with \( \zeta=0.4 \) and \( \omega_n=8 \) rad/s. Compute \( \%OS \) and \( t_p \).

Solution: First compute

\[ \omega_d = \omega_n\sqrt{1-\zeta^2} = 8\sqrt{1-0.16} = 8\sqrt{0.84} \approx 7.330. \]

\[ t_p = \frac{\pi}{\omega_d} \approx \frac{\pi}{7.330} \approx 0.429\;\text{s}. \]

\[ \%OS = 100\exp\!\left(-\frac{\zeta\pi}{\sqrt{1-\zeta^2}}\right) = 100\exp\!\left(-\frac{0.4\pi}{\sqrt{0.84}}\right) \approx 100 e^{-1.371} \approx 25.4\%. \]

Problem 3 (Derive the overshoot formula): Starting from \( y(t) = 1 - \frac{1}{\sqrt{1-\zeta^2}} e^{-\zeta\omega_n t}\sin(\omega_d t + \phi) \), show that the overshoot ratio is \( M_p = \exp\!\left(-\frac{\zeta\pi}{\sqrt{1-\zeta^2}}\right) \).

Solution:

At the first peak, the standard result is \( t_p = \pi/\omega_d \). Substitute into the response:

\[ y(t_p) = 1 - \frac{1}{\sqrt{1-\zeta^2}} e^{-\zeta\omega_n (\pi/\omega_d)}\sin(\pi + \phi). \]

Use \( \sin(\pi+\phi)=-\sin(\phi) \) and \( \sin(\phi)=\sqrt{1-\zeta^2} \) to obtain \( y(t_p)=1 + e^{-\zeta\omega_n(\pi/\omega_d)} \). Since \( \omega_d=\omega_n\sqrt{1-\zeta^2} \),

\[ y(t_p)-1 = \exp\!\left(-\frac{\zeta\pi}{\sqrt{1-\zeta^2}}\right) \equiv M_p. \]

Problem 4 (Unity-feedback steady-state step error): For a unity-feedback system with \( G(s)=\frac{K}{s+1} \), compute \( e_{ss} \) for a unit step reference and find \(K\) such that \( e_{ss} \le 0.05 \).

Solution:

The error transfer is \( E(s)=\frac{R(s)}{1+G(s)} \) with \( R(s)=1/s \). Therefore:

\[ e_{ss}=\lim_{s \to 0} sE(s) = \lim_{s \to 0} s\frac{1/s}{1+K/(s+1)} = \frac{1}{1+ \lim_{s \to 0} K/(s+1)} = \frac{1}{1+K}. \]

The constraint \( \frac{1}{1+K} \le 0.05 \) implies \( 1+K \ge 20 \), hence \( K \ge 19 \).

Problem 5 (Dominant pole settling-time estimate): A stable response satisfies \( y(t)=y_\infty + c_1 e^{-2t} + c_2 e^{-6t} \) with \( y_\infty \neq 0 \). Give a conservative bound on the 2% settling time in terms of \( |c_1|, |c_2|, |y_\infty| \).

Solution:

Here \( \sigma = 2 \) and \( C = |c_1|+|c_2| \). The envelope bound gives \( |y(t)-y_\infty| \le C e^{-2t} \). For 2% settling, require \( C e^{-2t} \le 0.02|y_\infty| \), so:

\[ t_s \ge \frac{1}{2}\ln\!\left(\frac{|c_1|+|c_2|}{0.02|y_\infty|}\right). \]

This shows explicitly that the slowest exponential mode dominates the settling-time scaling.

We defined step-response performance metrics rigorously and derived closed-form expressions for first-order and underdamped second-order systems. In particular, overshoot and settling time admit analytic formulas in terms of \( \zeta \) and \( \omega_n \), while rise time often depends on conventions and may require numerical evaluation. For higher-order systems, we formalized dominant-mode intuition via exponential envelope bounds. Finally, we implemented reliable metric extraction pipelines across Python, C++, Java, MATLAB/Simulink, and Mathematica.