Chapter 10: Time-Domain Response of Linear Systems

Lesson 4: Higher-Order Systems: Dominant Poles and Model Reduction

Real plants frequently yield third- and higher-order transfer functions. This lesson develops a rigorous time-domain view of why a small subset of poles often governs the observable transient, and how this motivates principled low-order approximations. We derive dominant-pole theorems, quantitative dominance criteria, and reduction recipes (first- and second-order) that preserve essential time-response features while enabling analysis and design.

1. Motivation and Scope

In prior lessons, we obtained closed-form step responses for first- and second-order systems and connected pole locations to rise/settling time and overshoot. For higher-order stable LTI systems, the step response is still a superposition of exponentials (and possibly exponentially modulated sinusoids). The key practical observation is:

\( \text{Slow modes decay last and often dominate the visible transient.} \)

This lesson formalizes that idea using partial fractions (transfer-function view) and modal decomposition (state-space view), then constructs reduced models that approximate the original response in the time domain.

flowchart TD
  T0["Early time: many modes contribute"] --> T1["Middle time: fast modes decay"]
  T1 --> T2["Late time: slow (dominant) modes remain"]
  T2 --> APP["Approximate model: keep dominant modes only"]
        

2. Higher-Order Step Response as a Sum of Modes

Consider a proper, stable SISO transfer function \( G(s)=\dfrac{N(s)}{D(s)} \) with \( \deg N \le \deg D \) and all poles in the open left half-plane. For a unit-step input, \( U(s)=\dfrac{1}{s} \), hence \( Y(s)=\dfrac{G(s)}{s} \).

If all poles of \( \dfrac{G(s)}{s} \) are simple, we may write a partial-fraction expansion

\[ \frac{G(s)}{s} = \frac{G(0)}{s} + \sum_{i=1}^{n} \frac{r_i}{s-p_i}, \quad \Re(p_i) < 0. \]

Taking inverse Laplace yields the step response

\[ y(t)=G(0) + \sum_{i=1}^{n} r_i e^{p_i t}, \quad t \ge 0. \]

When poles come as a complex-conjugate pair \( p=\sigma \pm j\omega_d \) (with \( \sigma<0 \)), the corresponding time term is real and can be written as an exponentially decaying sinusoid:

\[ r e^{(\sigma + j\omega_d)t} + \bar{r}\, e^{(\sigma - j\omega_d)t} = 2|r| e^{\sigma t}\cos(\omega_d t + \phi). \]

Thus, higher-order responses are sums of (i) decaying exponentials and (ii) decaying sinusoids, each tied to a pole. We define an associated (real) time constant for a real pole \( p \) as \( \tau = -\dfrac{1}{\Re(p)} \) (for stable poles where \( \Re(p)<0 \)).

3. Residues for Simple Poles and a Practical Formula

For reduction and error estimation, we need coefficients \( r_i \). A key identity for simple poles is:

Proposition (Residue at a simple pole): Let

