Consider a single-input single-output (SISO) LTI plant whose transfer function depends on a vector of real parameters \( \mathbf{p} = [p_1,\dots,p_q]^\top \) (masses, inertias, resistances, gains, etc.):

\[ G(s,\mathbf{p}) = \frac{B(s,\mathbf{p})}{A(s,\mathbf{p})} = \frac{\sum_{k=0}^{m} b_k(\mathbf{p}) s^k} {\sum_{k=0}^{n} a_k(\mathbf{p}) s^k}, \quad a_n(\mathbf{p}) \neq 0. \]

A nominal parameter vector \( \mathbf{p}_0 \) is obtained from identification, datasheets, or modeling assumptions. Manufacturing tolerances and operating-condition variations imply that the true parameter vector lies in a set \( \mathcal{P} \subset \mathbb{R}^q \):

\[ \mathbf{p} = \mathbf{p}_0 + \Delta \mathbf{p}, \qquad \Delta \mathbf{p} \in \mathcal{D}, \]

where \( \mathcal{D} \) is a bounded set of admissible deviations (typically an axis-aligned box or an ellipsoid). A common description for independent bounded relative errors is

\[ p_i = p_{i,0}(1 + \delta_i), \qquad |\delta_i| \le \bar{\delta}_i, \quad i = 1,\dots,q. \]

The corresponding family of plants is \( \{ G(s,\mathbf{p}) : \mathbf{p} \in \mathcal{P} \} \). Even though the model structure (order \(n\), input/output) is fixed, the pole set and zero set vary with \( \mathbf{p} \).

For classical linear control, we mainly need to understand how:

• Parametric uncertainty moves poles and zeros in the complex plane.
• These movements may push poles across the imaginary axis, creating potential loss of stability.
• Small parameter changes can induce large changes in closed-loop performance if the poles are close to the imaginary axis.

flowchart TD
  P0["Nominal parameters p0"] --> DSET["Uncertainty set D (delta p)"]
  DSET --> FAM["Family of models G(s,p)"]
  FAM --> ROOTS["Pole/zero sets for all p in P"]
  ROOTS --> STAB["Stability: do all poles stay in left half-plane?"]
  ROOTS --> PERF["Performance: damping, overshoot, bandwidth"]
        

Start from an ordinary differential equation whose coefficients depend on physical parameters. For a generic \(n\)-th order LTI system:

\[ a_n(\mathbf{p}) \, y^{(n)}(t) + a_{n-1}(\mathbf{p}) \, y^{(n-1)}(t) + \dots + a_0(\mathbf{p}) \, y(t) = b_m(\mathbf{p}) \, u^{(m)}(t) + \dots + b_0(\mathbf{p}) \, u(t), \]

the Laplace-domain transfer function from input \( u \) to output \( y \) is

\[ G(s,\mathbf{p}) = \frac{Y(s)}{U(s)} = \frac{b_m(\mathbf{p}) s^m + \dots + b_0(\mathbf{p})} {a_n(\mathbf{p}) s^n + \dots + a_0(\mathbf{p})}. \]

We assume that each coefficient is a smooth (often affine) function of parameters:

\[ a_k(\mathbf{p}) = a_{k,0} + \sum_{j=1}^{q} \alpha_{k j} p_j, \qquad b_\ell(\mathbf{p}) = b_{\ell,0} + \sum_{j=1}^{q} \beta_{\ell j} p_j. \]

For small deviations \( \Delta \mathbf{p} \) around the nominal \( \mathbf{p}_0 \), a first-order Taylor expansion gives

\[ G(s,\mathbf{p}_0 + \Delta \mathbf{p}) \approx G(s,\mathbf{p}_0) + \sum_{j=1}^{q} \frac{\partial G}{\partial p_j}(s,\mathbf{p}_0) \, \Delta p_j. \]

