For a unity-feedback system with loop transfer function \( L(s) = K\,\dfrac{N(s)}{D(s)} \), the closed-loop characteristic equation is

\[ 1 + L(s) = 0 \;\;\Longleftrightarrow\;\; D(s) + K N(s) = 0. \]

Let \( p(s,K) = D(s) + K N(s) \) be the characteristic polynomial in the complex variable \( s \) with real coefficients and real gain \( K \). Stability of the LTI closed-loop system in the time domain requires that all roots of \( p(s,K) \) satisfy

\[ \Re\{s_i(K)\} < 0 \quad\text{for all closed-loop poles } s_i(K). \]

The imaginary axis \( \Re\{s\} = 0 \) is therefore the boundary between asymptotically stable \( \Re\{s\}<0 \) and unstable \( \Re\{s\}>0 \) regions. On a root-locus plot

• Poles in the left-half plane correspond to decaying exponentials.
• Poles on the imaginary axis correspond to undamped oscillations.
• Poles in the right-half plane correspond to growing exponentials (instability).

An imaginary-axis crossing is a value of gain \( K = K_\mathrm{crit} \) and frequency \( \omega_0 > 0 \) such that \( s = \pm j\omega_0 \) are roots of the characteristic equation. These values form the stability boundary in the gain parameter.

flowchart TD
  A["Left-half plane (Re(s) < 0): stable"] --- B["Imaginary axis (Re(s) = 0): stability boundary"]
  B --- C["Right-half plane (Re(s) > 0): unstable"]
  D["Root locus branches"] --> B
  B -->|"K = Kcrit"| E["Undamped oscillation in time response"]
        

Assume \( D(s) \) and \( N(s) \) have real coefficients. To find imaginary-axis crossings, substitute \( s = j\omega \), \( \omega \in \mathbb{R} \), into the characteristic equation

\[ D(j\omega) + K N(j\omega) = 0. \]

Write \( D(j\omega) = D_R(\omega) + j D_I(\omega) \) and \( N(j\omega) = N_R(\omega) + j N_I(\omega) \). Then

\[ \big(D_R + K N_R\big) + j\big(D_I + K N_I\big) = 0 \quad\Longleftrightarrow\quad \begin{cases} D_R(\omega) + K N_R(\omega) = 0, \\ D_I(\omega) + K N_I(\omega) = 0. \end{cases} \]

Eliminating \( K \) by cross-multiplication gives a purely frequency-domain equation for \( \omega \):

\[ D_R(\omega) N_I(\omega) - D_I(\omega) N_R(\omega) = 0. \]

Any real solution \( \omega_0 > 0 \) of this equation produces a candidate crossing frequency. The corresponding gain is

\[ K_\mathrm{crit} = -\frac{D_R(\omega_0)}{N_R(\omega_0)} = -\frac{D_I(\omega_0)}{N_I(\omega_0)} \quad\text{(whichever quotient is numerically well-conditioned).} \]

By construction, \( s = \pm j\omega_0 \) are closed-loop poles for \( K = K_\mathrm{crit} \). On the root locus, this is exactly where a branch crosses the imaginary axis.

From Chapter 7, the Routh–Hurwitz criterion provides necessary and sufficient conditions for all roots of a polynomial to lie in the open left-half plane. Suppose the characteristic polynomial is

\[ p(s,K) = a_n(K) s^n + a_{n-1}(K) s^{n-1} + \cdots + a_0(K), \]

where certain coefficients \( a_i(K) \) depend linearly on the gain \( K \). Construct the Routh table, whose first column elements we denote by \( r_n(K), r_{n-1}(K),\dots,r_0(K) \). For a fixed value of \( K \), the closed-loop system is stable iff

\[ r_i(K) > 0 \quad\text{for all } i. \]

As \( K \) varies, these inequalities define an interval (or, more generally, a union of intervals) of stabilizing gains.

Imaginary-axis crossings and zero rows. For a polynomial with real coefficients:

• Any change in the sign pattern of the first column as \( K \) varies implies that a root has crossed the imaginary axis.
• In the generic case where a pair of complex-conjugate roots crosses the imaginary axis, a row of the Routh table becomes zero at \( K = K_\mathrm{crit} \).

Suppose we have a cubic

\[ p(s,K) = s^3 + \alpha s^2 + \beta s + \gamma(K), \]

with \( \alpha > 0 \), \( \beta > 0 \), and \( \gamma(K) > 0 \) for stabilizing gains. The Routh table is

\[ \begin{array}{c|cc} s^3 & 1 & \beta \\ s^2 & \alpha & \gamma(K) \\ s^1 & \dfrac{\alpha\beta - \gamma(K)}{\alpha} & 0 \\ s^0 & \gamma(K) & 0 \end{array} \]

Stability requires \( \gamma(K) > 0 \) and \( \alpha\beta - \gamma(K) > 0 \). The imaginary-axis crossing occurs when

\[ \alpha\beta - \gamma(K_\mathrm{crit}) = 0 \quad\Longrightarrow\quad \gamma(K_\mathrm{crit}) = \alpha\beta. \]

At this value of \( K \), the \( s^1 \) row becomes zero; the corresponding auxiliary polynomial is formed from the preceding row (here the \( s^2 \) row)

\[ A(s) = \alpha s^2 + \gamma(K_\mathrm{crit}), \]

whose roots lie exactly on the imaginary axis, giving the crossing frequencies.

Consider a unity-feedback loop with open-loop transfer function

\[ L(s) = \frac{K}{s(s+2)(s+4)}, \]

which arises, for example, as a simplified linear model of a DC motor driving a robotic joint where \( s \) represents the Laplace variable for joint position and the factors \( (s+2) \), \( (s+4) \) capture electrical and mechanical damping.

The closed-loop characteristic equation is

\[ 1 + \frac{K}{s(s+2)(s+4)} = 0 \quad\Longleftrightarrow\quad s(s+2)(s+4) + K = 0. \]

Expanding the polynomial:

\[ s(s+2)(s+4) = s^3 + 6s^2 + 8s, \quad\Rightarrow\quad p(s,K) = s^3 + 6s^2 + 8s + K. \]

The Routh table is

\[ \begin{array}{c|cc} s^3 & 1 & 8 \\ s^2 & 6 & K \\ s^1 & \dfrac{6\cdot 8 - 1\cdot K}{6} & 0 \\ s^0 & K & 0 \end{array} = \begin{array}{c|cc} s^3 & 1 & 8 \\ s^2 & 6 & K \\ s^1 & \dfrac{48 - K}{6} & 0 \\ s^0 & K & 0 \end{array}. \]

Stability requires \( K > 0 \) and \( 48 - K > 0 \), hence

\[ 0 < K < 48. \]

The imaginary-axis crossing occurs when the \( s^1 \) row becomes zero:

\[ \frac{48 - K_\mathrm{crit}}{6} = 0 \quad\Longrightarrow\quad K_\mathrm{crit} = 48. \]

The auxiliary polynomial is formed from the \( s^2 \) row at \( K = K_\mathrm{crit} \):

\[ A(s) = 6s^2 + 48 = 0 \quad\Longrightarrow\quad s^2 + 8 = 0 \quad\Longrightarrow\quad s = \pm j\sqrt{8} = \pm j 2\sqrt{2}. \]

Thus, on the root locus, a pair of poles crosses the imaginary axis at \( s = \pm j2\sqrt{2} \) when \( K = 48 \), and the closed-loop system is stable for \( 0 < K < 48 \).

The combined Routh–Hurwitz and root-locus view yields a systematic algorithm for finding imaginary-axis crossings and corresponding stabilizing gains:

1. Derive the closed-loop characteristic polynomial \( p(s,K) \).
2. Form the Routh table symbolically with parameter \( K \).
3. Solve for values of \( K \) where some row element or entire row becomes zero.
4. For each candidate \( K_\mathrm{crit} \), construct the auxiliary polynomial and compute its roots.
5. Retain the values where roots are purely imaginary; these yield the stability boundary.
6. Apply Routh inequalities to determine on which side of each boundary the system remains stable.

flowchart TD
  S["Start with G(s) and H(s)"] --> C["Form p(s,K) = D(s) + K N(s)"]
  C --> R["Build Routh table with parameter K"]
  R --> Z["Solve for K where a Routh element or row is zero"]
  Z --> A["Form auxiliary polynomial from row above the zero row"]
  A --> F["Find roots of auxiliary polynomial"]
  F --> J["Select roots on imag axis -> crossing frequencies"]
  J --> I["Use Routh inequalities to determine stable K intervals"]
        

In Python, the python-control and sympy libraries allow symbolic and numerical analysis of stability boundaries. For our example plant \( L(s) = \dfrac{K}{s(s+2)(s+4)} \), we can:

• Compute the Routh table algebraically for verification.
• Scan over \( K \) and numerically detect when a closed-loop pole crosses the imaginary axis.

import numpy as np
import sympy as sp

# Example: p(s,K) = s^3 + 6 s^2 + 8 s + K
s, K = sp.symbols("s K", real=True)
p = s**3 + 6*s**2 + 8*s + K

def routh_cubic_symbolic(poly, K_symbol):
    """
    Compute the first column of the Routh table for a cubic
    poly(s,K) = s^3 + a1 s^2 + a2 s + a3(K).
    """
    # Extract coefficients in descending powers of s
    a3, a2, a1, a0 = sp.Poly(poly, s).all_coeffs()
    # a3 should be 1 for this example
    row_s3 = [a3, a1]
    row_s2 = [a2, a0]
    b1 = (row_s2[0]*row_s3[1] - row_s3[0]*row_s2[1]) / row_s2[0]
    row_s1 = [sp.simplify(b1), 0]
    row_s0 = [row_s2[1], 0]
    return row_s3, row_s2, row_s1, row_s0

row_s3, row_s2, row_s1, row_s0 = routh_cubic_symbolic(p, K)
print("Routh rows (symbolic):")
print("s^3:", row_s3)
print("s^2:", row_s2)
print("s^1:", sp.simplify(row_s1[0]), 0)
print("s^0:", row_s0[0], 0)

# Solve for K where s^1 row becomes zero
Kcrit = sp.solve(sp.Eq(row_s1[0], 0), K)
print("Critical K (Routh zero row):", Kcrit)

# Auxiliary polynomial at Kcrit
Kcrit_val = Kcrit[0]
A = row_s2[0]*s**2 + row_s2[1].subs(K, Kcrit_val)
print("Auxiliary polynomial:", sp.factor(A))

# Evaluate imaginary-axis roots of A(s)
omega_sq = -sp.factor(A.subs(s, sp.I*sp.symbols("omega"))).as_real_imag()[0] / \
           sp.symbols("omega")**2
print("Crossing frequency squared (from A):", sp.simplify(omega_sq))

# Numerical confirmation with python-control (if installed)
try:
    import control as ctl
    # Define open-loop transfer function G(s) = K / (s (s+2) (s+4)), K will be applied later
    G = ctl.tf([1], [1, 6, 8, 0])

    def max_real_pole(Kval):
        T = ctl.feedback(Kval*G, 1)  # closed-loop transfer function
        poles = ctl.pole(T)
        return np.max(np.real(poles))

    K_vals = np.linspace(0, 80, 161)
    sign_changes = []
    previous = max_real_pole(K_vals[0])
    for Kv in K_vals[1:]:
        current = max_real_pole(Kv)
        if previous < 0 and current > 0:
            sign_changes.append(Kv)
        previous = current

    print("Approximate crossing K from numerical scan:", sign_changes[:3])
