In many control systems, especially robotics and mechatronics (e.g., joint servo loops with digital communication and centralized computation), control actions are applied with a non-negligible delay. A pure time delay of length \( L > 0 \) (seconds) is modeled in the Laplace domain as

\[ G_d(s) = e^{-sL}, \quad L > 0. \]

Evaluating on the imaginary axis \( s = j\omega \), we obtain the frequency response

\[ G_d(j\omega) = e^{-j\omega L} = \cos(\omega L) - j \sin(\omega L), \quad \omega \ge 0. \]

Thus a pure delay has:

• Unit magnitude: \( \left|G_d(j\omega)\right| = 1 \) for all \( \omega \).
• Frequency-dependent phase lag: \( \arg G_d(j\omega) = -\omega L \) (radians), which grows linearly in magnitude with frequency.

Consider a standard unity-feedback loop with controller \( C(s) \), plant \( P(s) \), and delay \( L \) in the forward path:

\[ L(s) = C(s)P(s)e^{-sL}, \qquad T(s) = \frac{L(s)}{1+L(s)}. \]

The delay \( e^{-sL} \) rotates the open-loop Nyquist curve clockwise, shrinking the phase margin and therefore reducing robustness. Since actual delays in robotic systems are rarely known exactly, we must analyze robustness to delay uncertainty.

flowchart TD
  A["Model nominal plant Ps and controller Cs"] --> B["Add pure delay model exp(-sLnom)"]
  B --> C["Compute open loop L0(s)"]
  C --> D["Measure phase margin pm and crossover omega_c"]
  D --> E["Compute delay margin Dm = pm / omega_c"]
  E --> F["Compare Dm with expected delay uncertainty"]
  F --> G["Is Dm large enough?"]
  G -->|yes| H["Accept design (robust to delay)"]
  G -->|no| I["Redesign controller (increase pm or reduce omega_c)"]
  

Let \( L_0 \) be a nominal delay and \( \Delta L \) an unknown perturbation:

\[ L = L_0 + \Delta L, \quad |\Delta L| \le \bar{L}. \]

The delay term becomes

\[ e^{-sL} = e^{-sL_0} e^{-s\Delta L}. \]

Grouping \( e^{-sL_0} \) into the nominal model and treating \( e^{-s\Delta L} \) as uncertainty, the open-loop can be written as

\[ L(s) = C(s)P(s)e^{-sL_0} e^{-s\Delta L} = L_0(s)\, e^{-s\Delta L}, \]

where \( L_0(s) = C(s)P(s)e^{-sL_0} \) is the nominal open-loop. In the frequency domain,

\[ L(j\omega) = L_0(j\omega)\, e^{-j\omega\Delta L}, \quad |e^{-j\omega\Delta L}| = 1,\; \arg e^{-j\omega\Delta L} = -\omega\Delta L. \]

Therefore, delay uncertainty appears as a pure phase perturbation whose magnitude grows linearly with frequency and is bounded by \( |\omega\Delta L| \le \omega \bar{L} \). Classical robustness questions become:

• For a given controller and nominal plant, what is the largest admissible \( |\Delta L| \) that preserves stability?
• Conversely, for a given range \( |\Delta L| \le \bar{L} \), how should we design \( C(s) \) to guarantee robust stability?

Consider the nominal open-loop \( L_0(s) = C(s)P(s) \) without delay. Its gain crossover frequency \( \omega_c \) is defined by

\[ |L_0(j\omega_c)| = 1. \]

Let \( \phi_0 = \arg L_0(j\omega_c) \) denote the phase at \( \omega_c \). The (classical) phase margin is

\[ \phi_m = \pi + \phi_0. \]

Now insert a pure delay of length \( D \) in the forward path. At \( \omega_c \), the new open-loop phase becomes

\[ \arg L(j\omega_c) = \arg\left(L_0(j\omega_c)e^{-j\omega_c D}\right) = \phi_0 - \omega_c D. \]

The closed loop reaches the brink of instability when the total phase at crossover is \( -\pi \):

\[ \phi_0 - \omega_c D_{\max} = -\pi. \]

Solving for the critical delay \( D_{\max} \) (sometimes called the delay margin in radians-per-second units):

\[ D_{\max} = \frac{\phi_0 + \pi}{\omega_c} = \frac{\phi_m}{\omega_c}. \]