This linearization in parameter space is fundamental: it allows us to quantify how sensitive poles and zeros are to each parameter without recomputing the roots for every possible \( \mathbf{p} \).

As a concrete example relevant to robotic joints or actuated mechanical links, consider the rotational mass-spring-damper model with torque input \( \tau(t) \) and angle output \( \theta(t) \):

\[ J \, \ddot{\theta}(t) + b \, \dot{\theta}(t) + k \, \theta(t) = \tau(t), \]

where \(J\) is the rotary inertia, \(b\) viscous friction, and \(k\) joint stiffness. The transfer function from torque to angle is

\[ G(s;J,b,k) = \frac{\Theta(s)}{\tau(s)} = \frac{1}{J s^2 + b s + k}. \]

Each of \(J, b, k\) is uncertain, leading to a family of second-order transfer functions with the same structure but different poles.

For a fixed parameter vector \( \mathbf{p} \), the poles of \( G(s,\mathbf{p}) \) are the roots of the characteristic polynomial

\[ A(s,\mathbf{p}) = a_n(\mathbf{p}) s^n + \dots + a_0(\mathbf{p}). \]

Let \( \lambda_i(\mathbf{p}) \) denote a simple pole (root) of \( A(s,\mathbf{p}) \), i.e.,

\[ A(\lambda_i(\mathbf{p}), \mathbf{p}) = 0, \qquad \frac{\partial A}{\partial s}(\lambda_i(\mathbf{p}),\mathbf{p}) \neq 0. \]

Applying the implicit function theorem and differentiating the identity \(A(\lambda_i(\mathbf{p}),\mathbf{p}) = 0\) with respect to a parameter \(p_j\), we obtain

\[ \frac{\partial A}{\partial s}(\lambda_i,\mathbf{p}) \, \frac{\partial \lambda_i}{\partial p_j}(\mathbf{p}) + \frac{\partial A}{\partial p_j}(\lambda_i,\mathbf{p}) = 0, \]

hence the pole sensitivity is

\[ \frac{\partial \lambda_i}{\partial p_j}(\mathbf{p}) = - \frac{\dfrac{\partial A}{\partial p_j}(\lambda_i,\mathbf{p})} {\dfrac{\partial A}{\partial s}(\lambda_i,\mathbf{p})}. \]

For small deviations \( \Delta p_j \), the pole variation is approximated by

\[ \lambda_i(\mathbf{p}_0 + \Delta \mathbf{p}) \approx \lambda_i(\mathbf{p}_0) + \sum_{j=1}^{q} \frac{\partial \lambda_i}{\partial p_j}(\mathbf{p}_0) \, \Delta p_j. \]

Zeros behave analogously, but they are roots of the numerator polynomial \(B(s,\mathbf{p})\), and their sensitivity is obtained by replacing \(A\) with \(B\) in the formula above.

Intuitively:

• If a pole is well inside the left half-plane, its real part is typically only mildly affected by small parameter changes.
• If a pole is close to the imaginary axis, even small \( \Delta p_j \) may push it into instability, which motivates conservative design and robustness margins.

flowchart TD
  DP["Parameter change delta p_j"] --> DA["Change in coefficients A(s,p)"]
  DA --> DL["Pole movement approx: d lambda_i / d p_j"]
  DL --> SPlane["New pole location in s-plane"]
  SPlane --> IMPL1["Stability check (Re(lambda) < 0?)"]
  SPlane --> IMPL2["Time response: damping, overshoot, speed"]
        

Returning to the rotational joint model \( J \ddot{\theta}(t) + b \dot{\theta}(t) + k \theta(t) = \tau(t) \), the poles of \( G(s;J,b,k) = 1/(J s^2 + b s + k) \) are the roots of

\[ J s^2 + b s + k = 0. \]

By the quadratic formula,

\[ \lambda_{1,2}(J,b,k) = \frac{-b \pm \sqrt{b^2 - 4 J k}}{2J}. \]

