In many engineered systems, information, mass, or energy does not appear at an output immediately after it is applied at an input. Instead, it must travel through space: along a pipe, conveyor, cable, chemical reactor, thermal duct, or fluid channel. If the dominant phenomenon is convection/propagation (with negligible diffusion and storage), the input is reproduced at the output after a fixed time shift \( T \). This is a transport delay (also called dead time).

The ideal delay operator maps input \( u(t) \) to output \( y(t) \) by \( y(t)=u(t-T) \). In the Laplace domain (zero initial conditions), \( Y(s)=e^{-sT}U(s) \). The factor \( e^{-sT} \) is not a rational function of \( s \), which means it cannot be represented exactly by a finite-dimensional LTI state-space model.

flowchart TD
  A["Physical transport (pipe / belt / signal path)"] --> B["Distributed model (transport PDE)"]
  B --> C["Ideal delay: y(t)=u(t-T)"]
  C --> D["Laplace form: Gd(s)=exp(-sT) (non-rational)"]
  D --> E["Need finite model for analysis/simulation"]
  E --> F["Rational approximation (Pade)"]
  E --> G["Stable pseudo-distributed \n(N lags / method-of-lines)"]
  F --> H["Good low-frequency phase match; \nmay introduce RHP zeros"]
  G --> I["Minimum-phase, stable; \nimproves as N increases"]
        

Consider a scalar field \( q(x,t) \) transported with constant velocity \( v>0 \) along a 1D domain \( x\in[0,L] \) with negligible diffusion and reaction. The simplest transport PDE is

\[ \frac{\partial q}{\partial t}(x,t) + v\,\frac{\partial q}{\partial x}(x,t) = 0,\quad 0 \le x \le L,\; t \ge 0. \]

Impose an inlet boundary condition \( q(0,t)=u(t) \) and define the measured outlet as \( y(t)=q(L,t) \). We solve the PDE using characteristics.

Characteristic curves. Let \( x(t) \) satisfy \( \dot{x}(t)=v \). Along such a curve, the total derivative of \( q(x(t),t) \) is

\[ \frac{d}{dt}q(x(t),t) = \frac{\partial q}{\partial t} + \dot{x}\,\frac{\partial q}{\partial x} = \frac{\partial q}{\partial t} + v\,\frac{\partial q}{\partial x} = 0. \]

Hence \( q \) is constant along characteristics. Since \( \dot{x}=v \), we have \( x(t)=vt + c \), or equivalently the invariant \( t - x/v = \text{constant} \). Therefore, for some function \( \phi \),

\[ q(x,t)=\phi\!\left(t-\frac{x}{v}\right). \]

Use the inlet boundary condition: \( q(0,t)=\phi(t)=u(t) \), so \( \phi(\cdot)=u(\cdot) \). Thus the exact solution is

\[ q(x,t)=u\!\left(t-\frac{x}{v}\right),\qquad y(t)=q(L,t)=u\!\left(t-\frac{L}{v}\right). \]

Define the transport time \( T = \frac{L}{v} \). Then the PDE input-output map is exactly the delay: \( y(t)=u(t-T) \).

Taking Laplace transforms (zero initial conditions) yields the canonical transfer function \( G_d(s)=\frac{Y(s)}{U(s)}=e^{-sT} \).

\[ Y(s) = \int_0^\infty u(t-T)e^{-st}\,dt = e^{-sT}\int_0^\infty u(\tau)e^{-s\tau}\,d\tau = e^{-sT}U(s). \]

Many tools in system dynamics and classical control assume rational transfer functions \( G(s)=\frac{N(s)}{D(s)} \) (finite-degree polynomials), corresponding to finite-dimensional ODE state-space models. The delay element \( e^{-sT} \) is transcendental and leads to characteristic equations of the form \( D(s) + N(s)e^{-sT}=0 \), which generally have infinitely many roots.

A key frequency-domain identity is that the ideal delay has unit magnitude and linear phase:

\[ e^{-j\omega T} = \cos(\omega T) - j\sin(\omega T),\qquad \left|e^{-j\omega T}\right|=1,\qquad \angle e^{-j\omega T} = -\omega T. \]