If \( \phi_m \) is measured in radians, this formula is exact. If \( \phi_m^{\circ} \) is in degrees, then

\[ D_{\max} = \frac{\phi_m^{\circ}\,\pi/180}{\omega_c}. \]

Interpretation:

• The larger the phase margin \( \phi_m \), the more delay we can tolerate at the gain crossover.
• The higher the crossover frequency \( \omega_c \), the smaller the admissible delay.

In many robotic servo designs, one therefore avoids excessively high bandwidth if the communication or computation delays cannot be kept sufficiently small.

Suppose the actual delay is \( L = L_0 + \Delta L \) with \( |\Delta L| \le \bar{L} \). Using the delay margin \( D_{\max} = \phi_m / \omega_c \) computed from the nominal design, a simple sufficient condition for robust stability is

\[ \bar{L} < D_{\max} = \frac{\phi_m}{\omega_c}. \]

This ensures that, even with worst-case delay perturbation, the additional phase lag at \( \omega_c \) does not exhaust the phase margin. More refined analyses examine all unit-circle crossings of the Nyquist plot, but the single-crossover approximation captures the essential idea:

• At each frequency where \( |L_0(j\omega)| = 1 \), delay uncertainty contributes a phase shift of magnitude \( |\omega\Delta L| \).
• If this extra lag does not drive the Nyquist curve through the \( -1 \) point, the loop remains stable.

Design rule of thumb:

• For required robustness to a delay uncertainty bound \( \bar{L} \), choose the controller such that the crossover frequency satisfies \( \omega_c \bar{L} \le \phi_m^{\text{target}} \), where \( \phi_m^{\text{target}} \) includes a safety margin (e.g., \( 10^\circ\text{–}20^\circ \) beyond minimal stability).

Classical tools like root locus and Bode plots are most convenient when the open-loop is a rational transfer function. The exponential term \( e^{-sL} \) makes the characteristic equation transcendental. A standard workaround is the Padé approximation, which replaces \( e^{-sL} \) by a stable, proper rational function.

The first-order Padé approximation of \( e^{-sL} \) is

\[ e^{-sL} \approx \frac{1 - \tfrac{1}{2}sL}{1 + \tfrac{1}{2}sL}. \]

Proof (series matching): Expand both sides around \( s = 0 \). The delay:

\[ e^{-sL} = 1 - sL + \frac{1}{2}(sL)^2 - \frac{1}{6}(sL)^3 + \cdots. \]

The Padé approximant:

\[ \frac{1 - \tfrac{1}{2}sL}{1 + \tfrac{1}{2}sL} = (1 - \tfrac{1}{2}sL)\left(1 - \tfrac{1}{2}sL + \frac{1}{4}(sL)^2 - \cdots\right) = 1 - sL + \frac{1}{2}(sL)^2 + \mathcal{O}\!\left((sL)^3\right). \]

Hence the Taylor coefficients up to order \( (sL)^2 \) match, giving a good low-frequency approximation.

The second-order Padé approximation is

\[ e^{-sL} \approx \frac{1 - \tfrac{1}{2}sL + \tfrac{1}{12}(sL)^2}{1 + \tfrac{1}{2}sL + \tfrac{1}{12}(sL)^2}. \]

Important properties:

• The approximant has poles and zeros symmetric with respect to the imaginary axis, approximating the trajectory of \( e^{-j\omega L} \) on the unit circle.
• Higher-order Padé approximants improve accuracy over a larger frequency range but introduce additional poles and zeros, which can distort robustness estimates if used incautiously.

In classical practice, one often uses a first- or second-order Padé approximation to include delay in Bode and root-locus designs, but validates robustness using the exact exponential model via Nyquist or simulation.

flowchart TD
  A["Exact delay exp(-sL)"] --> B["Choose Padé order N (N=1,2,...)"]
  B --> C["Compute rational approximation H_N(s)"]
  C --> D["Form rational open loop L_N(s)=C(s)P(s)H_N(s)"]
  D --> E["Use root locus / Bode on L_N(s)"]
  E --> F["Validate with exact exp(-sL) model (Nyquist, simulation)"]
        

Suppose a nominal open loop (without delay) has phase margin \( \phi_m^{\circ} = 60^\circ \) at crossover frequency \( \omega_c = 10 \) rad/s. The corresponding delay margin is