Introducing the standard second-order parameters \( \omega_n = \sqrt{k/J} \) and \( \zeta = \dfrac{b}{2\sqrt{kJ}} \), we can rewrite the poles as

\[ \lambda_{1,2}(\zeta,\omega_n) = -\zeta \,\omega_n \pm j \,\omega_n \sqrt{1 - \zeta^2}. \]

Suppose \(J, b, k\) are uncertain but small relative variations keep the system underdamped (\(\zeta < 1\)). Then:

• Increasing \(J\) (heavier link) decreases \(\omega_n\) and thus slows down the response.
• Increasing \(b\) (more damping) increases \(\zeta\), making poles move further into the left half-plane and reducing overshoot.
• Increasing \(k\) (stiffer joint) increases both \(\omega_n\) and slightly decreases \(\zeta\) (for fixed \(b,J\)), pulling poles farther from the origin but closer to the imaginary axis (potentially increasing overshoot).

For a single parameter, say \(J\), we can explicitly differentiate:

\[ \omega_n(J) = \sqrt{\frac{k}{J}}, \qquad \frac{d\omega_n}{dJ} = -\frac{1}{2} \sqrt{\frac{k}{J^3}}. \]

Keeping \(b,k\) fixed, the damping ratio is

\[ \zeta(J) = \frac{b}{2\sqrt{kJ}}, \qquad \frac{d\zeta}{dJ} = -\frac{b}{4} \frac{1}{\sqrt{k}} J^{-3/2}. \]

Combining these derivatives with \( \lambda_{1,2}(\zeta,\omega_n) \), we can compute \( \dfrac{d\lambda_{1,2}}{dJ} \) via the chain rule, obtaining an analytic description of how the joint poles move with inertia uncertainty. Similar expressions can be derived for \(b\) and \(k\).

In practical robotic control, one often has intervals \(J \in [J_{\min}, J_{\max}]\), \(b \in [b_{\min}, b_{\max}]\), etc. The goal is then to design a controller such that all closed-loop poles remain in a desirable region of the \(s\)-plane for every admissible combination of these parameters.

We now illustrate how to explore parametric uncertainty and pole variations with different software environments commonly used in control and robotics.

5.1 Python Example (with python-control and Robotics Context)

The following Python script uses numpy and python-control to sample uncertain parameters for a robotic joint and scatter the corresponding poles. Such scripts are typically integrated into larger robotics projects (for example, ROS-based simulation using ros_control and roboticstoolbox-python), but we keep the example minimal here.

import numpy as np
import matplotlib.pyplot as plt
from control import TransferFunction, step_response

# Nominal joint parameters (e.g., small robot arm link)
J0 = 0.01   # kg m^2
b0 = 0.05   # N m s/rad
k0 = 2.0    # N m/rad

# Relative uncertainty bounds (+/- 20%)
rel_unc = 0.2
N = 200  # number of Monte Carlo samples

def sample_params():
    def sample_one(p0):
        # uniform in [p0*(1-rel_unc), p0*(1+rel_unc)]
        return p0 * (1.0 + rel_unc * (2.0*np.random.rand() - 1.0))
    J = sample_one(J0)
    b = sample_one(b0)
    k = sample_one(k0)
    return J, b, k

poles_real = []
poles_imag = []

for _ in range(N):
    J, b, k = sample_params()
    # Second-order denominator: J s^2 + b s + k
    den = [J, b, k]
    num = [1.0]
    G = TransferFunction(num, den)
    lam = np.roots(den)  # open-loop poles of the joint dynamics
    for p in lam:
        poles_real.append(np.real(p))
        poles_imag.append(np.imag(p))

# Plot pole cloud
plt.figure()
plt.scatter(poles_real, poles_imag, marker='x')
plt.axvline(x=0.0, linestyle='--')  # imaginary axis
plt.xlabel("Re(s)")
plt.ylabel("Im(s)")
plt.title("Joint pole locations under parametric uncertainty")