So delay does not attenuate amplitude, but it adds phase lag proportional to frequency, which can strongly reduce stability margins when delay is inside feedback loops (detailed loop-shaping is treated in later chapters).

A common way to obtain a finite-dimensional approximation is to replace \( e^{-sT} \) by a rational function that matches the Taylor series of the exponential around \( s=0 \). The \([n/n]\) Padé approximant is defined as the rational function \( \frac{P_n(s)}{Q_n(s)} \) such that

\[ e^{-sT} - \frac{P_n(s)}{Q_n(s)} = \mathcal{O}\!\left(s^{2n+1}\right)\quad\text{as}\quad s \to 0. \]

This means the first \( 2n \) derivatives at \( s=0 \) match exactly (maximal local accuracy). Two low-order cases are especially important:

(a) First-order Padé \([1/1]\):

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

(b) Second-order Padé \([2/2]\):

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

Important structural caution. The \([1/1]\) approximation has a zero at \( s=\frac{2}{T} \), i.e., in the right half-plane, so it is non-minimum-phase even though the true delay is all-pass in frequency. Higher-order Padé approximations can also introduce right-half-plane zeros. This matters in feedback design and in time-domain simulation (non-minimum-phase approximations can generate undershoot/overshoot artifacts).

Padé matching proof sketch. Let the Maclaurin series be \( e^{-sT}=\sum_{k=0}^{\infty}\frac{(-T)^k}{k!}s^k \). Padé constructs \( P_n,Q_n \) so that the series of \( Q_n(s)e^{-sT}-P_n(s) \) has zero coefficients up to degree \( 2n \), yielding the stated remainder order \( \mathcal{O}(s^{2n+1}) \).

A different approximation philosophy is to preserve physical realizability and stability by using a cascade of stable low-order elements. The most widely used form is the N-lag chain:

\[ G_N(s) \;=\; \left(\frac{1}{1 + \frac{sT}{N}}\right)^N. \]

This is always stable and minimum-phase for any finite \( N \) because it is a product of stable first-order factors. It is called pseudo-distributed because it mimics a distributed transport line by discretizing space into \( N \) compartments.

Limit theorem (convergence to the true delay). For any complex \( s \),

\[ \lim_{N\to\infty}\left(1+\frac{sT}{N}\right)^N = e^{sT} \quad\Longrightarrow\quad \lim_{N\to\infty}\left(1+\frac{sT}{N}\right)^{-N} = e^{-sT}. \]

Proof. Use the standard exponential limit \( \lim_{N\to\infty}\left(1+\frac{z}{N}\right)^N=e^z \) with \( z=sT \), then invert both sides. This shows that increasing \( N \) improves fidelity, especially at low-to-moderate frequencies.

State-space realization (continuous-time chain). Let \( x_k(t) \) be the state of compartment \( k \). With time constant \( \tau = T/N \), a standard chain is:

\[ \dot{x}_1 = \frac{1}{\tau}(u - x_1),\qquad \dot{x}_k = \frac{1}{\tau}(x_{k-1} - x_k),\; k=2,\dots,N,\qquad y = x_N. \]

Taking Laplace transforms gives \( X_1/U = \frac{1}{1+s\tau} \) and recursively \( X_k/U = \left(\frac{1}{1+s\tau}\right)^k \), hence \( Y/U = G_N(s) \).

flowchart TD
  U["u(t)"] --> X1["Lag 1: (1/(1+s*T/N))"]
  X1 --> X2["Lag 2"]
  X2 --> X3["..."]
  X3 --> XN["Lag N"]
  XN --> Y["y(t) approx u(t-T)"]
  NOTE["As N increases: \nbetter delay fidelity (still stable)"] --> Y
        

Connection to spatial discretization. If you discretize the transport PDE in space (method-of-lines), you obtain a large ODE system. Many stable discretizations introduce numerical diffusion, and the resulting finite model resembles a chain that approximates a pure shift. Thus, lag chains are a principled bridge between distributed PDEs and lumped ODE models.