\[ D_{\max} = \frac{\phi_m^{\circ}\,\pi/180}{\omega_c} = \frac{60\cdot\pi/180}{10} = \frac{\pi/3}{10} \approx 0.105 \text{ s}. \]

If we expect delay uncertainty \( |\Delta L| \le 0.03 \text{ s} \), the design is comfortably robust. However, if network-induced delays may reach \( 0.12 \text{ s} \), the delay margin is insufficient and we must:

• Increase phase margin (e.g., via lead compensation),
• Lower the crossover frequency (e.g., reducing controller gain), or
• Introduce prediction or feedforward methods (beyond this classical lesson).

In robotic joint control, a typical requirement might be \( \phi_m^{\circ} \approx 50^\circ\text{–}60^\circ \) and \( \omega_c \) chosen such that the delay margin comfortably exceeds all known sources of delay (sensor filtering, communication, OS scheduling jitter, and actuator driver latencies).

We illustrate delay-margin computation using the python-control library for a simple joint-like plant \( P(s) = \frac{1}{s(s+1)} \) with a PD controller. This is a crude model of a lightly damped rotational servo axis.

import numpy as np
import control as ct  # python-control library
import matplotlib.pyplot as plt

# Joint-like plant: P(s) = 1 / (s (s + 1))
s = ct.TransferFunction.s
P = 1 / (s * (s + 1))

# Simple PD controller (chosen arbitrarily for demonstration)
Kp = 40.0
Kd = 2.0
C = Kp + Kd * s

# Nominal open loop without explicit delay term
L0 = C * P

# Classical margins (gain margin, phase margin, crossover frequencies)
gm, pm_deg, wgc, wpc = ct.margin(L0)
print("Phase margin [deg]:", pm_deg)
print("Gain crossover [rad/s]:", wgc)

# Delay margin from phase margin and crossover frequency
pm_rad = pm_deg * np.pi / 180.0
D_max = pm_rad / wgc
print("Approximate delay margin [s]:", D_max)

# Now include an explicit delay via Padé approximation
# Example nominal delay for communication/computation:
L_nom = 0.03  # 30 ms
num_d, den_d = ct.pade(L_nom, 1)  # first-order Padé
D_pade = ct.TransferFunction(num_d, den_d)

L_with_delay = C * P * D_pade
T_with_delay = ct.feedback(L_with_delay, 1)

# Compare step responses for different additional delays
delays = [0.0, 0.02, 0.05, 0.08]  # added to nominal delay
t = np.linspace(0, 5, 1000)

plt.figure()
for d_add in delays:
    # total delay approximated as Padé of (L_nom + d_add)
    num_d2, den_d2 = ct.pade(L_nom + d_add, 1)
    D_tot = ct.TransferFunction(num_d2, den_d2)
    Ld = C * P * D_tot
    T_d = ct.feedback(Ld, 1)
    t_out, y_out = ct.step_response(T_d, T=t)
    label = f"delay = {(L_nom + d_add):.3f} s"
    plt.plot(t_out, y_out, label=label)

plt.xlabel("Time [s]")
plt.ylabel("Joint angle (normalized)")
plt.title("Step response vs. time delay (Padé approximation)")
plt.legend()
plt.grid(True)
plt.show()

# Note: In a robotics stack (e.g. ROS), P and C may come
# from identified models of a DC motor + load, and L_nom
# corresponds to sensing and communication delays.
      

The code uses ct.margin to extract phase margin and gain crossover frequency, then computes the approximate delay margin. Using Padé approximations, we visualize how increasing delay degrades transient performance and eventually destabilizes the loop.

In embedded robotic controllers (e.g., running under a real-time OS or within a ROS node), one may want a lightweight C++ routine to approximate delay margin from a numerically sampled frequency response of a nominal open loop \( L_0(j\omega) \). The following example uses the standard library <complex>. For more complex plants, matrix computations can be handled with Eigen.

#include <iostream>
#include <complex>
#include <cmath>

// Plant P(s) = 1 / (s (s + 1)), controller C(s) = Kp + Kd s
// Evaluate L0(jw) on a log-spaced grid and estimate phase margin and delay margin.

std::complex<double> L0_of_jw(double w, double Kp, double Kd) {
    std::complex<double> j(0.0, 1.0);
    std::complex<double> s = j * w;
    std::complex<double> P = 1.0 / (s * (s + 1.0));
    std::complex<double> C = Kp + Kd * s;
    return C * P;
}