# Example step response for nominal model
G_nom = TransferFunction([1.0], [J0, b0, k0])
t, y = step_response(G_nom)
plt.figure()
plt.plot(t, y)
plt.xlabel("t [s]")
plt.ylabel("theta(t) [rad]")
plt.title("Nominal joint step response")
plt.show()
      

5.2 C++ Example (Analytic Second-Order Poles, Robotics-Oriented)

The C++ snippet below computes the poles of the same second-order model for randomly perturbed parameters. In a robotics application, such logic could be embedded in a ROS node (e.g., within a ros_control controller) to assess robustness of joint-level controllers.

#include <iostream>
#include <complex>
#include <random>

std::pair<std::complex<double>, std::complex<double>>
joint_poles(double J, double b, double k)
{
    // Solve J s^2 + b s + k = 0
    std::complex<double> disc = std::complex<double>(b*b - 4.0*J*k, 0.0);
    std::complex<double> sqrt_disc = std::sqrt(disc);
    std::complex<double> s1 = (-b + sqrt_disc) / (2.0*J);
    std::complex<double> s2 = (-b - sqrt_disc) / (2.0*J);
    return {s1, s2};
}

int main()
{
    double J0 = 0.01;
    double b0 = 0.05;
    double k0 = 2.0;
    double rel_unc = 0.2;

    std::mt19937 gen(42);
    std::uniform_real_distribution<double> uni(-rel_unc, rel_unc);

    for (int i = 0; i < 10; ++i) {
        double J = J0 * (1.0 + uni(gen));
        double b = b0 * (1.0 + uni(gen));
        double k = k0 * (1.0 + uni(gen));
        auto poles = joint_poles(J, b, k);
        std::cout << "Sample " << i
                  << ": poles = " << poles.first
                  << ", " << poles.second << std::endl;
    }

    return 0;
}
      

5.3 Java Example (Using Apache Commons Math Complex Numbers)

Java is used in some robotics platforms (for example, educational robotics or Android-based control interfaces). The following Java class uses org.apache.commons.math3.complex.Complex to compute poles.

import org.apache.commons.math3.complex.Complex;

public class JointPoleUncertainty {

    public static Complex[] poles(double J, double b, double k) {
        // J s^2 + b s + k = 0
        double disc = b*b - 4.0*J*k;
        Complex sqrtDisc = new Complex(disc, 0.0).sqrt();
        Complex s1 = new Complex(-b, 0.0).add(sqrtDisc).divide(2.0*J);
        Complex s2 = new Complex(-b, 0.0).subtract(sqrtDisc).divide(2.0*J);
        return new Complex[]{s1, s2};
    }

    public static void main(String[] args) {
        double J0 = 0.01;
        double b0 = 0.05;
        double k0 = 2.0;
        double rel_unc = 0.2;

        java.util.Random rng = new java.util.Random(42L);

        for (int i = 0; i < 10; ++i) {
            double J = J0 * (1.0 + rel_unc * (2.0*rng.nextDouble() - 1.0));
            double b = b0 * (1.0 + rel_unc * (2.0*rng.nextDouble() - 1.0));
            double k = k0 * (1.0 + rel_unc * (2.0*rng.nextDouble() - 1.0));
            Complex[] ps = poles(J, b, k);
            System.out.println("Sample " + i
                + " poles: " + ps[0] + ", " + ps[1]);
        }
    }
}
      

5.4 MATLAB/Simulink Example (Control System Toolbox and Robotics Context)

MATLAB and Simulink are widely used in robotics and mechatronics. Below we use the Control System Toolbox to sweep parameters and observe pole movement. In a robotics workflow, one could combine this with models from the Robotics System Toolbox (e.g. linearizing a rigidBodyTree model around an operating point).

J0 = 0.01;  b0 = 0.05;  k0 = 2.0;
rel_unc = 0.2;
N = 200;

poles_all = [];