Practical guideline: Padé approximations are attractive when you need a compact rational model that matches phase well near \( \omega \approx 0 \). Lag chains are attractive when you need guaranteed stability/minimum-phase behavior (e.g., robust simulation blocks), accepting that larger \( N \) may be required.

In sampled-data simulation, an exact integer delay of \( k \) samples is simply \( y[n]=u[n-k] \), which corresponds to \( G(z)=z^{-k} \). This is implemented exactly by a ring buffer.

When a delay is not an integer multiple of the sample time (fractional delay), one uses interpolation or IIR all-pass designs whose group delay is approximately flat at low frequencies. A classical closed-form approach is the maximally-flat group-delay all-pole/all-pass construction introduced by Thiran (details are treated more deeply in digital modeling chapters).

In this lesson’s code, the C++/Java/MATLAB examples demonstrate: (i) exact integer-sample delay via ring buffer, and (ii) a continuous-time lag-chain approximation integrated in discrete time.

The following script constructs Padé and lag-chain approximations for \( e^{-sT} \), compares phase, and simulates a first-order plant \( G_p(s)=\frac{1}{s+1} \) driven by a delayed unit step.

# Chapter19_Lesson4.py
"""System Dynamics — Chapter 19, Lesson 4
Transport delay models and pseudo-distributed approximations.

This script demonstrates:
  (1) Pure delay: y(t) = u(t - T)  <=>  G(s) = exp(-s T)
  (2) Padé rational approximations of exp(-sT)
  (3) Stable 'pseudo-distributed' approximation via an N-lag chain:
          G_N(s) = (1 / (1 + s T/N))^N
  (4) Time-domain simulation of a first-order plant driven by delayed input.

Dependencies: numpy, scipy, matplotlib
"""

import numpy as np
from math import ceil

try:
    from scipy import signal
except Exception as e:
    raise RuntimeError("This script requires SciPy (scipy.signal).") from e

import matplotlib.pyplot as plt


def chain_delay_tf(T: float, N: int):
    """Return (num, den) of (1/(1+sT/N))^N."""
    if T <= 0:
        raise ValueError("T must be positive.")
    if N < 1:
        raise ValueError("N must be >= 1.")
    num = np.array([1.0])
    den = np.array([T / N, 1.0])
    denN = den.copy()
    for _ in range(N - 1):
        denN = np.convolve(denN, den)
    return num, denN


def pade_delay_tf(T: float, n: int):
    """Return (num, den) for n-th order Padé approximation of exp(-sT)."""
    if T <= 0:
        raise ValueError("T must be positive.")
    if n < 1:
        raise ValueError("n must be >= 1.")
    # SciPy provides a standard Padé routine for delays.
    num, den = signal.pade(T, n)
    return np.array(num, dtype=float), np.array(den, dtype=float)


def series_tf(num1, den1, num2, den2):
    """Multiply two transfer functions."""
    return np.convolve(num1, num2), np.convolve(den1, den2)


def step_response(num, den, t):
    sys = signal.TransferFunction(num, den)
    tout, y = signal.step(sys, T=t)
    return tout, y


def freq_response(num, den, w):
    sys = signal.TransferFunction(num, den)
    w, H = signal.freqresp(sys, w=w)
    return w, H


def simulate_true_delay_first_order(T=1.0, a=1.0, b=1.0, dt=1e-3, t_end=12.0):
    """Simulate ydot = -a y + b u(t-T) with u(t)=1 (unit step) using a buffer."""
    if dt <= 0:
        raise ValueError("dt must be positive.")
    n_steps = int(ceil(t_end / dt)) + 1
    delay_steps = int(round(T / dt))
    buf = np.zeros(delay_steps + 1)
    y = 0.0
    ys = np.zeros(n_steps)
    ts = np.arange(n_steps) * dt

    for k in range(n_steps):
        u = 1.0
        u_del = buf[0]
        buf[:-1] = buf[1:]
        buf[-1] = u
        y = y + dt * (-a * y + b * u_del)
        ys[k] = y

    return ts, ys