except ImportError:
    print("python-control not installed; skipping numeric confirmation.")
      

In a robotics context (e.g., joint position control), one often chooses \( K \) at a fraction of \( K_\mathrm{crit} \) (e.g., \( K \approx 0.3 K_\mathrm{crit} \)) to maintain a safety margin against unmodeled dynamics and parameter uncertainty.

In embedded robotic controllers written in C++ (for instance on top of ROS and Eigen), one often needs a fast stability check for a fixed structure. For the cubic example \( p(s,K) = s^3 + 6s^2 + 8s + K \), the Routh first column is

\[ \begin{bmatrix} 1 \\ 6 \\ \dfrac{48 - K}{6} \\ K \end{bmatrix}, \]

so stability reduces to \( 0 < K < 48 \). A C++ routine can enforce this constraint when tuning a joint-position loop:

#include <iostream>
#include <cmath>

bool isStableGain(double K) {
    // Routh column for p(s,K) = s^3 + 6 s^2 + 8 s + K
    double r3 = 1.0;
    double r2 = 6.0;
    double r1 = (48.0 - K) / 6.0;
    double r0 = K;

    return (r3 > 0.0) && (r2 > 0.0) && (r1 > 0.0) && (r0 > 0.0);
}

int main() {
    for (double K = 0.0; K <= 60.0; K += 10.0) {
        std::cout << "K = " << K
                  << (isStableGain(K) ? " : stable\n" : " : unstable\n");
    }

    // Example policy for a robotic joint:
    double Kcrit = 48.0;
    double Kdesign = 0.4 * Kcrit; // design at 40% of critical gain
    std::cout << "Suggested design gain Kdesign = " << Kdesign
              << (isStableGain(Kdesign) ? " (stable)\n" : " (unstable)\n");
    return 0;
}
      

In a full robotic stack, this code would coexist with Eigen-based state-space models and ROS control loops, but the scalar stability test remains identical.

Java is used in some robotics frameworks (e.g., ROSJava, Java-based simulation). While a full polynomial root solver requires a numerical library (such as EJML or Apache Commons Math), for low-degree polynomials we can exploit the analytic Routh result and still perform a numerical scan for verification.

public class StabilityBoundary {

    // First-column Routh entries for p(s,K) = s^3 + 6 s^2 + 8 s + K
    public static boolean isStable(double K) {
        double r3 = 1.0;
        double r2 = 6.0;
        double r1 = (48.0 - K) / 6.0;
        double r0 = K;
        return (r3 > 0.0) && (r2 > 0.0) && (r1 > 0.0) && (r0 > 0.0);
    }

    public static void main(String[] args) {
        double Kcrit = 48.0;
        System.out.println("Critical gain from Routh: Kcrit = " + Kcrit);

        // Scan around Kcrit to illustrate loss of stability
        for (double K = 30.0; K <= 60.0; K += 5.0) {
            System.out.println("K = " + K + (isStable(K) ? " : stable" : " : unstable"));
        }

        // In a Java-based robotic control framework, K could be tuned online
        // while enforcing this stability constraint.
    }
}
      

In combination with a numerical linear algebra package, one could construct the closed-loop state matrix of a robotic servo and directly check eigenvalues for any chosen gain, while Routh gives symbolic understanding of the stability boundaries.

9.1 MATLAB/Simulink

MATLAB provides direct tools for root locus and stability boundary visualization, which integrate smoothly with Simulink models of robotic joints and actuators.

% Define transfer function G(s) = 1 / (s (s + 2) (s + 4))
s = tf('s');
G = 1 / (s * (s + 2) * (s + 4));

% Root locus of K G(s)
figure;
rlocus(G);
title('Root locus of K / (s (s+2) (s+4))');

% Analytical critical gain from Routh
Kcrit = 48;

% Closed-loop transfer function for Kcrit
K = Kcrit;
Tcrit = feedback(K * G, 1);
poles_crit = pole(Tcrit)