int main() {
    double Kp = 40.0;
    double Kd = 2.0;

    double w_min = 0.1;
    double w_max = 100.0;
    int N = 200;

    double wgc = 0.0;      // estimated gain crossover frequency
    double phi_deg = 0.0;  // phase at crossover

    // Scan over logarithmic grid to find |L0(jw)| ~ 1
    double prev_mag = 0.0;
    double prev_w = w_min;
    for (int k = 0; k <= N; ++k) {
        double alpha = static_cast<double>(k) / N;
        double w = w_min * std::pow(w_max / w_min, alpha);
        std::complex<double> L = L0_of_jw(w, Kp, Kd);
        double mag = std::abs(L);

        if (k > 0) {
            // Look for a sign change in (mag - 1), indicating a crossing of unity gain
            double f_prev = prev_mag - 1.0;
            double f_curr = mag - 1.0;
            if (f_prev * f_curr <= 0.0) {
                // Linear interpolation for a better estimate of wgc
                double t = f_prev / (f_prev - f_curr + 1e-12);
                wgc = prev_w + t * (w - prev_w);
                break;
            }
        }
        prev_mag = mag;
        prev_w = w;
    }

    if (wgc > 0.0) {
        std::complex<double> Lc = L0_of_jw(wgc, Kp, Kd);
        phi_deg = std::arg(Lc) * 180.0 / M_PI;
        double pm_deg = 180.0 + phi_deg;     // phase margin in degrees
        double pm_rad = pm_deg * M_PI / 180.0;
        double D_max = pm_rad / wgc;         // delay margin in seconds

        std::cout << "Gain crossover wgc = " << wgc << " rad/s\n";
        std::cout << "Phase at wgc = " << phi_deg << " deg\n";
        std::cout << "Phase margin = " << pm_deg << " deg\n";
        std::cout << "Approximate delay margin D_max = "
                  << D_max << " s\n";
    } else {
        std::cout << "No gain crossover found in the scanned range.\n";
    }

    return 0;
}
      

This code can be integrated into a robotic controller configuration tool that checks whether measured or specified worst-case delays (sensing, control computation, actuator drivers, and communication) are below the computed delay margin.

In some Java-based robotics frameworks (e.g., FIRST robotics, custom Java middleware), classical analysis is done offline (Python/MATLAB), and only scalar robustness metrics such as phase margin, crossover frequency, and delay margin are carried into the real-time controller. The following code illustrates a simple check that disables aggressive commands when the estimated delay exceeds a chosen fraction of the delay margin.

public class DelayRobustnessMonitor {

    // Phase margin and crossover frequency obtained from offline Bode analysis
    private final double phaseMarginDeg; // e.g., 55.0
    private final double omegaC;         // rad/s, e.g., 12.0
    private final double delayMargin;    // seconds

    // Safety factor: maximum allowed delay as a fraction of delay margin
    private final double safetyFactor;   // e.g., 0.6 (60%)

    public DelayRobustnessMonitor(double phaseMarginDeg, double omegaC,
                                  double safetyFactor) {
        this.phaseMarginDeg = phaseMarginDeg;
        this.omegaC = omegaC;
        this.safetyFactor = safetyFactor;

        double phaseMarginRad = Math.toRadians(phaseMarginDeg);
        this.delayMargin = phaseMarginRad / omegaC;
    }

    public double getDelayMargin() {
        return delayMargin;
    }

    public boolean isDelayAcceptable(double estimatedDelay) {
        // estimatedDelay: current worst-case I/O + computation delay (seconds)
        double maxAllowed = safetyFactor * delayMargin;
        return estimatedDelay <= maxAllowed;
    }

    public static void main(String[] args) {
        DelayRobustnessMonitor monitor =
            new DelayRobustnessMonitor(55.0, 12.0, 0.6);

        double estimatedDelay = 0.020; // 20 ms, for example
        System.out.println("Delay margin [s] = " + monitor.getDelayMargin());
        if (monitor.isDelayAcceptable(estimatedDelay)) {
            System.out.println("Delay is acceptable for high-performance mode.");
        } else {
            System.out.println("Delay too large: switch to conservative gains.");
        }
    }
}
      

Such logic can be combined with runtime delay estimators and robot state monitors to automatically adapt controller aggressiveness to current network and computation conditions while retaining classical robustness guarantees.

