Chapter 24: Robustness Analysis with Classical Tools

Lesson 1: Robust Stability via Gain and Phase Margins

This lesson develops a quantitative understanding of how gain margin (GM) and phase margin (PM) measure distance to instability for linear feedback systems. We interpret GM/PM both geometrically via the Nyquist plot and analytically via inequalities that bound gain, phase, and time-delay uncertainties. We then implement margin computations in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, with particular emphasis on their use in robotic motion and joint control.

1. Robust Stability – Conceptual Overview

Consider a unity-feedback SISO loop with controller \( C(s) \), plant \( G(s) \), and loop transfer function \( L(s) = C(s)G(s) \). The closed-loop complementary sensitivity and sensitivity functions (introduced in earlier chapters) are

\[ T(s) = \frac{L(s)}{1+L(s)}, \quad S(s) = \frac{1}{1+L(s)}, \quad L(s) = C(s)G(s). \]

A nominal design is internally stable if all poles of \( T(s) \) (equivalently of \(1+L(s)\)) lie in the left half-plane. Robust stability goes further: the closed loop should remain stable for all plants in an uncertainty set \( \mathcal{P} \) (e.g., due to parameter variations, unmodeled dynamics, and time delays).

Classical gain and phase margins provide a compact way to certify robust stability against:

  • Multiplicative or constant gain changes (actuator gain drift, load variations).
  • Phase errors (unmodeled time delay, sensor/actuator lags, approximate models).

Geometrically, GM and PM measure how close the Nyquist plot of \( L(j\omega) \) comes to the critical point \( -1 + j0 \); analytically, they induce inequalities guaranteeing stability under bounded gain and phase perturbations.

flowchart TD
  P["Nominal plant model G(s)"] --> Cblk["Controller C(s) (e.g. PID)"]
  Cblk --> Lblk["Open-loop L(s) = C(s) G(s)"]
  Lblk --> B["Bode / Nyquist analysis"]
  B --> M["Compute gain margin GM and phase margin PM"]
  M --> R["Compare with uncertainty (gain, phase, delay)"]
  R --> OK["If margins adequate: \nimplement on robot/plant"]
  R --> REDESIGN["If margins small: \nretune or redesign controller"]
        

2. Formal Definitions of Gain and Phase Margins

Let \( L(j\omega) \) be the loop transfer function evaluated on the imaginary axis. Assume that the closed loop with the nominal plant is stable and that the Nyquist plot of \( L(j\omega) \) does not pass exactly through \( -1 \).

The gain crossover frequency \( \omega_{\mathrm{gc}} \) and the phase crossover frequency \( \omega_{\mathrm{pc}} \) are defined by

\[ \left| L(j\omega_{\mathrm{gc}})\right| = 1, \qquad \arg L(j\omega_{\mathrm{pc}}) = -\pi. \]

The phase margin is

\[ \varphi_m = \pi + \arg L(j\omega_{\mathrm{gc}}), \]

i.e., how many radians (or degrees) of additional phase lag at unit magnitude can be tolerated before the phase reaches \( -\pi \).

The (upper) gain margin is

\[ G_m = \frac{1}{\left| L(j\omega_{\mathrm{pc}})\right|}, \qquad G_{m,\mathrm{dB}} = 20\log_{10} G_m. \]

Informally, the loop remains stable if the open-loop gain is multiplied by any constant factor less than \( G_m \), i.e., if the effective gain is increased up to \( G_m \) times its nominal value.

It is often convenient to introduce the minimal Nyquist distance

\[ d_{\min} := \min_{\omega \in \mathbb{R}} \left| 1 + L(j\omega) \right|, \]

which is the shortest distance between the Nyquist curve and the critical point \( -1 \). This distance directly quantifies robustness to additive perturbations of the loop transfer function.

3. Additive and Multiplicative Loop Uncertainty

3.1 Additive loop uncertainty

Suppose the true loop transfer function is

\[ \tilde{L}(s) = L(s) + \Delta L(s), \qquad \left| \Delta L(j\omega) \right| \leq \bar{\delta}_a(\omega). \]

The perturbed closed-loop characteristic equation is \( 1 + \tilde{L}(s) = 0 \). On the imaginary axis

\[ 1 + \tilde{L}(j\omega) = 1 + L(j\omega) + \Delta L(j\omega). \]

A sufficient Nyquist-based condition for robust stability is

\[ \left| \Delta L(j\omega) \right| < \left| 1 + L(j\omega)\right| \quad \forall \omega \in \mathbb{R}. \]

Proof sketch. For any frequency \( \omega \),