\[ F(s)=\frac{P(s)}{Q(s)}, \quad Q(p)=0,\quad Q'(p)\neq 0, \]

then the residue of \( F \) at \( s=p \) is

\[ \operatorname{Res}(F,p)=\frac{P(p)}{Q'(p)}. \]

Proof: Since \( p \) is a simple root, factor

\[ Q(s)=(s-p)\,Q_1(s), \quad Q_1(p)=Q'(p). \]

Then

\[ F(s)=\frac{P(s)}{(s-p)Q_1(s)} = \frac{1}{s-p}\cdot \frac{P(s)}{Q_1(s)}. \]

The residue is the coefficient of \( \dfrac{1}{s-p} \), i.e.

\[ \operatorname{Res}(F,p)=\lim_{s\,\to\,p}(s-p)F(s)=\lim_{s\,\to\,p}\frac{P(s)}{Q_1(s)}=\frac{P(p)}{Q_1(p)}=\frac{P(p)}{Q'(p)}. \]

This completes the proof. □

For the step response, apply this to \( F(s)=\dfrac{G(s)}{s}=\dfrac{N(s)}{sD(s)} \). For each system pole \( p_i \) (root of \( D \)), the corresponding step residue is \( r_i=\dfrac{N(p_i)}{p_i D'(p_i)} \) (simple poles).

4. Dominant Poles and an Asymptotic Dominance Theorem

Intuitively, among stable exponentials, the term with real part closest to zero decays slowest. This motivates a precise notion of dominance.

Definition (Dominant pole set): For stable poles \( \{p_i\}_{i=1}^n \), define \( \alpha^\star = \max_i \Re(p_i) \) (note \( \alpha^\star<0 \)). Any pole \( p_k \) satisfying \( \Re(p_k)=\alpha^\star \) is called dominant. When the dominant poles form a complex-conjugate pair, we speak of a dominant second-order mode.

Theorem (Asymptotic dominance of the slowest decay): Let

\[ f(t)=\sum_{i=1}^{n} c_i e^{p_i t}, \quad \Re(p_i)<0, \]

and assume there is a unique pole \( p_d \) such that \( \Re(p_d)=\alpha^\star \), with \( c_d\neq 0 \). Then

\[ \lim_{t\,\to\,\infty}\frac{f(t)}{c_d e^{p_d t}}=1. \]

Proof: Write

\[ \frac{f(t)}{c_d e^{p_d t}} = 1 + \sum_{i\neq d}\frac{c_i}{c_d} e^{(p_i-p_d)t}. \]

For each \( i\neq d \), we have \( \Re(p_i-p_d)=\Re(p_i)-\Re(p_d)<0 \). Hence \( e^{(p_i-p_d)t} \to 0 \) as \( t→\infty \). Therefore the sum tends to zero and the limit equals 1. □

Interpretation: at sufficiently large times, the transient is essentially the dominant term(s). This is the mathematical basis for dominant-pole approximations and simple model reduction.

5. Quantitative Dominance Criteria and Error Bounds

Practical reduction requires a quantitative rule: when is it safe to neglect a pole? Let the step response be \( y(t)=G(0)+\sum r_i e^{p_i t} \). Suppose we split poles into a retained (slow) set \( \mathcal{S} \) and discarded (fast) set \( \mathcal{F} \).

Define a reduced response by keeping only \( \mathcal{S} \):

\[ y_{\text{red}}(t)=G(0)+\sum_{i\in\mathcal{S}} r_i e^{p_i t}. \]

The error is

\[ e(t)=y(t)-y_{\text{red}}(t)=\sum_{i\in\mathcal{F}} r_i e^{p_i t}. \]

Simple bound (envelope): Using the triangle inequality, for \( t\ge 0 \):

\[ |e(t)| \le \sum_{i\in\mathcal{F}} |r_i|\, e^{\Re(p_i)t}. \]

If all discarded poles satisfy \( \Re(p_i)\le -\beta \) for some \( \beta>0 \), then

\[ |e(t)| \le e^{-\beta t}\sum_{i\in\mathcal{F}} |r_i|. \]

This yields a useful dominance guideline: if the fastest neglected pole real parts are much more negative than the retained pole real parts, then at time scales relevant to the slow poles, the neglected contribution is exponentially small.

Heuristic (time-scale separation rule): Let the dominant real part be \( \alpha^\star \). If all neglected poles satisfy \( \Re(p_i) \le m\,\alpha^\star \) with \( m \ge 5 \) (remember \( \alpha^\star<0 \)), then the neglected exponentials decay at least \( e^{(m-1)\alpha^\star t} \) relative to the dominant decay, often making their effect negligible after a short initial transient.

6. Constructing Low-Order Approximations

6.1 First-order approximation (single dominant real pole)

If the dominant pole is real \( p_d=\alpha^\star \), a natural approximation is a first-order model

\[ G_1(s)=\frac{K}{\tau s + 1}, \quad \tau=-\frac{1}{p_d}, \quad K=G(0). \]

This preserves the steady-state gain and the dominant time constant. If a more accurate low-frequency match is desired, match the first two Maclaurin coefficients of \( G(s) \) at \( s=0 \). Let \( g_0=G(0) \) and \( g_1=G'(0) \). Expanding \( \dfrac{K}{\tau s+1}=K(1-\tau s+\cdots) \) gives

\[ g_0=K,\quad g_1=-K\tau \;\;\Rightarrow\;\; \tau=-\frac{g_1}{g_0}. \]

6.2 Second-order approximation (dominant complex pair)

If the dominant poles are \( p_{1,2}=\sigma \pm j\omega_d \) with \( \sigma<0 \), define the standard second-order parameters (from Lesson 2–3):

\[ \omega_n=\sqrt{\sigma^2+\omega_d^2}, \quad \zeta=-\frac{\sigma}{\omega_n}. \]

A standard reduced model that preserves DC gain is

\[ G_2(s)=K\frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}, \quad K=G(0). \]

This is appropriate when the remaining poles are significantly faster and do not introduce comparable oscillatory content.

6.3 Modal truncation viewpoint in state space

Using the state-space form (Chapter 8), \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}u \), \( y=\mathbf{C}\mathbf{x}+\mathbf{D}u \). If \( \mathbf{A} \) is diagonalizable, \( \mathbf{A}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1} \) with \( \mathbf{\Lambda}=\operatorname{diag}(\lambda_i) \), then

\[ e^{\mathbf{A}t}=\mathbf{V}e^{\mathbf{\Lambda}t}\mathbf{V}^{-1} =\mathbf{V}\operatorname{diag}(e^{\lambda_i t})\mathbf{V}^{-1}. \]

Therefore, the state and output naturally decompose into modes \( e^{\lambda_i t} \). Modal truncation keeps modes with \( \Re(\lambda_i) \) closest to zero and discards fast-decaying modes.

7. Reduction Workflow

The following workflow is appropriate for the time-domain focus of this chapter: compute poles, identify dominant subset, build a low-order approximation preserving steady-state gain, and validate by comparing step responses and performance metrics (Lesson 5).

flowchart TD
  M0["Start from G(s)=N(s)/D(s) (stable)"] --> P0["Compute poles p_i (roots of D)"]
  P0 --> S0["Sort poles by real part: closest to 0 are slowest"]
  S0 --> C0["Choose retained set S: dominant pole or dominant complex pair"]
  C0 --> G0["Set reduced order: first-order or second-order"]
  G0 --> K0["Preserve DC gain: set K = G(0)"]
  K0 --> V0["Validate: compare step responses y(t) and y_red(t)"]
  V0 --> E0["If error too large: retain more poles or match more moments"]
        

8. Python Implementation (dominant pole / dominant pair reduction)

Libraries commonly used for system dynamics in Python include control (Control Systems Library) and scipy.signal. The example below creates a 4th-order stable system with a dominant complex pair and two fast real poles, then builds a 2nd-order approximation from the dominant pair and compares step responses.


import numpy as np
import control as ct

# Example higher-order plant:
# dominant pair at -1 ± j1.5, fast poles at -8 and -12, DC gain chosen by numerator scaling
poles = [-1 + 1.5j, -1 - 1.5j, -8, -12]
D = np.poly(poles)              # denominator coefficients
K_dc = 2.0                      # desired DC gain G(0)
# Choose numerator N(s) as constant to set DC gain: G(0) = N(0)/D(0)
N0 = K_dc * np.polyval(D, 0.0)
N = [N0]

G = ct.tf(N, D)

# Compute poles and select dominant pair (closest real part to 0)
p = ct.pole(G)
p_sorted = sorted(p, key=lambda z: np.real(z), reverse=True)  # largest real part first (closest to 0)
p_dom = p_sorted[:2]  # dominant pair here
sigma = np.real(p_dom[0])
omega_d = np.imag(p_dom[0])

omega_n = np.sqrt(sigma**2 + omega_d**2)
zeta = -sigma / omega_n

# Reduced 2nd-order model preserving DC gain
G2 = K_dc * ct.tf([omega_n**2], [1, 2*zeta*omega_n, omega_n**2])

# Step response comparison
t = np.linspace(0, 8, 2000)
t1, y1 = ct.step_response(G, T=t)
t2, y2 = ct.step_response(G2, T=t)

# Quick error measures (time-domain)
err_inf = np.max(np.abs(y1 - y2))
err_rel = err_inf / max(1e-12, np.max(np.abs(y1)))

print("Original poles:", p)
print("Dominant pair:", p_dom)
print("Reduced params: zeta=%.4f, omega_n=%.4f" % (zeta, omega_n))
print("Max abs error:", err_inf)
print("Max rel error:", err_rel)

# If you want performance metrics (rise/settling/overshoot), use ct.step_info
info_full = ct.step_info(G)
info_red  = ct.step_info(G2)
print("Full step info:", info_full)
print("Reduced step info:", info_red)
      

Notes: (1) Dominance is assessed here by real part ordering; in practice, residues also matter (Section 5). (2) If zeros strongly shape the transient, preserving them (or matching moments) may be necessary.

9. MATLAB and Simulink Implementation

MATLAB’s Control System Toolbox provides direct access to poles and step responses. The script below builds the same structure: a higher-order transfer function, then a dominant-pair second-order approximation preserving DC gain.


% Higher-order example: dominant pair -1 ± j1.5, fast poles -8 and -12
p = [-1+1.5j, -1-1.5j, -8, -12];
D = poly(p);              % denominator coefficients
Kdc = 2.0;                % desired DC gain
N0 = Kdc * polyval(D, 0); % choose constant numerator for DC gain
G = tf(N0, D);

% Poles and dominant pair selection (closest real part to 0)
pp = pole(G);
[~, idx] = sort(real(pp), 'descend');
pdom = pp(idx(1:2));
sigma = real(pdom(1));
omega_d = imag(pdom(1));

omega_n = sqrt(sigma^2 + omega_d^2);
zeta = -sigma/omega_n;

G2 = Kdc * tf([omega_n^2], [1, 2*zeta*omega_n, omega_n^2]);

% Step response comparison
t = linspace(0, 8, 2000);
[y_full, t_full] = step(G, t);
[y_red,  t_red ] = step(G2, t);

err_inf = max(abs(y_full - y_red));
err_rel = err_inf / max(1e-12, max(abs(y_full)));

disp('Original poles:'), disp(pp.')
disp('Dominant pair:'), disp(pdom.')
fprintf('Reduced params: zeta=%.4f, omega_n=%.4f\n', zeta, omega_n);
fprintf('Max abs error: %.6g\n', err_inf);
fprintf('Max rel error: %.6g\n', err_rel);

% Step metrics (rise/settling/overshoot) if needed:
info_full = stepinfo(G);
info_red  = stepinfo(G2);
disp(info_full), disp(info_red)
      

Simulink workflow (conceptual): Use two Transfer Fcn blocks (one for \( G(s) \), one for \( G_2(s) \)), feed both by a Step block, and compare in a Scope. Ensure identical simulation time and solver settings. This directly visualizes how fast poles affect only early-time transients.

10. C++ Implementation (residue-based step response + reduction)

In C++, common building blocks for system dynamics include Eigen (linear algebra) and standard complex arithmetic. Below is a compact demonstration that: (i) computes poles by eigenvalues of a companion matrix, (ii) computes step-response residues using \( r_i=\dfrac{N(p_i)}{p_i D'(p_i)} \) (simple poles), (iii) truncates to the two dominant poles (largest real parts), and (iv) evaluates \( y(t)=G(0)+\sum r_i e^{p_i t} \).


#include <iostream>
#include <vector>
#include <complex>
#include <algorithm>
#include <Eigen/Dense>

using cd = std::complex<double>;

static cd poly_eval(const std::vector<double>& a, cd s) {
  // a[0] s^n + a[1] s^(n-1) + ... + a[n]
  cd y = 0.0;
  for (size_t i = 0; i < a.size(); ++i) y = y * s + a[i];
  return y;
}

static std::vector<double> poly_derivative(const std::vector<double>& a) {
  // derivative of a[0] s^n + ... + a[n]
  int n = (int)a.size() - 1;
  std::vector<double> d;
  d.reserve((size_t)std::max(0, n));
  for (int i = 0; i < n; ++i) {
    d.push_back(a[i] * (n - i));
  }
  return d;
}

static std::vector<cd> companion_roots(const std::vector<double>& den) {
  // den = [1, a1, a2, ..., an] for s^n + a1 s^(n-1) + ... + an
  int n = (int)den.size() - 1;
  Eigen::MatrixXd C = Eigen::MatrixXd::Zero(n, n);
  for (int i = 0; i < n - 1; ++i) C(i + 1, i) = 1.0;
  for (int j = 0; j < n; ++j) C(0, j) = -den[j + 1];

  Eigen::EigenSolver<Eigen::MatrixXd> es(C);
  std::vector<cd> roots;
  roots.reserve((size_t)n);
  for (int i = 0; i < n; ++i) {
    auto ev = es.eigenvalues()(i);
    roots.emplace_back(ev.real(), ev.imag());
  }
  return roots;
}

int main() {
  // Example: poles = -1±j1.5, -8, -12
  // Denominator coefficients from expansion:
  // D(s) = (s^2 + 2s + 3.25)(s+8)(s+12) = s^4 + 22 s^3 + 167.25 s^2 + 487 s + 312
  std::vector<double> D = {1.0, 22.0, 167.25, 487.0, 312.0};

  // Numerator constant set for DC gain Kdc = 2: N0 = Kdc * D(0) = 2*312 = 624
  std::vector<double> N = {624.0};

  // Compute poles
  auto poles = companion_roots(D);

  // Sort by real part descending (closest to 0 first)
  std::sort(poles.begin(), poles.end(), [](cd a, cd b){
    return a.real() > b.real();
  });

  // Residues for step response: r_i = N(p_i)/(p_i D'(p_i))
  auto Dp = poly_derivative(D);

  std::vector<cd> r;
  r.reserve(poles.size());
  for (auto p : poles) {
    cd num = poly_eval(N, p);
    cd denp = poly_eval(Dp, p);
    r.push_back( num / (p * denp) );
  }

  // DC gain G(0) = N(0)/D(0)
  double G0 = N.back() / D.back();

  // Keep only 2 dominant poles (dominant complex pair)
  std::vector<cd> poles_keep = {poles[0], poles[1]};
  std::vector<cd> r_keep     = {r[0],     r[1]};

  // Evaluate step response
  double T = 8.0;
  int M = 2000;
  for (int k = 0; k <= M; ++k) {
    double t = T * k / M;
    cd y_full = G0;
    cd y_red  = G0;
    for (size_t i = 0; i < poles.size(); ++i) y_full += r[i] * std::exp(poles[i] * t);
    for (size_t i = 0; i < poles_keep.size(); ++i) y_red  += r_keep[i] * std::exp(poles_keep[i] * t);

    if (k % 200 == 0) {
      std::cout << "t=" << t << "  y_full=" << y_full.real()
                << "  y_red="  << y_red.real()
                << "  err="    << std::abs(y_full - y_red) << "\n";
    }
  }

  std::cout << "Poles (sorted):\n";
  for (auto p : poles) std::cout << "  " << p << "\n";
  std::cout << "G(0)=" << G0 << "\n";
  return 0;
}
      

This implementation follows the theory directly: pole computation, residue formula (Section 3), truncation, and time-domain reconstruction. For production-grade work, you would add robustness checks (repeated poles, numerical conditioning, and real-valued reconstruction).

11. Java Implementation (Apache Commons Math)

In Java, typical mathematical dependencies for system dynamics tasks include Apache Commons Math (polynomials, eigenvalues, complex arithmetic). The snippet below shows the same companion-matrix approach and residue-based step response evaluation. (For brevity, plotting is omitted; you can export samples to CSV.)


import java.util.*;
import org.apache.commons.math3.complex.Complex;
import org.apache.commons.math3.linear.*;

public class DominantPoleReduction {

  static Complex polyEval(double[] a, Complex s) {
    // a[0] s^n + ... + a[n]
    Complex y = Complex.ZERO;
    for (double coeff : a) y = y.multiply(s).add(coeff);
    return y;
  }

  static double[] polyDerivative(double[] a) {
    int n = a.length - 1;
    double[] d = new double[n];
    for (int i = 0; i < n; i++) d[i] = a[i] * (n - i);
    return d;
  }

  static List<Complex> companionRoots(double[] den) {
    // den = [1, a1, ..., an] for s^n + a1 s^(n-1) + ... + an
    int n = den.length - 1;
    RealMatrix C = MatrixUtils.createRealMatrix(n, n);
    for (int i = 0; i < n - 1; i++) C.setEntry(i + 1, i, 1.0);
    for (int j = 0; j < n; j++) C.setEntry(0, j, -den[j + 1]);

    EigenDecomposition ed = new EigenDecomposition(C);
    List<Complex> roots = new ArrayList<>();
    for (int i = 0; i < n; i++) {
      roots.add(new Complex(ed.getRealEigenvalue(i), ed.getImagEigenvalue(i)));
    }
    return roots;
  }

  public static void main(String[] args) {
    // D(s)=s^4 + 22 s^3 + 167.25 s^2 + 487 s + 312
    double[] D = new double[]{1.0, 22.0, 167.25, 487.0, 312.0};
    double[] N = new double[]{624.0}; // DC gain 2: N0 = 2*312

    List<Complex> poles = companionRoots(D);
    poles.sort((a,b) -> Double.compare(b.getReal(), a.getReal())); // descending real part

    double[] Dp = polyDerivative(D);

    List<Complex> residues = new ArrayList<>();
    for (Complex p : poles) {
      Complex num = polyEval(N, p);
      Complex denp = polyEval(Dp, p);
      residues.add(num.divide(p.multiply(denp))); // r_i = N(p_i)/(p_i D'(p_i))
    }

    double G0 = N[N.length - 1] / D[D.length - 1];

    // keep 2 dominant poles
    List<Complex> polesKeep = poles.subList(0, 2);
    List<Complex> resKeep   = residues.subList(0, 2);

    double T = 8.0;
    int M = 2000;
    for (int k = 0; k <= M; k += 200) {
      double t = T * k / M;
      Complex yFull = new Complex(G0, 0.0);
      Complex yRed  = new Complex(G0, 0.0);

      for (int i = 0; i < poles.size(); i++) {
        yFull = yFull.add(residues.get(i).multiply(poles.get(i).multiply(t).exp()));
      }
      for (int i = 0; i < polesKeep.size(); i++) {
        yRed = yRed.add(resKeep.get(i).multiply(polesKeep.get(i).multiply(t).exp()));
      }

      System.out.println("t=" + t + "  y_full=" + yFull.getReal() +
                         "  y_red=" + yRed.getReal() +
                         "  err=" + yFull.subtract(yRed).abs());
    }

    System.out.println("Poles:");
    for (Complex p : poles) System.out.println("  " + p);
    System.out.println("G(0)=" + G0);
  }
}
      

If you need visualization, export \( (t,y) \) samples to CSV and plot in a separate tool. The core reduction idea remains residue-based truncation guided by pole real parts and residue magnitudes.

12. Wolfram Mathematica Implementation

Mathematica has native transfer-function manipulation via TransferFunctionModel, pole computation, and symbolic/analytic response functions. The code below computes poles, extracts the dominant pair, forms the second-order approximation preserving DC gain, and compares step responses.


(* Define higher-order transfer function *)
poles = {-1 + I*1.5, -1 - I*1.5, -8, -12};
D[s_] := Expand[Times @@ (s - # & /@ poles)];
Kdc = 2.0;
N0 = Kdc * (D[0] /. s -> 0);  (* D(0) *)
G[s_] := N0 / D[s];

(* Poles of G(s) are poles list; confirm numerically *)
p = N[NSolve[D[s] == 0, s][[All, 1, 2]]];

(* Select dominant pair by sorting real parts descending *)
pSorted = SortBy[p, -Re[#] &];
pDom = Take[pSorted, 2];
sigma = Re[pDom[[1]]];
omegaD = Im[pDom[[1]]];

omegaN = Sqrt[sigma^2 + omegaD^2];
zeta = -sigma/omegaN;

G2[s_] := Kdc * omegaN^2/(s^2 + 2*zeta*omegaN*s + omegaN^2);

(* Step responses: y(t) = InverseLaplaceTransform(G(s)/s, s, t) *)
yFull[t_] := Evaluate[InverseLaplaceTransform[G[s]/s, s, t]];
yRed[t_]  := Evaluate[InverseLaplaceTransform[G2[s]/s, s, t]];

(* Compare numerically on a grid *)
tgrid = N@Subdivide[0, 8, 2000];
errInf = Max[Abs[(yFull /@ tgrid) - (yRed /@ tgrid)]];
{zeta, omegaN, errInf}
      

Mathematica can keep responses symbolic for additional analysis (e.g., asymptotic behavior as \( t→\infty \)).

13. Problems and Solutions

Problem 1 (Dominant complex pair & equivalent \(\zeta,\omega_n\)): A stable system has dominant poles at \( p_{1,2}=-1 \pm j2 \). Compute \( \omega_n \) and \( \zeta \) for the corresponding second-order approximation.

Solution: Here \( \sigma=-1 \), \( \omega_d=2 \), hence

\[ \omega_n=\sqrt{\sigma^2+\omega_d^2}=\sqrt{1+4}=\sqrt{5},\quad \zeta=-\frac{\sigma}{\omega_n}=\frac{1}{\sqrt{5}}. \]

Problem 2 (Asymptotic dominance ratio): Let \( f(t)=3e^{-t}+2e^{-5t} \). Show that \( \lim_{t\,\to\,\infty}\dfrac{f(t)}{3e^{-t}}=1 \) and interpret the result.

Solution:

\[ \frac{f(t)}{3e^{-t}} = 1 + \frac{2}{3}e^{-4t}. \]

Since \( e^{-4t}\to 0 \) as \( t→\infty \), the ratio tends to 1. The slow mode \( e^{-t} \) dominates late-time behavior; the fast mode \( e^{-5t} \) matters only at early times.

Problem 3 (First-order reduction by low-frequency moment matching): Consider \( G(s)=\dfrac{s+3}{(s+1)(s+5)(s+20)} \). Find a first-order approximation \( G_1(s)=\dfrac{K}{\tau s +1} \) that matches \( G(0) \) and \( G'(0) \).

Solution: First compute

\[ G(0)=\frac{3}{(1)(5)(20)}=\frac{3}{100}=0.03 \;\Rightarrow\; K=0.03. \]

Write \( G(s)=\dfrac{N(s)}{D(s)} \) with \( N(s)=s+3 \), \( D(s)=(s+1)(s+5)(s+20) \). Then

\[ G'(s)=\frac{N'(s)D(s)-N(s)D'(s)}{D(s)^2}. \]

Evaluate at \( s=0 \). We have \( N'(0)=1 \), \( N(0)=3 \), \( D(0)=100 \). Next compute \( D'(0) \). Expand:

\[ D(s)=(s+1)(s+5)(s+20)=(s^2+6s+5)(s+20)=s^3+26s^2+125s+100, \]

hence \( D'(s)=3s^2+52s+125 \), so \( D'(0)=125 \). Therefore

\[ G'(0)=\frac{(1)(100)-(3)(125)}{100^2}=\frac{100-375}{10000}=-\frac{275}{10000}=-0.0275. \]

For \( G_1(s)=\dfrac{K}{\tau s+1}=K(1-\tau s+\cdots) \), we have \( G_1'(0)=-K\tau \). Matching gives

\[ -K\tau = -0.0275 \;\Rightarrow\; \tau=\frac{0.0275}{0.03}=0.916\overline{6}. \]

Thus \( G_1(s)=\dfrac{0.03}{0.9167\,s+1} \) (approximately).

Problem 4 (Residue formula for step response): Let \( G(s)=\dfrac{N(s)}{D(s)} \) have a simple pole at \( s=p \). Show that the residue of \( \dfrac{G(s)}{s} \) at \( s=p \) equals \( \dfrac{N(p)}{pD'(p)} \).

Solution: Apply the result of Section 3 to

\[ F(s)=\frac{G(s)}{s}=\frac{N(s)}{sD(s)}. \]

Here the pole at \( s=p \) comes from \( D(p)=0 \). Define \( Q(s)=sD(s) \), \( P(s)=N(s) \). Then \( Q'(s)=D(s)+sD'(s) \), and at \( s=p \) we have \( D(p)=0 \), hence \( Q'(p)=pD'(p) \). Therefore

\[ \operatorname{Res}\!\left(\frac{G(s)}{s},p\right)=\frac{P(p)}{Q'(p)}=\frac{N(p)}{pD'(p)}. \]

Problem 5 (Error bound from neglected poles): Suppose the discarded poles satisfy \( \Re(p_i)\le -10 \) and \( \sum_{i\in\mathcal{F}}|r_i| \le 0.5 \). Bound \( |e(t)| \) for \( t\ge 0.3 \).

Solution: From Section 5,

\[ |e(t)| \le e^{-10t}\sum_{i\in\mathcal{F}}|r_i| \le 0.5\,e^{-10t}. \]

For \( t\ge 0.3 \), \( |e(t)| \le 0.5 e^{-3} \approx 0.0249 \). This quantifies how quickly fast poles become negligible.

14. Summary

Higher-order stable LTI step responses decompose into sums of modal terms \( r_i e^{p_i t} \). Poles with real parts closest to zero decay slowest and dominate late-time transients. This lesson proved an asymptotic dominance theorem, provided residue formulas for constructing time responses, and developed practical reduction workflows that keep a dominant real pole (first-order) or dominant complex pair (second-order) while preserving DC gain. The resulting reduced models enable time-domain analysis using the first/second-order tools developed earlier in Chapter 10.

15. References

  1. Moore, B.C. (1981). Principal Component Analysis in Linear Systems: Controllability, Observability, and Model Reduction. IEEE Transactions on Automatic Control, 26(1), 17–32. :contentReference[oaicite:0]{index=0}
  2. Pernebo, L., & Silverman, L.M. (1982). Model Reduction by Balanced State Space Representations. IEEE Transactions on Automatic Control, 27(2), 382–387. :contentReference[oaicite:1]{index=1}
  3. Enns, D.F. (1984). Model Reduction with Balanced Realizations: An Error Bound and a Frequency Weighted Generalization. Proceedings of the 23rd IEEE Conference on Decision and Control. :contentReference[oaicite:2]{index=2}
  4. Glover, K. (1984). All Optimal Hankel-Norm Approximations of Linear Multivariable Systems and Their L∞-Error Bounds. International Journal of Control, 39(6), 1115–1193. :contentReference[oaicite:3]{index=3}
  5. Antoulas, A.C., Sorensen, D.C., & Gugercin, S. (2001). A Survey of Model Reduction Methods for Large-Scale Systems. Contemporary Mathematics, 280, 193–219. :contentReference[oaicite:4]{index=4}
  6. Rommes, J. (2008). Computing Transfer Function Dominant Poles of Large-Scale Second-Order Dynamical Systems. SIAM Journal on Scientific Computing. :contentReference[oaicite:5]{index=5}