def main():
    T = 1.0
    n_pade = 2
    N_chain = 8

    # Plant: Gp(s) = 1/(s+1)
    num_p = np.array([1.0])
    den_p = np.array([1.0, 1.0])

    num_pade, den_pade = pade_delay_tf(T, n_pade)
    num_chain, den_chain = chain_delay_tf(T, N_chain)

    num_pd, den_pd = series_tf(num_p, den_p, num_pade, den_pade)
    num_cd, den_cd = series_tf(num_p, den_p, num_chain, den_chain)

    t = np.linspace(0, 12, 2000)
    t1, y_pade = step_response(num_pd, den_pd, t)
    t2, y_chain = step_response(num_cd, den_cd, t)

    t3, y_true = simulate_true_delay_first_order(T=T, a=1.0, b=1.0, dt=2e-3, t_end=12.0)

    plt.figure()
    plt.plot(t3, y_true, label="True delay (buffer sim)")
    plt.plot(t1, y_pade, label=f"Padé [{n_pade}/{n_pade}]")
    plt.plot(t2, y_chain, label=f"Chain N={N_chain}")
    plt.xlabel("t [s]")
    plt.ylabel("y(t)")
    plt.title("First-order plant with transport delay: step response")
    plt.grid(True)
    plt.legend()

    w = np.logspace(-2, 2, 800)
    _, H_pade = freq_response(num_pade, den_pade, w)
    _, H_chain = freq_response(num_chain, den_chain, w)
    H_ideal = np.exp(-1j * w * T)

    plt.figure()
    plt.semilogx(w, np.unwrap(np.angle(H_ideal)), label="Ideal delay phase")
    plt.semilogx(w, np.unwrap(np.angle(H_pade)), label=f"Padé [{n_pade}/{n_pade}] phase")
    plt.semilogx(w, np.unwrap(np.angle(H_chain)), label=f"Chain N={N_chain} phase")
    plt.xlabel("ω [rad/s]")
    plt.ylabel("Phase [rad]")
    plt.title("Delay phase: ideal vs approximations")
    plt.grid(True)
    plt.legend()

    plt.show()


if __name__ == "__main__":
    main()
      

// Chapter19_Lesson4.cpp
// System Dynamics — Chapter 19, Lesson 4
// Discrete-time simulation of a transport delay using:
//   (1) a ring-buffer delay block (exact integer-sample delay)
//   (2) a stable pseudo-distributed N-lag chain (approximation of e^{-sT})
//
// Build (example):
//   g++ -O2 -std=c++17 Chapter19_Lesson4.cpp -o lesson4
//
// Run:
//   ./lesson4

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>

class DelayLine {
public:
    explicit DelayLine(int delaySamples)
        : buf_(std::max(1, delaySamples + 1), 0.0), idx_(0) {}

    double step(double u) {
        double y = buf_[idx_];
        buf_[idx_] = u;
        idx_ = (idx_ + 1) % static_cast<int>(buf_.size());
        return y;
    }

private:
    std::vector<double> buf_;
    int idx_;
};

class LagChain {
public:
    LagChain(double T, int N)
        : T_(T), N_(N), x_(std::max(1, N), 0.0) {
        if (T_ <= 0.0) throw std::runtime_error("T must be positive.");
        if (N_ < 1) throw std::runtime_error("N must be >= 1.");
    }

    double step(double u, double dt) {
        double alpha = (static_cast<double>(N_) / T_);
        std::vector<double> xnew = x_;
        xnew[0] = x_[0] + dt * alpha * (u - x_[0]);
        for (int k = 1; k < N_; ++k) {
            xnew[k] = x_[k] + dt * alpha * (x_[k - 1] - x_[k]);
        }
        x_.swap(xnew);
        return x_.back();
    }

private:
    double T_;
    int N_;
    std::vector<double> x_;
};