\[ \left| 1 + L(j\omega) + \Delta L(j\omega) \right| \geq \left| \left|1+L(j\omega)\right| - \left|\Delta L(j\omega)\right| \right| > 0 \quad \forall \omega. \]

Hence the perturbed Nyquist curve does not pass through \( -1 \), and (under the standard assumptions on the number of open-loop right-half-plane poles) the closed loop remains stable.

3.2 Multiplicative loop uncertainty

In many control and robotics models, a multiplicative description is more natural:

\[ \tilde{L}(s) = L(s)\bigl(1+\Delta_m(s)\bigr), \qquad \left|\Delta_m(j\omega)\right| \leq \bar{\delta}_m(\omega). \]

Then \( \Delta L(s) = L(s)\Delta_m(s) \), and the sufficient robust stability condition above becomes

\[ \bar{\delta}_m(\omega) < \frac{\left|1+L(j\omega)\right|} {\left|L(j\omega)\right|} \quad \forall \omega. \]

The right-hand side is large when the Nyquist curve is far from the critical point and small near the gain and phase crossover frequencies. In practice, the worst case typically occurs near \( \omega_{\mathrm{gc}} \).

3.3 Relation between phase margin and multiplicative uncertainty

At the gain crossover frequency, we have \( \left|L(j\omega_{\mathrm{gc}})\right| = 1 \) and \( \arg L(j\omega_{\mathrm{gc}}) = \varphi_{\mathrm{gc}} \). By definition,

\[ \varphi_m = \pi + \varphi_{\mathrm{gc}}. \]

Thus we can write

\[ L(j\omega_{\mathrm{gc}}) = \mathrm{e}^{j(-\pi + \varphi_m)} = -\mathrm{e}^{j\varphi_m}, \quad 1 + L(j\omega_{\mathrm{gc}}) = 1 - \mathrm{e}^{j\varphi_m}. \]

Then

\[ \left|1-\mathrm{e}^{j\varphi_m}\right|^2 = (1-\mathrm{e}^{j\varphi_m})(1-\mathrm{e}^{-j\varphi_m}) = 2 - 2\cos\varphi_m = 4\sin^2\frac{\varphi_m}{2}, \]

\[ \left|1+L(j\omega_{\mathrm{gc}})\right| = 2\left|\sin\frac{\varphi_m}{2}\right|. \]

Since \( \left|L(j\omega_{\mathrm{gc}})\right| = 1 \), the multiplicative robust stability condition at \( \omega_{\mathrm{gc}} \) becomes

\[ \left|\Delta_m(j\omega_{\mathrm{gc}})\right| < 2\left|\sin\frac{\varphi_m}{2}\right| =: \delta_{\max}. \]

Hence for small perturbations concentrated near \( \omega_{\mathrm{gc}} \), the admissible magnitude of multiplicative uncertainty is approximately

\[ \delta_{\max} \approx 2\sin\frac{\varphi_m}{2}. \]

For example, if \( \varphi_m = 60^\circ \), then \( \delta_{\max} \approx 2\sin 30^\circ = 1 \): the loop can tolerate roughly \( 100\% \) multiplicative model error at the crossover frequency without losing stability.

4. Phase Margin, Delay Margin, and Time-Delay Uncertainty

Unmodeled time delay is a major source of phase uncertainty in robotics and mechatronics (e.g., computation, communication, or sensor delays). Suppose an additional delay \( \tau_d \) appears in the loop:

\[ L_d(s) = L(s)\mathrm{e}^{-s \tau_d}. \]

The delay does not change the magnitude, but contributes an extra phase lag of \( -\omega \tau_d \) at frequency \( \omega \). In particular, at the nominal gain crossover frequency \( \omega_{\mathrm{gc}} \),

\[ \arg L_d(j\omega_{\mathrm{gc}}) = \arg L(j\omega_{\mathrm{gc}}) - \omega_{\mathrm{gc}} \tau_d. \]

The Nyquist curve remains at unit magnitude at \( \omega_{\mathrm{gc}} \), but its phase is rotated by \( -\omega_{\mathrm{gc}} \tau_d \). A sufficient condition for stability is that, at \( \omega_{\mathrm{gc}} \),

\[ \arg L_d(j\omega_{\mathrm{gc}}) > -\pi, \]

i.e.

\[ \arg L(j\omega_{\mathrm{gc}}) - \omega_{\mathrm{gc}} \tau_d > -\pi. \]

Adding \( \pi \) to both sides,

\[ \pi + \arg L(j\omega_{\mathrm{gc}}) > \omega_{\mathrm{gc}} \tau_d. \]

The left-hand side is precisely the phase margin \( \varphi_m \), therefore