10.1 MATLAB / Simulink

MATLAB® Control System Toolbox provides direct computation of delay margins via margin and allmargin. For the same joint-like example:

% Joint-like plant and PD controller
s = tf('s');
P = 1 / (s * (s + 1));
Kp = 40;
Kd = 2;
C = Kp + Kd * s;

L0 = C * P;

% Classical stability margins
[gm, pm_deg, wgc, wpc] = margin(L0);
pm_rad = pm_deg * pi / 180;
D_max = pm_rad / wgc;

fprintf('Phase margin = %.2f deg\n', pm_deg);
fprintf('Gain crossover wgc = %.2f rad/s\n', wgc);
fprintf('Approximate delay margin D_max = %.4f s\n', D_max);

% Use first-order Padé approximation for a nominal delay
L_nom = 0.03;  % 30 ms
[num_d, den_d] = pade(L_nom, 1);
D_pade = tf(num_d, den_d);
L = C * P * D_pade;
T = feedback(L, 1);

figure;
step(T);
title('Step response with 30 ms nominal delay (Padé approx)');
grid on;

% Simulink: implement P, C, and a Transport Delay block:
%  - Use "Transfer Fcn" blocks for P and C
%  - Insert a "Transport Delay" block with DelayTime = L_nom
%  - Close unity feedback and simulate response
      

10.2 Wolfram Mathematica

In Mathematica, we can similarly compute delay margins and visualize how Padé approximations affect the frequency response.

(* Joint-like plant and PD controller *)
s = ComplexExpand[I*ω] /. ω -> s; (* symbolic placeholder, not used directly *)

P[s_] := 1/(s (s + 1));
C[s_] := 40 + 2 s;

(* Nominal open loop *)
L0[s_] := C[s] P[s];

(* Use built-in frequency response tools *)
ωlist = Table[10.^x, {x, -1, 2, 0.01}];
resp = Table[{ω, L0[I ω]}, {ω, ωlist}];

(* Extract magnitude and phase (in degrees) *)
magPhase = Table[
   {ω, Abs[L0[I ω]], Arg[L0[I ω]] 180./Pi},
   {ω, ωlist}
];

(* Find crossover frequency approximately where magnitude crosses 1 *)
crossIndex = First@FirstPosition[
   Partition[magPhase[[All, 2]] - 1, 2, 1],
   {a_, b_} /; a b <= 0
];
ωc = ωlist[[crossIndex]];