int main() {
    const double dt = 0.001;
    const double tEnd = 8.0;

    const double a = 1.0, b = 1.0;

    const double T = 1.0;
    const int delaySamples = static_cast<int>(std::round(T / dt));

    DelayLine delay(delaySamples);
    LagChain chain(T, 10);

    double y_delay = 0.0;
    double y_chain = 0.0;

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "# t, y_trueDelay, y_chainApprox\\n";

    for (int k = 0; k <= static_cast<int>(tEnd / dt); ++k) {
        double t = k * dt;
        double u = (t >= 0.0) ? 1.0 : 0.0;

        double u_del = delay.step(u);
        double u_chain = chain.step(u, dt);

        y_delay = y_delay + dt * (-a * y_delay + b * u_del);
        y_chain = y_chain + dt * (-a * y_chain + b * u_chain);

        if (k % 50 == 0) {
            std::cout << t << ", " << y_delay << ", " << y_chain << "\\n";
        }
    }

    std::cout << "\\n# Note: The chain approximation is stable and minimum-phase,\\n"
              << "# but it does not reproduce an exact dead-time unless N→∞.\\n";
    return 0;
}
      

// Chapter19_Lesson4.java
// System Dynamics — Chapter 19, Lesson 4
// Discrete-time delay line (ring buffer) + N-lag chain approximation.

import java.util.Locale;

class DelayLine {
    private final double[] buf;
    private int idx;

    public DelayLine(int delaySamples) {
        int n = Math.max(1, delaySamples + 1);
        buf = new double[n];
        idx = 0;
    }

    public double step(double u) {
        double y = buf[idx];
        buf[idx] = u;
        idx = (idx + 1) % buf.length;
        return y;
    }
}

class LagChain {
    private final double T;
    private final int N;
    private final double[] x;

    public LagChain(double T, int N) {
        if (T <= 0.0) throw new IllegalArgumentException("T must be positive.");
        if (N < 1) throw new IllegalArgumentException("N must be >= 1.");
        this.T = T;
        this.N = N;
        this.x = new double[N];
    }

    public double step(double u, double dt) {
        double alpha = ((double) N) / T;
        double[] xnew = x.clone();
        xnew[0] = x[0] + dt * alpha * (u - x[0]);
        for (int k = 1; k < N; k++) {
            xnew[k] = x[k] + dt * alpha * (x[k - 1] - x[k]);
        }
        System.arraycopy(xnew, 0, x, 0, N);
        return x[N - 1];
    }
}

public class Chapter19_Lesson4 {
    public static void main(String[] args) {
        Locale.setDefault(Locale.US);

        double dt = 0.001;
        double tEnd = 8.0;

        double a = 1.0, b = 1.0;

        double T = 1.0;
        int delaySamples = (int) Math.round(T / dt);

        DelayLine delay = new DelayLine(delaySamples);
        LagChain chain = new LagChain(T, 10);

        double yDelay = 0.0;
        double yChain = 0.0;

        System.out.println("# t, y_trueDelay, y_chainApprox");
        int steps = (int) Math.round(tEnd / dt);

        for (int k = 0; k <= steps; k++) {
            double t = k * dt;
            double u = (t >= 0.0) ? 1.0 : 0.0;

            double uDel = delay.step(u);
            double uChain = chain.step(u, dt);

            yDelay = yDelay + dt * (-a * yDelay + b * uDel);
            yChain = yChain + dt * (-a * yChain + b * uChain);

            if (k % 50 == 0) {
                System.out.printf("%.6f, %.6f, %.6f%n", t, yDelay, yChain);
            }
        }

        System.out.println("\\n# The delay line is an exact integer-sample delay.");
        System.out.println("# The chain is a stable approximation of a transport delay.");
    }
}
      

% Chapter19_Lesson4.m
% System Dynamics — Chapter 19, Lesson 4
% Transport delay models and pseudo-distributed approximations.
%
% Demonstrates:
%  (1) Gd(s) = exp(-sT)
%  (2) Padé approximation (built-in)
%  (3) Stable N-lag chain approximation: (1/(1+sT/N))^N
%  (4) Step response of Gp(s)=1/(s+1) with delay.

clear; clc;

T = 1.0;          % delay [s]
nPade = 2;        % Padé order
Nchain = 8;       % chain order

s = tf('s');
Gp = 1/(s+1);