\[ \varphi_m > \omega_{\mathrm{gc}} \tau_d \quad \Rightarrow \quad \tau_d < \frac{\varphi_m}{\omega_{\mathrm{gc}}} =: \tau_{\max}. \]

This is the classical delay margin approximation: a loop with phase margin \( \varphi_m \) and gain crossover frequency \( \omega_{\mathrm{gc}} \) can tolerate an additional delay of roughly \( \tau_{\max} \) before losing stability. In robotics, \( \tau_{\max} \) provides an upper bound on admissible sampling or communication delays.

5. Gain Margin and Constant Gain Uncertainty

Let the controller contain a tunable scalar gain \( k \), and write the loop transfer function as

\[ L(s;k) = k\,L_0(s), \]

where \( L_0(s) \) captures the remaining frequency dependence. Assume that the Nyquist plot of \( L(s;k_{\mathrm{nom}}) \) is known and that the corresponding closed loop is stable.

If the (upper) gain margin is \( G_m \), then any gain \( k \) in the interval

\[ 0 < k < G_m \, k_{\mathrm{nom}} \]

yields a stable closed loop (in the simplest case where there is no lower gain limit and only one relevant phase crossover). More generally, Bode/Nyquist tools may report upper and lower gain margins \( G_m^+ \) and \( G_m^- \), and the loop remains stable for

\[ \frac{1}{G_m^-} < \frac{k}{k_{\mathrm{nom}}} < G_m^+. \]

For example, if the nominal gain margin is \( G_m^+ = 6 \) (about \( 15.6 \) dB), then, ignoring the lower margin, the true gain in operation may be up to \( 6 \) times larger than the design value without destabilizing the loop. If model uncertainty suggests that the actuator gain may vary by at most a factor of \( 3 \), the condition \( G_m^+ \geq 3 \) is sufficient for robust stability against that particular source of uncertainty.

6. Algorithmic Use of Margins in Design

Combining the previous sections, a standard robust design workflow in classical control (for example, a robot joint position controller) is:

  1. Identify a nominal plant model \( G(s) \) from physics or system identification.
  2. Estimate ranges of parameter and delay uncertainty (e.g., load inertia, friction, computation and communication delays).
  3. Choose a controller structure \( C(s) \) (PID, lead, lead&lag, etc.).
  4. Shape the loop \( L(s) = C(s)G(s) \) to meet performance goals (bandwidth, overshoot, disturbance rejection).
  5. Compute \( G_m, \varphi_m, \omega_{\mathrm{gc}} \) and check:
    • delay margin \( \tau_{\max} = \varphi_m/\omega_{\mathrm{gc}} \)
    • approximate multiplicative uncertainty bound \( \delta_{\max} \approx 2\sin(\varphi_m/2) \)
  6. If the margins are insufficient compared to uncertainty estimates, increase phase margin (e.g., adding a lead compensator) or reduce bandwidth.
  7. Validate the design in simulation and, eventually, in robot experiments.

7. Python Implementation – Margins for a Robot Joint

As a simple robotics-motivated example, consider a rigid robot joint modeled as

\[ G(s) = \frac{K_t}{Js^2 + bs}, \]

where \( J \) is the joint inertia, \( b \) the viscous friction, and \( K_t \) a torque constant. We apply a PI controller

\[ C(s) = K_p + \frac{K_i}{s} = \frac{K_p s + K_i}{s}. \]

Using python-control (and optionally interfacing later with roboticstoolbox and ROS for full robot models), we can compute GM, PM, and delay margin:


import numpy as np
import control as ctl
import matplotlib.pyplot as plt

# Robot joint parameters (nominal)
J = 0.01     # kg m^2
b = 0.1      # N m s/rad
K_t = 0.5    # N m/A

# PI controller gains (to be tuned)
Kp = 20.0
Ki = 10.0

# Plant and controller transfer functions
G = ctl.tf([K_t], [J, b, 0.0])          # G(s) = K_t / (J s^2 + b s)
C = ctl.tf([Kp, Ki], [1.0, 0.0])        # C(s) = (Kp s + Ki) / s

L = C * G                               # open-loop
T = ctl.feedback(L, 1)                  # closed-loop

# Compute classical stability margins
gm, pm, w_gc, w_pc = ctl.margin(L)
gm_db = 20.0 * np.log10(gm) if gm is not None and gm > 0 else np.inf

print(f"Gain margin: {gm:.3f} ({gm_db:.2f} dB) at w_pc = {w_pc:.3f} rad/s")
print(f"Phase margin: {pm:.3f} deg at w_gc = {w_gc:.3f} rad/s")