φ0deg = Last@Select[magPhase, #[[1]] == ωc &][[1]];
φmdeg = 180. + φ0deg;
Dm = (φmdeg Pi/180.)/ωc;

Print["Phase margin [deg] = ", φmdeg];
Print["Gain crossover [rad/s] = ", ωc];
Print["Approximate delay margin [s] = ", Dm];

(* First-order Padé approximation of exp(-s L) *)
Lnom = 0.03;
Hd[s_] := (1 - (Lnom s)/2)/(1 + (Lnom s)/2);
L[s_] := L0[s] Hd[s];

(* BodePlot of approximate delayed loop *)
BodePlot[{L0[s] /. s -> I ω, L[s] /. s -> I ω},
 {ω, 0.1, 100},
 PlotLegends -> {"No delay", "Padé delay"},
 GridLines -> Automatic]
      

These scripts demonstrate how classical delay-margin concepts are embedded in software tools used routinely in robotics and mechatronics for loop-shaping and robustness validation.

Problem 1 (Derivation of Delay Margin): Let \( L_0(s) = C(s)P(s) \) be the nominal open-loop transfer function (no explicit delay). Its gain crossover frequency \( \omega_c \) is defined by \( |L_0(j\omega_c)| = 1 \), and its phase at \( \omega_c \) is \( \phi_0 = \arg L_0(j\omega_c) \). Show that the maximum tolerable pure delay \( D_{\max} \) in the forward path, before losing stability, satisfies \( D_{\max} = \phi_m / \omega_c \), where \( \phi_m \) is the phase margin in radians.

Solution:

By definition, the phase margin is

\[ \phi_m = \pi + \phi_0. \]

Adding a delay \( D \) multiplies the open-loop by \( e^{-j\omega D} \), so the phase at the gain crossover becomes

\[ \arg L(j\omega_c) = \phi_0 - \omega_c D. \]

At the stability boundary, the phase at gain crossover must be \( -\pi \):

\[ \phi_0 - \omega_c D_{\max} = -\pi. \]

Solving for \( D_{\max} \):

\[ D_{\max} = \frac{\phi_0 + \pi}{\omega_c} = \frac{\phi_m}{\omega_c}. \]

This exactly matches the relationship used in classical delay-margin calculations.

Problem 2 (Numerical Delay Margin and Robustness Check): A feedback system has nominal phase margin \( \phi_m^{\circ} = 50^\circ \) at crossover frequency \( \omega_c = 5 \) rad/s. The total delay (sensing + computation + actuation) is uncertain but satisfies \( |\Delta L| \le 0.05 \) s. Is the system robustly stable against this delay uncertainty according to the classical delay-margin approximation?

Solution:

Compute the delay margin:

\[ D_{\max} = \frac{\phi_m^{\circ}\,\pi/180}{\omega_c} = \frac{50\cdot\pi/180}{5} = \frac{5\pi}{90} \approx 0.1745 \text{ s}. \]

Since \( \bar{L} = 0.05 \text{ s} \) and \( \bar{L} < D_{\max} \), the system is predicted to be robustly stable with respect to this delay uncertainty. In practice, some additional safety margin is often added.

Problem 3 (Padé Approximation Accuracy Near the Origin): Show that the first-order Padé approximation \( e^{-sL} \approx (1 - \tfrac{1}{2}sL)/(1 + \tfrac{1}{2}sL) \) matches the Taylor series of \( e^{-sL} \) up to terms of order \( (sL)^2 \).

Solution:

The Taylor series of the exact delay is

\[ e^{-sL} = 1 - sL + \frac{1}{2}(sL)^2 - \frac{1}{6}(sL)^3 + \cdots. \]

The Padé approximant is

\[ \frac{1 - \tfrac{1}{2}sL}{1 + \tfrac{1}{2}sL} = (1 - \tfrac{1}{2}sL)\bigl(1 - \tfrac{1}{2}sL + \tfrac{1}{4}(sL)^2 - \cdots\bigr) = 1 - sL + \frac{1}{2}(sL)^2 + \mathcal{O}\!\left((sL)^3\right). \]

Thus, coefficients of \( 1 \), \( sL \), and \( (sL)^2 \) are identical in the two series, so the approximation is second-order accurate around \( s = 0 \).

Problem 4 (Crossover Frequency Constraint for Delay Robustness): You must design a controller for a robotic axis subject to an uncertain delay with bound \( \bar{L} = 0.02 \) s. You require a minimal safety factor of \( 2 \), i.e., \( 2\bar{L} \le D_{\max} \). Assuming you will design for a phase margin of \( \phi_m^{\circ} = 60^\circ \), derive an upper bound on the allowable crossover frequency \( \omega_c \).

Solution:

The requirement is

\[ 2\bar{L} \le D_{\max} = \frac{\phi_m^{\circ}\,\pi/180}{\omega_c}. \]

Rearranging for \( \omega_c \) gives

\[ \omega_c \le \frac{\phi_m^{\circ}\,\pi/180}{2\bar{L}}. \]

Plugging in \( \phi_m^{\circ} = 60^\circ \) and \( \bar{L} = 0.02 \) s:

\[ \omega_c \le \frac{60\cdot\pi/180}{2\cdot 0.02} = \frac{\pi/3}{0.04} \approx 26.18 \text{ rad/s}. \]

Thus, the controller should be designed such that the crossover frequency does not exceed approximately \( 26 \) rad/s to ensure the desired robustness to delay uncertainty.

In this lesson we analyzed the impact of pure time delays and delay uncertainty on classical SISO feedback systems. We showed that delays contribute a frequency-dependent phase lag while leaving the magnitude unchanged, and derived the fundamental relation between phase margin and delay margin, \( D_{\max} = \phi_m / \omega_c \). This provides a simple, powerful rule for ensuring robust stability in the presence of uncertain delays.

We introduced Padé approximations as rational surrogates for delay, enabling classical root-locus and Bode-based design while emphasizing the need to validate robustness using exact exponential delay models. Finally, we demonstrated how to compute and exploit delay margins with Python, C++, Java, MATLAB/Simulink, and Mathematica, particularly in the context of robotic servo systems where delays are unavoidable. These ideas form a crucial bridge between classical robustness analysis and more advanced treatments of time-delay systems.