% Padé approximation of exp(-sT)
[numP, denP] = pade(T, nPade);
GdPade = tf(numP, denP);

% Stable pseudo-distributed chain approximation
GdChain = 1;
for k = 1:Nchain
    GdChain = series(GdChain, 1/(1 + s*(T/Nchain)));
end

% Combine with plant
G_pade  = series(Gp, GdPade);
G_chain = series(Gp, GdChain);

% Ideal delay (MATLAB can represent dead-time using InputDelay)
G_ideal = Gp;
G_ideal.InputDelay = T;

t = linspace(0, 12, 2000);

figure;
step(G_ideal, G_pade, G_chain, t);
grid on;
legend('Ideal delay','Padé','Chain');
title('First-order plant with transport delay: step response');

% Frequency-domain: delay phase
w = logspace(-2, 2, 800);
[magI, phI] = bode(tf(1,1,'InputDelay',T), w); %#ok<ASGLU>
[magP, phP] = bode(GdPade, w); %#ok<ASGLU>
[magC, phC] = bode(GdChain, w); %#ok<ASGLU>

figure;
semilogx(w, squeeze(phI), w, squeeze(phP), w, squeeze(phC));
grid on;
xlabel('\omega [rad/s]'); ylabel('Phase [deg]');
legend('Ideal delay','Padé','Chain');
title('Delay phase: ideal vs approximations');

% Discrete-time exact delay via buffer (optional)
dt = 0.001;
K = round(T/dt);
buf = zeros(K+1,1);
y = 0;
a = 1; b = 1;
tEnd = 8;
steps = round(tEnd/dt);

for k = 1:steps
    u = 1;                 % step
    uDel = buf(1);
    buf(1:end-1) = buf(2:end);
    buf(end) = u;
    y = y + dt*(-a*y + b*uDel);
end
disp('Done.');
      

Simulink note (conceptual): The ideal dead-time is represented by the Transport Delay block. The stable pseudo-distributed model can be built by cascading Transfer Fcn blocks of the form \( \frac{1}{1+s(T/N)} \) repeated \( N \) times.

(* Chapter19_Lesson4.nb
   System Dynamics — Chapter 19, Lesson 4
   Transport delay models and pseudo-distributed approximations.
*)

ClearAll["Global`*"];

T = 1.0;
n = 2;
Nchain = 8;

(* Padé approximation of Exp[-s T] about s=0 *)
pade = PadeApproximant[Exp[-s T], {s, 0, n}];

(* Stable chain approximation: (1/(1+sT/N))^N *)
chain = (1/(1 + s T/Nchain))^Nchain;

(* Build transfer function models (delay approximations only) *)
Gp = TransferFunctionModel[1/(s + 1), s];
GdPade = TransferFunctionModel[pade, s];
GdChain = TransferFunctionModel[chain, s];

G_pade = SystemsModelSeriesConnect[Gp, GdPade];
G_chain = SystemsModelSeriesConnect[Gp, GdChain];

(* Step responses *)
tmax = 12;
p1 = OutputResponse[G_pade, UnitStep[t], {t, 0, tmax}];
p2 = OutputResponse[G_chain, UnitStep[t], {t, 0, tmax}];

Plot[{p1, p2}, {t, 0, tmax},
 PlotLegends -> {"Padé [2/2]", "Chain N=8"},
 GridLines -> Automatic,
 PlotRange -> All,
 AxesLabel -> {"t", "y(t)"}]

(* Frequency-domain phase comparison (delay only) *)
wmin = 10^-2; wmax = 10^2;
phIdeal[w_] := Arg[Exp[-I w T]];
phPade[w_] := Arg[Evaluate[pade /. s -> I w]];
phChain[w_] := Arg[Evaluate[chain /. s -> I w]];

LogLinearPlot[
 {phIdeal[w], phPade[w], phChain[w]},
 {w, wmin, wmax},
 PlotLegends -> {"Ideal", "Padé [2/2]", "Chain N=8"},
 GridLines -> Automatic,
 AxesLabel -> {"ω", "phase (rad)"}
]
      