# Approximate delay margin (radians for pm, then convert)
pm_rad = np.deg2rad(pm)
tau_max = pm_rad / w_gc
print(f"Approximate additional delay margin tau_max = {tau_max:.4f} s")

# Approximate admissible multiplicative uncertainty near crossover
delta_mult = 2.0 * np.sin(0.5 * pm_rad)
print(f"Approximate multiplicative uncertainty at crossover: delta_max = {delta_mult:.3f}")

# Bode plot with margin annotations
mag, phase, omega = ctl.bode(L, dB=True, Hz=False, omega_limits=(1e-1, 1e3), Plot=False)

plt.figure()
plt.subplot(2, 1, 1)
plt.semilogx(omega, 20.0 * np.log10(mag))
plt.ylabel("Magnitude (dB)")
plt.grid(True, which="both")

plt.subplot(2, 1, 2)
plt.semilogx(omega, phase * 180.0 / np.pi)
plt.xlabel("Frequency (rad/s)")
plt.ylabel("Phase (deg)")
plt.grid(True, which="both")

plt.tight_layout()
plt.show()

# NOTE (robotics context):
# In a robotics workflow, G and C are often derived using a library like
# 'roboticstoolbox' to obtain a linearized joint model:
#
#   from roboticstoolbox import DHRobot
#   # define robot, linearize around a configuration, then extract G(s)
#
# The same margin analysis on L(s) is then used to certify robust stability of
# the joint controller before deployment in ROS/ros2 control loops.
      

8. C++ Implementation – Evaluating Margins Numerically

In embedded robot controllers (e.g., running under ROS 2), the control loop is often implemented in C++. The following example shows how to evaluate the loop transfer function on a frequency grid and approximate gain and phase margins. Libraries such as Eigen (for linear algebra) and ros_control/ros2_control are commonly used in practice.


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

using std::complex;
using std::vector;

complex<double> eval_tf(const vector<double>& num,
                         const vector<double>& den,
                         double w)
{
    complex<double> s(0.0, w);
    complex<double> N(0.0, 0.0), D(0.0, 0.0);

    // Horner evaluation
    for (std::size_t i = 0; i < num.size(); ++i)
    {
        N = N * s + num[i];
    }
    for (std::size_t i = 0; i < den.size(); ++i)
    {
        D = D * s + den[i];
    }
    return N / D;
}

int main()
{
    // Robot joint plant G(s) = K_t / (J s^2 + b s)
    double J  = 0.01;
    double b  = 0.1;
    double Kt = 0.5;

    // Numerator and denominator for G(s)
    vector<double> G_num{Kt};
    vector<double> G_den{J, b, 0.0};

    // PI controller C(s) = (Kp s + Ki) / s
    double Kp = 20.0;
    double Ki = 10.0;
    vector<double> C_num{Kp, Ki};
    vector<double> C_den{1.0, 0.0};

    // Frequency grid (rad/s)
    double w_min = 0.1;
    double w_max = 100.0;
    int N = 2000;
    double log_w_min = std::log10(w_min);
    double log_w_max = std::log10(w_max);

    double best_gc_diff = 1e9;
    double best_gc = 0.0;
    double best_pc_diff = 1e9;
    double best_pc = 0.0;
    double phase_at_gc = 0.0;
    double mag_at_pc = 0.0;

    for (int k = 0; k < N; ++k)
    {
        double logw = log_w_min + (log_w_max - log_w_min) * k / (N - 1);
        double w = std::pow(10.0, logw);
        complex<double> s(0.0, w);

        complex<double> G = eval_tf(G_num, G_den, w);
        complex<double> C = eval_tf(C_num, C_den, w);
        complex<double> L = C * G;

        double mag = std::abs(L);
        double phase = std::arg(L); // radians

        // Track approximate gain crossover |L| = 1
        double gc_diff = std::fabs(mag - 1.0);
        if (gc_diff < best_gc_diff)
        {
            best_gc_diff = gc_diff;
            best_gc = w;
            phase_at_gc = phase;
        }

        // Track approximate phase crossover angle = -pi
        double pc_diff = std::fabs(phase + M_PI);
        if (pc_diff < best_pc_diff)
        {
            best_pc_diff = pc_diff;
            best_pc = w;
            mag_at_pc = mag;
        }
    }

    // Phase margin (deg)
    double pm = (M_PI + phase_at_gc) * 180.0 / M_PI;

    // Gain margin factor
    double gm = (mag_at_pc > 0.0) ? (1.0 / mag_at_pc) : INFINITY;
    double gm_db = 20.0 * std::log10(gm);

    // Delay margin approximation
    double pm_rad = pm * M_PI / 180.0;
    double tau_max = pm_rad / best_gc;

    std::cout << "Approx gain crossover omega_gc = " << best_gc << " rad/s\n";
    std::cout << "Approx phase crossover omega_pc = " << best_pc << " rad/s\n";
    std::cout << "Phase margin pm = " << pm << " deg\n";
    std::cout << "Gain margin gm = " << gm << " (" << gm_db << " dB)\n";
    std::cout << "Approx delay margin tau_max = " << tau_max << " s\n";

    // In a ROS2-based robot, the resulting tau_max can be compared against
    // control loop sampling, communication, and computation delays.

    return 0;
}
      