for i = 1:N
    J = J0 * (1 + rel_unc*(2*rand - 1));
    b = b0 * (1 + rel_unc*(2*rand - 1));
    k = k0 * (1 + rel_unc*(2*rand - 1));

    sys = tf(1, [J b k]);   % G(s) = 1/(J s^2 + b s + k)
    poles_all = [poles_all; pole(sys).'];
end

figure; plot(real(poles_all), imag(poles_all), 'x');
xlabel('Re(s)'); ylabel('Im(s)');
title('Pole cloud under parametric uncertainty'); grid on;

% Nominal step response (e.g., joint angle under torque command)
sys_nom = tf(1, [J0 b0 k0]);
figure; step(sys_nom);
title('Nominal joint step response');
      

5.5 Wolfram Mathematica Example (Symbolic Pole Sensitivities)

Mathematica can symbolically compute pole sensitivities for our second-order model.

(* Define symbolic parameters and characteristic polynomial *)
Clear[J, b, k, s];
charPoly = J s^2 + b s + k;

(* Solve for poles symbolically *)
lambdaSol = Solve[charPoly == 0, s];
lambda1 = s /. lambdaSol[[1]];
lambda2 = s /. lambdaSol[[2]];

(* Sensitivity of the first pole with respect to J *)
dLambda1dJ = D[lambda1, J] // FullSimplify

(* Substitute nominal values to obtain numerical sensitivity *)
lambda1Num = lambda1 /. {J -> 0.01, b -> 0.05, k -> 2.0};
dLambda1dJNum = dLambda1dJ /. {J -> 0.01, b -> 0.05, k -> 2.0} // N
      

These symbolic expressions can be used to approximate pole movement due to small parameter changes, in line with the Taylor expansion discussed earlier.

Problem 1 (Pole Sensitivity of a Simple Third-Order Polynomial): Consider the characteristic polynomial \( A(s,a) = s^3 + (3 + a) s^2 + 4 s + 1 \) where \(a\) is a real parameter. Let \( \lambda(a) \) be a simple pole satisfying \( A(\lambda(a),a) = 0 \). Derive an expression for \( \dfrac{d\lambda}{da} \) at a generic value of \(a\).

Solution:

We have

\[ A(s,a) = s^3 + (3 + a) s^2 + 4 s + 1. \]

The partial derivatives are

\[ \frac{\partial A}{\partial s}(s,a) = 3 s^2 + 2(3 + a) s + 4, \qquad \frac{\partial A}{\partial a}(s,a) = s^2. \]

Using the general pole sensitivity formula,

\[ \frac{d\lambda}{da} = - \frac{\dfrac{\partial A}{\partial a}(\lambda,a)} {\dfrac{\partial A}{\partial s}(\lambda,a)} = - \frac{\lambda(a)^2} {3 \lambda(a)^2 + 2(3 + a)\lambda(a) + 4}. \]

This gives the first-order change of the pole with respect to the parameter \(a\).

Problem 2 (Second-Order System with Uncertain Natural Frequency): Consider the standard second-order transfer function \( G(s) = \dfrac{\omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2} \) with fixed damping ratio \(0 < \zeta < 1\) and uncertain natural frequency \( \omega_n \in [\omega_{\min}, \omega_{\max}] \). Determine the region of possible pole locations in the complex plane.

Solution:

The poles are

\[ \lambda_{1,2}(\omega_n) = -\zeta \omega_n \pm j \,\omega_n \sqrt{1 - \zeta^2}. \]

As \( \omega_n \) varies in the interval \( [\omega_{\min}, \omega_{\max}] \), the real and imaginary parts of the poles are

\[ \operatorname{Re}(\lambda_{1,2}) = -\zeta \omega_n \in [-\zeta \omega_{\max}, -\zeta \omega_{\min}], \] \[ |\operatorname{Im}(\lambda_{1,2})| = \omega_n \sqrt{1 - \zeta^2} \in [\omega_{\min}\sqrt{1 - \zeta^2}, \, \omega_{\max}\sqrt{1 - \zeta^2} ]. \]

Thus the two complex-conjugate poles move along a straight line through the origin with slope \( \pm \sqrt{1 - \zeta^2} \), between the points corresponding to \( \omega_{\min} \) and \( \omega_{\max} \). The possible pole locations form a line segment in the left half-plane for each of the two complex-conjugate poles.

Problem 3 (Robotic Joint with Inertia Uncertainty): For the joint model \( G(s;J,b,k) = 1/(J s^2 + b s + k) \) with fixed \(b,k > 0\), suppose \(J\) is uncertain in the interval \( [J_0 (1 - \alpha), J_0 (1 + \alpha)] \) with \(0 < \alpha < 1\). Show that the real parts of the poles satisfy \( \operatorname{Re}(\lambda_{1,2}) \le -\dfrac{b}{2 J_0 (1 + \alpha)} \).

Solution:

The poles are

\[ \lambda_{1,2}(J) = \frac{-b \pm \sqrt{b^2 - 4 J k}}{2J}. \]

For real \(J > 0\), both poles have real part

\[ \operatorname{Re}(\lambda_{1,2}(J)) = -\frac{b}{2J}. \]

Since \(J\) is in the interval \( [J_0 (1 - \alpha), J_0 (1 + \alpha)] \), the smallest (most negative) real part occurs at the smallest \(J\), and the largest (least negative) real part occurs at the largest \(J\). Therefore

\[ \operatorname{Re}(\lambda_{1,2}(J)) \in \left[ -\frac{b}{2 J_0 (1 - \alpha)}, \; -\frac{b}{2 J_0 (1 + \alpha)} \right]. \]

Hence each pole always satisfies \( \operatorname{Re}(\lambda_{1,2}) \le -\dfrac{b}{2 J_0 (1 + \alpha)} \), which provides a simple conservative lower bound on how close the poles can approach the imaginary axis.

Problem 4 (Closed-Loop Sensitivity to Plant Gain): A unit-feedback loop has plant \( G(s,K) = \dfrac{K}{s(s + 1)} \) and controller \(C(s) = 1\). The closed-loop characteristic polynomial is \( s^2 + s + K \). Using the sensitivity formula, compute \( \dfrac{d\lambda}{dK} \) for a closed-loop pole \( \lambda(K) \) and discuss qualitatively how increasing \(K\) moves the poles.

Solution:

The characteristic polynomial is \( A(s,K) = s^2 + s + K \). Then

\[ \frac{\partial A}{\partial s}(s,K) = 2s + 1, \qquad \frac{\partial A}{\partial K}(s,K) = 1. \]

For a simple root \( \lambda(K) \), the sensitivity is

\[ \frac{d\lambda}{dK} = - \frac{\dfrac{\partial A}{\partial K}(\lambda,K)} {\dfrac{\partial A}{\partial s}(\lambda,K)} = - \frac{1}{2\lambda(K) + 1}. \]

When the pole lies in the left half-plane with negative real part, the denominator \(2\lambda(K) + 1\) has negative real part, so \(d\lambda/dK\) tends to have positive real part. This indicates that increasing \(K\) tends to move the poles toward the right half-plane, consistent with the classical root-locus behavior: higher loop gain pulls poles toward zeros, potentially reducing stability margins.

In this lesson we modeled parametric uncertainty in LTI systems by allowing physical parameters to vary within bounded sets, inducing families of transfer functions with varying poles and zeros. Using polynomial-based sensitivity analysis, we derived analytic expressions for how poles move with respect to parameters and applied these to second-order mechanical/robotic joint models. We then demonstrated numerical exploration techniques across Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, showing how engineers can visualize pole clouds and assess robustness of candidate control designs. These concepts form the basis for more advanced robustness analysis in subsequent lessons.