% Check that poles lie on imaginary axis
damp(Tcrit)

% In Simulink:
% - Create a plant block implementing 1 / (s (s+2) (s+4)).
% - Connect a Gain block K, and a unity feedback loop.
% - Use a mask parameter for K and vary it up to Kcrit to observe transition
%   from overdamped to sustained oscillations at K = Kcrit.
      

9.2 Wolfram Mathematica

Mathematica can symbolically derive stability boundaries and visualize root loci for control of electromechanical systems and robot joints.

(* Define symbols *)
Clear[s, k];
s = ComplexExpand@Symbol["s"];
k = Symbol["k"];

(* Characteristic polynomial for example *)
p[s_, k_] := s^3 + 6 s^2 + 8 s + k;

(* Routh-based critical gain: solve (48 - k)/6 == 0 *)
kcrit = Solve[(48 - k)/6 == 0, k][[1, 1, 2]]

(* Auxiliary polynomial at kcrit *)
aux[s_] = 6 s^2 + k /. k -> kcrit

(* Imaginary-axis roots *)
omega = Symbol["omega"];
eq = aux[I omega] == 0;
solOmega = Solve[eq, omega]

(* Transfer function model and root locus *)
Needs["ControlSystems`"];
G = TransferFunctionModel[1/(s (s + 2) (s + 4)), s];

RootLocusPlot[G, {k, 0, 80},
  PlotLabel -> "Root locus of K / (s (s+2) (s+4))"];

(* Evaluate closed-loop poles at kcrit *)
Tcrit = SystemsModelFeedback[kcrit G, 1];
Eigenvalues[Tcrit]
      

These symbolic tools are particularly useful for verifying hand-derived stability boundaries before embedding numerical controllers in a robotic system.

Problem 1 (Routh-based stability boundary for a cubic): Consider the characteristic polynomial \( p(s,K) = s^3 + 3s^2 + 2s + K \).
(a) Use the Routh–Hurwitz criterion to find the range of \( K \) for which the system is asymptotically stable.
(b) Determine the value(s) of \( K \) at which the closed-loop poles lie on the imaginary axis and find the corresponding oscillation frequency.

Solution:

(a) Construct the Routh table:

\[ \begin{array}{c|cc} s^3 & 1 & 2 \\ s^2 & 3 & K \\ s^1 & \dfrac{3\cdot 2 - 1\cdot K}{3} & 0 \\ s^0 & K & 0 \end{array} = \begin{array}{c|cc} s^3 & 1 & 2 \\ s^2 & 3 & K \\ s^1 & \dfrac{6 - K}{3} & 0 \\ s^0 & K & 0 \end{array}. \]

Stability requires all first-column elements positive:

\[ 1 > 0,\quad 3 > 0,\quad \frac{6 - K}{3} > 0,\quad K > 0 \;\;\Longrightarrow\;\; 0 < K < 6. \]

(b) Imaginary-axis crossing occurs when the \( s^1 \) row becomes zero:

\[ \frac{6 - K_\mathrm{crit}}{3} = 0 \quad\Longrightarrow\quad K_\mathrm{crit} = 6. \]

The auxiliary polynomial is obtained from the \( s^2 \) row:

\[ A(s) = 3s^2 + 6 = 0 \quad\Longrightarrow\quad s^2 + 2 = 0 \quad\Longrightarrow\quad s = \pm j\sqrt{2}. \]

Therefore, the crossing frequency is \( \omega_0 = \sqrt{2} \), and the system is stable for \( 0 < K < 6 \).

Problem 2 (Imaginary-axis crossing via \( s = j\omega \) substitution): For the same polynomial \( p(s,K) = s^3 + 3s^2 + 2s + K \), verify the result of Problem 1 using the substitution \( s = j\omega \).

Solution:

Substitute \( s = j\omega \):