9. Java Implementation – Margins with Apache Commons Math

Java is frequently used in educational robotics (e.g., FRC robots with WPILib) and some industrial control platforms. Below is a sketch using org.apache.commons.math3.complex.Complex to compute margins for the same loop transfer function.


import org.apache.commons.math3.complex.Complex;
import java.util.function.DoubleFunction;

public class MarginExample {

    static Complex evalTF(double[] num, double[] den, double w) {
        Complex s = new Complex(0.0, w);
        Complex N = Complex.ZERO;
        Complex D = Complex.ZERO;

        for (double v : num) {
            N = N.multiply(s).add(v);
        }
        for (double v : den) {
            D = D.multiply(s).add(v);
        }
        return N.divide(D);
    }

    public static void main(String[] args) {

        // Plant G(s) = Kt / (J s^2 + b s)
        double J = 0.01;
        double b = 0.1;
        double Kt = 0.5;
        double[] Gnum = {Kt};
        double[] Gden = {J, b, 0.0};

        // PI controller C(s) = (Kp s + Ki) / s
        double Kp = 20.0;
        double Ki = 10.0;
        double[] Cnum = {Kp, Ki};
        double[] Cden = {1.0, 0.0};

        DoubleFunction<Complex> L = (double w) -> {
            Complex G = evalTF(Gnum, Gden, w);
            Complex C = evalTF(Cnum, Cden, w);
            return G.multiply(C);
        };

        double wMin = 0.1;
        double wMax = 100.0;
        int N = 2000;

        double bestGcDiff = 1e9;
        double bestGc = 0.0;
        double bestPcDiff = 1e9;
        double bestPc = 0.0;
        double phaseAtGc = 0.0;
        double magAtPc = 0.0;

        for (int k = 0; k < N; ++k) {
            double logw = Math.log10(wMin) + (Math.log10(wMax) - Math.log10(wMin)) * k / (N - 1);
            double w = Math.pow(10.0, logw);

            Complex Lw = L.apply(w);
            double mag = Lw.abs();
            double phase = Lw.getArgument(); // radians

            double gcDiff = Math.abs(mag - 1.0);
            if (gcDiff < bestGcDiff) {
                bestGcDiff = gcDiff;
                bestGc = w;
                phaseAtGc = phase;
            }

            double pcDiff = Math.abs(phase + Math.PI);
            if (pcDiff < bestPcDiff) {
                bestPcDiff = pcDiff;
                bestPc = w;
                magAtPc = mag;
            }
        }

        double pm = (Math.PI + phaseAtGc) * 180.0 / Math.PI;
        double gm = (magAtPc > 0.0) ? (1.0 / magAtPc) : Double.POSITIVE_INFINITY;
        double gmDb = 20.0 * Math.log10(gm);

        double pmRad = pm * Math.PI / 180.0;
        double tauMax = pmRad / bestGc;
        double deltaMult = 2.0 * Math.sin(0.5 * pmRad);

        System.out.println("omega_gc ~ " + bestGc + " rad/s");
        System.out.println("omega_pc ~ " + bestPc + " rad/s");
        System.out.println("Phase margin ~ " + pm + " deg");
        System.out.println("Gain margin ~ " + gm + " (" + gmDb + " dB)");
        System.out.println("Delay margin tau_max ~ " + tauMax + " s");
        System.out.println("delta_max (multiplicative) ~ " + deltaMult);

        // Integration with a robot framework like WPILib:
        // the gains Kp, Ki can be tuned so that the analytically computed margins
        // remain acceptable for the range of robot operating conditions.
    }
}
      

10. MATLAB/Simulink Implementation – Margins and Delay Margin

MATLAB and Simulink are standard tools for classical control and robotics. Assuming the Control System Toolbox (and optionally the Robotics System Toolbox) are available, one can compute margins and visualize Nyquist/Bode plots concisely:


J  = 0.01;
b  = 0.1;
Kt = 0.5;

Kp = 20;
Ki = 10;