Problem 1 (Transport PDE → Delay): Solve \( q_t + v q_x = 0 \) on \( x\in[0,L] \) with inlet \( q(0,t)=u(t) \). Show that the outlet \( y(t)=q(L,t) \) equals \( u(t-T) \) with \( T=L/v \).

Solution: Along characteristics \( \dot{x}=v \),

\[ \frac{d}{dt}q(x(t),t) = q_t + \dot{x} q_x = q_t + v q_x = 0, \]

hence \( q \) is constant on curves \( t-x/v=\text{const} \). Therefore \( q(x,t)=u(t-x/v) \), and \( y(t)=q(L,t)=u(t-L/v)=u(t-T) \).

Problem 2 (Derive Padé \([1/1]\)): Find coefficients \( a_0,a_1,b_0,b_1 \) such that \( e^{-sT}\approx \frac{a_0+a_1 s}{b_0+b_1 s} \) matches the Taylor series up to order \( s^2 \).

Solution: Normalize \( b_0=1 \). Expand:

\[ \frac{a_0+a_1 s}{1+b_1 s} = (a_0+a_1 s)\left(1 - b_1 s + b_1^2 s^2 + \mathcal{O}(s^3)\right) = a_0 + (a_1-a_0 b_1)s + (a_0 b_1^2 - a_1 b_1)s^2 + \mathcal{O}(s^3). \]

Match to \( e^{-sT}=1 - Ts + \frac{T^2}{2}s^2 + \mathcal{O}(s^3) \): \( a_0=1 \), \( a_1-b_1=-T \), \( b_1^2-a_1 b_1=\frac{T^2}{2} \). Solving gives \( b_1=\frac{T}{2} \), \( a_1=-\frac{T}{2} \), hence

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

Problem 3 (Lag-Chain Convergence): Prove that \( \left(1+\frac{sT}{N}\right)^{-N}\to e^{-sT} \) as \( N\to\infty \).

Solution: Use the exponential limit \( \lim_{N\to\infty}(1+z/N)^N=e^z \) with \( z=sT \), then invert both sides to obtain the result.

Problem 4 (Phase Lag and Crossover): Let \( G_p(s)=\frac{1}{s+1} \) and delay \( e^{-sT} \). The open-loop is \( L(j\omega)=K\frac{e^{-j\omega T}}{j\omega+1} \). Show that the phase is \( \angle L(j\omega) = -\arctan(\omega) - \omega T \). For a given \( T \), explain qualitatively why increasing \( T \) reduces the maximum stabilizing gain \( K \).

Solution: Since \( \angle\frac{1}{j\omega+1}=-\arctan(\omega) \) and \( \angle e^{-j\omega T}=-\omega T \), we get \( \angle L=-\arctan(\omega)-\omega T \). Increasing \( T \) adds extra phase lag growing linearly with \( \omega \), so the frequency where the phase approaches \( -\pi \) shifts to lower values; at those lower frequencies the magnitude \( |1/(j\omega+1)| \) is larger, so a smaller \( K \) is needed to keep \( |L(j\omega)| \) below 1 near the critical phase condition.

Problem 5 (Stability/Realizability Trade-off): Explain why the Padé \([1/1]\) approximation can create a right-half-plane zero, while the lag-chain approximation cannot.

Solution: Padé is constructed to match a local power series of the exponential and may require numerator/denominator polynomials that place zeros anywhere in the complex plane, including the right half-plane (non-minimum phase). The lag-chain is a product of stable first-order factors \( 1/(1+s\tau) \) and therefore has poles only in the left half-plane and zeros only at infinity (minimum phase).

We derived the transport delay element directly from the 1D convection PDE using characteristics, obtaining the exact shift \( y(t)=u(t-T) \) and the Laplace factor \( e^{-sT} \). Because this operator is non-rational, we built two finite-dimensional approximation families: Padé rational models (high local accuracy but potential non-minimum-phase artifacts) and stable pseudo-distributed lag chains (physically realizable, minimum-phase, converging to the delay as \( N \) grows). These approximations serve as the bridge from distributed transport physics to lumped simulation blocks used throughout system dynamics.