\[ p(j\omega,K) = -j\omega^3 - 3\omega^2 + 2j\omega + K = (K - 3\omega^2) + j(-\omega^3 + 2\omega). \]

Separating real and imaginary parts,

\[ \begin{cases} K - 3\omega^2 = 0, \\ -\omega^3 + 2\omega = 0. \end{cases} \]

From the imaginary part, \( \omega(-\omega^2 + 2) = 0 \), so either \( \omega = 0 \) or \( \omega^2 = 2 \). The nontrivial oscillation corresponds to \( \omega_0^2 = 2 \). Substituting into the real equation,

\[ K_\mathrm{crit} - 3\omega_0^2 = 0 \quad\Longrightarrow\quad K_\mathrm{crit} = 3\cdot 2 = 6. \]

This matches the Routh result, confirming that the imaginary-axis crossing with \( s = \pm j\sqrt{2} \) occurs at \( K = 6 \).

Problem 3 (Stability range for a robotic joint-like plant): Consider the characteristic polynomial \( p(s,K) = s^3 + 14s^2 + 40s + K \), which arises from a plant \( L(s) = \dfrac{K}{s(s+4)(s+10)} \).
(a) Determine the range of \( K \) for which the closed-loop system is stable.
(b) Find the critical gain and imaginary-axis crossing frequency.

Solution:

The Routh table is

\[ \begin{array}{c|cc} s^3 & 1 & 40 \\ s^2 & 14 & K \\ s^1 & \dfrac{14\cdot 40 - K}{14} & 0 \\ s^0 & K & 0 \end{array} = \begin{array}{c|cc} s^3 & 1 & 40 \\ s^2 & 14 & K \\ s^1 & \dfrac{560 - K}{14} & 0 \\ s^0 & K & 0 \end{array}. \]

Stability requires \( K > 0 \) and \( 560 - K > 0 \), hence

\[ 0 < K < 560. \]

The imaginary-axis crossing occurs when the \( s^1 \) row is zero:

\[ 560 - K_\mathrm{crit} = 0 \quad\Longrightarrow\quad K_\mathrm{crit} = 560. \]

The auxiliary polynomial is

\[ A(s) = 14s^2 + 560 = 0 \quad\Longrightarrow\quad s^2 + 40 = 0 \quad\Longrightarrow\quad s = \pm j\sqrt{40}. \]

Thus the oscillation frequency at the stability boundary is \( \omega_0 = \sqrt{40} \), and the system is stable for \( 0 < K < 560 \).

Problem 4 (Decision flow for selecting a safe gain): Sketch a decision process for choosing a gain \( K \) that satisfies: (i) asymptotic stability, and (ii) a safety margin against uncertainty (e.g., choose \( K \) not too close to \( K_\mathrm{crit} \)). Illustrate the logic in a flow diagram.

Solution (flow):

flowchart TD
  ST["Start with characteristic p(s,K)"] --> R["Use Routh to find stability range (0 < K < Kcrit)"]
  R --> P1["Pick nominal K within range"]
  P1 --> M["Is K far enough from Kcrit? (e.g. K < 0.5*Kcrit)"]
  M -->|yes| OK["Accept K for implementation"]
  M -->|no| ADJ["Reduce K to increase margin"]
  ADJ --> OK
        

In this lesson we treated the imaginary axis as the stability boundary in the complex plane and showed how root-locus branches cross this boundary as the loop gain varies. Analytically, the condition for imaginary-axis crossing is obtained either by substitution \( s = j\omega \) and separating real and imaginary parts, or via the Routh–Hurwitz criterion by detecting zero rows in the Routh table. We demonstrated how to compute critical gains and crossing frequencies for cubic characteristic equations and linked these results to software implementations in Python, C++, Java, MATLAB/Simulink, and Mathematica, with an eye toward robotic joint and servo control. These stability boundaries will be crucial in subsequent design lessons, where we shape the root locus to satisfy transient and steady-state specifications while maintaining adequate robustness margins.