% Plant and controller
G = tf(Kt, [J b 0]);          % G(s) = Kt / (J s^2 + b s)
C = tf([Kp Ki], [1 0]);       % C(s) = (Kp s + Ki) / s
L = series(C, G);
T = feedback(L, 1);

% Classical stability margins
[gm, pm, w_gc, w_pc] = margin(L);
gm_db = 20*log10(gm);

fprintf("Gain margin: %.3f (%.2f dB) at w_pc = %.3f rad/s\n", gm, gm_db, w_pc);
fprintf("Phase margin: %.3f deg at w_gc = %.3f rad/s\n", pm, w_gc);

% Delay and multiplicative uncertainty margins
pm_rad   = deg2rad(pm);
tau_max  = pm_rad / w_gc;
deltaMax = 2 * sin(0.5 * pm_rad);
fprintf("Approx delay margin tau_max = %.4f s\n", tau_max);
fprintf("Approx multiplicative uncertainty delta_max = %.3f\n", deltaMax);

% Bode and Nyquist plots
figure;
margin(L); grid on;           % Bode with margins highlighted

figure;
nyquist(L); grid on;

% Simulink integration (example sketch):
% 1. Build a Simulink model with 'Plant', 'Controller', and 'Feedback' blocks.
% 2. Use 'linearize' or 'linmod' to obtain the linearized open-loop model L.
% 3. Use 'margin' on the linearized L to verify GM/PM even when the plant
%    is part of a more complex robot model (e.g., using Robotics System Toolbox).
      

In a robotics context, one can build a rigid-body model via rigidBodyTree (Robotics System Toolbox), linearize around a configuration, and then apply the same margin analysis to the joint or Cartesian control loops.

11. Wolfram Mathematica Implementation – Symbolic and Numeric View

Mathematica is well suited for combining symbolic and numeric analysis. Here we implement the same PI-controlled robot joint and compute margins directly from the frequency response.


J  = 0.01;
b  = 0.1;
Kt = 0.5;

Kp = 20;
Ki = 10;

G[s_] := Kt/(J s^2 + b s);
C[s_] := (Kp s + Ki)/s;
L[s_] := G[s] C[s];

wmin = 0.1;
wmax = 100.0;

(* Bode-like data *)
bodeData = Table[
   {w, 20 Log10[Abs[L[I w]]], Arg[L[I w]] 180/Pi},
   {w, wmin, wmax, 0.1}
];

magPlot = ListLinePlot[
   bodeData[[All, {1, 2}]],
   Joined -> True,
   AxesLabel -> {"omega (rad/s)", "Magnitude (dB)"}
];

phasePlot = ListLinePlot[
   bodeData[[All, {1, 3}]],
   Joined -> True,
   AxesLabel -> {"omega (rad/s)", "Phase (deg)"}
];

(* Gain crossover: |L(j w_gc)| = 1 *)
wgcSol = FindRoot[Abs[L[I w]] == 1, {w, 10.0}];
wgc = w /. wgcSol;
pm = 180 + Arg[L[I wgc]] 180/Pi;

(* Phase crossover: Arg L(j w_pc) = -Pi *)
wpcSol = FindRoot[Arg[L[I w]] == -Pi, {w, 30.0}];
wpc = w /. wpcSol;
gm = 1/Abs[L[I wpc]];
gmDb = 20 Log10[gm];

pmRad = pm Pi/180;
tauMax = pmRad/wgc;
deltaMax = 2 Sin[pmRad/2];

{wgc, wpc, pm, gm, gmDb, tauMax, deltaMax}
      

Mathematica can also be used to manipulate the symbolic expressions for \( L(s) \), \( S(s) \), and \( T(s) \), enabling analytic exploration of how GM/PM and robust stability depend on physical parameters \( J \), \( b \), and \( K_t \).

12. Problems and Solutions

Problem 1 (Additive loop uncertainty and Nyquist distance). Let the nominal loop transfer function \( L(s) \) yield a stable closed loop. Consider the perturbed loop \( \tilde{L}(s) = L(s) + \Delta L(s) \). Show that if

\[ \left| \Delta L(j\omega) \right| < \left| 1 + L(j\omega)\right| \quad \forall \omega \in \mathbb{R}, \]

then the perturbed closed loop is also stable.

Solution. On the imaginary axis the perturbed characteristic function is \( 1 + \tilde{L}(j\omega) = 1 + L(j\omega) + \Delta L(j\omega) \). By the reverse triangle inequality,

\[ \left| 1 + L(j\omega) + \Delta L(j\omega) \right| \geq \left| \left|1+L(j\omega)\right| - \left|\Delta L(j\omega)\right| \right|. \]

If \( \left|\Delta L(j\omega)\right| < \left|1+L(j\omega)\right| \) for all \( \omega \), then the right-hand side is strictly positive, so \( 1 + \tilde{L}(j\omega) \neq 0 \) on the imaginary axis. Since the Nyquist plot of \( L \) does not encircle \( -1 \) and the perturbation cannot cross the circle of radius \( \left|1+L(j\omega)\right| \) about \( -1 \), the perturbed Nyquist curve also avoids \( -1 \) and preserves the closed-loop stability.

Problem 2 (Phase margin and multiplicative uncertainty for \( L(s)=1/(s(s+1)) \)). Consider \( L(s) = 1/(s(s+1)) \). Compute the gain crossover frequency \( \omega_{\mathrm{gc}} \), the phase margin \( \varphi_m \), and the approximate multiplicative uncertainty bound \( \delta_{\max} \approx 2\sin(\varphi_m/2) \).

Solution. The frequency response is

\[ L(j\omega) = \frac{1}{j\omega(1+j\omega)}. \]

The magnitude satisfies

\[ \left|L(j\omega)\right|^2 = \frac{1}{\omega^2(1+\omega^2)}. \]

The gain crossover frequency solves \( \left|L(j\omega_{\mathrm{gc}})\right|=1 \), hence

\[ \omega_{\mathrm{gc}}^2(1+\omega_{\mathrm{gc}}^2) = 1. \]

Setting \( x = \omega_{\mathrm{gc}}^2 \) gives \( x(1+x) = 1 \), i.e. \( x^2 + x - 1 = 0 \), so

\[ \omega_{\mathrm{gc}}^2 = \frac{-1 + \sqrt{5}}{2}, \quad \omega_{\mathrm{gc}} \approx 0.786. \]

The phase of \( L(j\omega) \) can be written as

\[ \arg L(j\omega) = -\left(\frac{\pi}{2} + \arctan\omega\right), \]

since \( j\omega \) contributes \( +\pi/2 \) and \( 1+j\omega \) contributes \( \arctan\omega \). Therefore,

\[ \varphi_m = \pi + \arg L(j\omega_{\mathrm{gc}}) = \pi - \left(\frac{\pi}{2} + \arctan\omega_{\mathrm{gc}}\right) = \frac{\pi}{2} - \arctan\omega_{\mathrm{gc}}. \]

Numerically, \( \omega_{\mathrm{gc}}^2 = (-1+\sqrt{5})/2 \approx 0.618 \), \( \omega_{\mathrm{gc}} \approx 0.786 \), \( \varphi_m \approx 0.905 \,\text{rad} \approx 51.8^\circ \). Hence,

\[ \delta_{\max} \approx 2\sin\frac{\varphi_m}{2} \approx 2\sin\left(\frac{0.905}{2}\right) \approx 0.87. \]

This indicates that the loop can tolerate roughly \( 87\% \) multiplicative uncertainty near the crossover frequency without losing stability.

Problem 3 (Delay margin from phase margin). Let a loop have phase margin \( \varphi_m \) (in radians) at gain crossover frequency \( \omega_{\mathrm{gc}} \). Show that an approximate delay margin is \( \tau_{\max} = \varphi_m / \omega_{\mathrm{gc}} \).

Solution. Adding an extra delay \( \tau_d \) changes the phase at \( \omega_{\mathrm{gc}} \) from \( \arg L(j\omega_{\mathrm{gc}}) \) to \( \arg L(j\omega_{\mathrm{gc}}) - \omega_{\mathrm{gc}} \tau_d \). The loop becomes unstable once the phase reaches \( -\pi \) at unit magnitude. Thus,

\[ \arg L(j\omega_{\mathrm{gc}}) - \omega_{\mathrm{gc}} \tau_d = -\pi. \]

Using the definition \( \varphi_m = \pi + \arg L(j\omega_{\mathrm{gc}}) \), we obtain

\[ \varphi_m = \omega_{\mathrm{gc}} \tau_d, \quad \Rightarrow \quad \tau_d = \frac{\varphi_m}{\omega_{\mathrm{gc}}}. \]

Interpreting this at the boundary of instability yields the approximate delay margin \( \tau_{\max} = \varphi_m / \omega_{\mathrm{gc}} \).

Problem 4 (Required phase margin for a given multiplicative uncertainty). Using the approximation \( \delta_{\max} \approx 2\sin(\varphi_m/2) \), determine the minimal phase margin \( \varphi_m \) (in degrees) that allows a multiplicative uncertainty bound \( \delta_{\mathrm{req}} = 0.5 \) (i.e. \( 50\% \) uncertainty) at the gain crossover frequency.

Solution. We require

\[ 2\sin\frac{\varphi_m}{2} \geq \delta_{\mathrm{req}} = 0.5, \]

i.e.

\[ \sin\frac{\varphi_m}{2} \geq 0.25, \quad \frac{\varphi_m}{2} \geq \arcsin(0.25), \quad \varphi_m \geq 2\arcsin(0.25). \]

Numerically, \( \arcsin(0.25) \approx 14.48^\circ \), so \( \varphi_m \geq 28.96^\circ \). Thus a phase margin of about \( 30^\circ \) is the minimum compatible with this level of uncertainty according to the local approximation.

Problem 5 (Design decision flow from uncertainty bounds). A robot joint controller must tolerate: (i) a possible factor-of-\( 3 \) variation in effective gain due to changing loads; and (ii) up to \( 5 \,\mathrm{ms} \) of additional computational and communication delay. Sketch a decision flow that links these uncertainty bounds to required gain and phase margins, and explain qualitatively how you would retune a PI or lead–lag compensator if the initially computed margins are too small.

Solution (qualitative with flow).

flowchart TD
  S["Start: identify gain and delay uncertainty"] --> G1["Gain variation factor up to 3"]
  S --> D1["Extra delay up to 0.005 s"]
  G1 --> GREQ["Set required gain margin GM_req >= 3"]
  D1 --> PMREQ["From tau_max = phi_m/omega_c infer \nrequired phase margin phi_m_req"]
  GREQ --> DESIGN["Tune controller (e.g. increase lead, reduce bandwidth)"]
  PMREQ --> DESIGN
  DESIGN --> VERIFY["Recompute GM, PM; run simulations"]
  VERIFY --> DONE["If GM, PM >= requirements, implement on robot"]
        

Quantitatively, the gain variation factor of \( 3 \) implies \( G_m \geq 3 \). The delay margin condition \( \tau_{\max} \geq 5 \,\mathrm{ms} \) yields \( \varphi_m \geq \omega_{\mathrm{gc}} \tau_{\max} \). If the initial design fails these checks, one can:

  • Introduce a phase-lead compensator to increase \( \varphi_m \) while keeping bandwidth acceptable.
  • Reduce the loop gain to decrease \( \omega_{\mathrm{gc}} \), trading some speed for more delay margin.
  • Adjust PI gains to moderate integral action, which often reduces high-frequency phase lag.

After retuning, margins are recomputed and evaluated against the uncertainty bounds before implementation on the robot hardware.

13. Summary

In this lesson we interpreted gain and phase margins as geometric distances from the Nyquist locus to the critical point and connected them to explicit inequalities guaranteeing robust stability under loop gain, phase, and time-delay uncertainty. We derived approximate formulas linking phase margin to admissible multiplicative uncertainty and delay margin and showed how gain margins bound constant gain variations. Finally, we implemented margin computations in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, with an eye toward their use in robotic joint control. These constructions provide the classical foundation for the more general robustness and performance analysis of sensitivity functions in subsequent lessons.

14. References

  1. Nyquist, H. (1932). Regeneration theory. Bell System Technical Journal, 11(1), 126–147.
  2. Bode, H. W. (1945). Network analysis and feedback amplifier design. Van Nostrand. (Monograph; foundational for gain/phase margin concepts.)
  3. Truxal, J. G. (1955). Automatic Feedback Control System Synthesis. McGraw–Hill. (Early rigorous treatment of classical margins and design.)
  4. Zames, G. (1966). On the input–output stability of time-varying nonlinear feedback systems. Part I: Conditions using concepts of loop gain, conicity, and positivity. IEEE Transactions on Automatic Control, 11(2), 228–238.
  5. Safonov, M. G., & Athans, M. (1977). Gain and phase margin for multiloop LQG regulators. IEEE Transactions on Automatic Control, 22(2), 173–179.
  6. Doyle, J. C., Stein, G. (1981). Multivariable feedback design: Concepts for a classical/modern synthesis. IEEE Transactions on Automatic Control, 26(1), 4–16.
  7. Maciejowski, J. M. (1989). Multivariable feedback design using classical techniques. IEE Proceedings D (Control Theory and Applications), 136(4), 159–168.
  8. Middleton, R. H., & Goodwin, G. C. (1988). Improved finite word length characteristics in digital control using delta operators. IEEE Transactions on Automatic Control, 33(10), 947–958.
  9. Skogestad, S., & Postlethwaite, I. (1997). Multivariable Feedback Control. Wiley. (Contains rigorous links between classical margins and modern robust control metrics.)
  10. Chen, C.-T. (1984). Linear System Theory and Design. Holt, Rinehart and Winston. (Chapters on frequency-domain stability margins and robustness.)