Consider a unity-feedback loop with controller \( C(s,\boldsymbol{\alpha}) \) and plant \( G(s,\boldsymbol{\theta}) \) parameterized by vectors \( \boldsymbol{\alpha} \in \mathbb{R}^{p} \) (controller parameters) and \( \boldsymbol{\theta} \in \mathbb{R}^{q} \) (plant parameters). The loop transfer function is

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

Robustness questions are posed over a set of admissible parameters \( \mathcal{P} \subset \mathbb{R}^{p+q} \), for example a hyper-rectangle:

\[ \mathcal{P} = \Theta_{\text{plant}} \times \Theta_{\text{ctrl}} = \Theta_1 \times \cdots \times \Theta_{p+q}, \quad \Theta_i = [\theta_i^{\min},\theta_i^{\max}]. \]

For each \( (\boldsymbol{\alpha},\boldsymbol{\theta}) \in \mathcal{P} \), stability is determined by the roots of the characteristic polynomial \( D(s,\boldsymbol{\alpha},\boldsymbol{\theta}) \) in the denominator of \( T(s,\boldsymbol{\alpha},\boldsymbol{\theta}) \). Denote its roots by \( s_i(\boldsymbol{\alpha},\boldsymbol{\theta}) \). Internal stability requires

\[ \operatorname{Re}\big(s_i(\boldsymbol{\alpha},\boldsymbol{\theta})\big) < 0 \quad \forall i. \]

Directly verifying this for all points in a continuous set \( \mathcal{P} \) is typically intractable. Parameter sweeps approximate this by evaluating stability and performance on a finite grid or finite sample \( \{(\boldsymbol{\alpha}^{(k)},\boldsymbol{\theta}^{(k)})\}_{k=1}^M \subset \mathcal{P} \).

A single revolute robot joint with inertia \( J \) and viscous friction \( B \) can be modeled (neglecting gravity) as

\[ G(s;J,B) = \frac{\Theta(s)}{U(s)} = \frac{1}{J s^2 + B s}, \]

where \( U(s) \) is the motor torque input and \( \Theta(s) \) is the joint angle. Consider a PD controller

\[ C(s;K_p,K_d) = K_p + K_d s. \]

With unity feedback, the loop transfer function is \( L(s) = C(s)\,G(s) \) and the closed-loop characteristic polynomial is obtained from \( 1 + L(s) = 0 \). After clearing denominators,

\[ J s^2 + B s + K_p + K_d s = 0 \quad \Rightarrow \quad J s^2 + (B + K_d) s + K_p = 0. \]

This is a standard second-order system. Identifying with \( s^2 + 2\zeta\omega_n s + \omega_n^2 = 0 \) gives

\[ \omega_n(J,K_p) = \sqrt{\frac{K_p}{J}}, \qquad \zeta(J,B,K_p,K_d) = \frac{B + K_d}{2 \sqrt{J K_p}}. \]

For this second-order case, Routh–Hurwitz yields the simple stability conditions

\[ J > 0, \quad B + K_d > 0, \quad K_p > 0. \]

Many robot arms operate in ranges where \( J \) and \( B \) are uncertain but bounded. A PD design is robustly stable over a box \( [J^{\min},J^{\max}] \times [B^{\min},B^{\max}] \) if the above inequalities hold throughout the box. However, robustness in performance (e.g., overshoot, settling time) requires additional analysis, well suited to software sweeps.

Given a finite sample \( (\boldsymbol{\alpha}^{(k)},\boldsymbol{\theta}^{(k)}) \), software can evaluate:

• Stability: all poles of \( T(s) \) in the left half-plane.
• Time-domain metrics for step response: damping ratio \( \zeta \), overshoot \( M_p \), settling time \( t_s \), etc.
• Frequency-domain metrics: gain margin, phase margin, crossover frequencies, bandwidth.

For an underdamped second-order closed-loop system (\( 0 < \zeta < 1 \)), the standard formulas (see previous chapters) give

\[ M_p(\zeta) = \exp\!\left(-\frac{\pi \zeta}{\sqrt{1 - \zeta^2}}\right), \qquad t_s(\zeta,\omega_n) \approx \frac{4}{\zeta\,\omega_n}. \]

When these are treated as functions of the parameters via \( \zeta(J,B,K_p,K_d) \) and \( \omega_n(J,K_p) \), we can define a robust performance region:

\[ \mathcal{R} = \left\{ (J,B,K_p,K_d) \in \mathcal{P} : \\ \zeta(J,B,K_p,K_d) \ge \zeta_{\min},\; M_p\big(\zeta(J,B,K_p,K_d)\big) \le M_p^{\max},\; t_s(J,B,K_p,K_d) \le t_s^{\max} \right\}. \]

A parameter sweep numerically approximates \( \mathcal{R} \) by evaluating these inequalities on a discrete grid.

At a high level, a parameter sweep and robustness study can be summarized as:

flowchart TD
  A["Define plant G(s,theta) and controller C(s,alpha)"] --> B["Specify parameter ranges for theta and alpha"]
  B --> C["Choose sampling strategy: grid or random"]
  C --> D["For each sample: build closed-loop model"]
  D --> E["Compute poles, step metrics, margins"]
  E --> F["Classify: stable / unstable, spec satisfied?"]
  F --> G["Aggregate results: maps, histograms, worst cases"]
  G --> H["Refine controller design or parameter ranges"]
        

Software environments differ mainly in how closed-loop systems are constructed and how performance metrics are extracted. The mathematical concepts are identical across tools.

Parameter sweeps give discrete information. For local analysis, it is useful to consider derivatives of closed-loop poles with respect to parameters. Let \( p(s,\boldsymbol{\eta}) \) be the characteristic polynomial with parameter vector \( \boldsymbol{\eta} \in \mathbb{R}^{m} \):

\[ p(s,\boldsymbol{\eta}) = a_n(\boldsymbol{\eta}) s^n + \cdots + a_1(\boldsymbol{\eta}) s + a_0(\boldsymbol{\eta}). \]

For each root \( s_i(\boldsymbol{\eta}) \) satisfying \( p(s_i(\boldsymbol{\eta}),\boldsymbol{\eta}) = 0 \), differentiating implicitly with respect to parameter \( \eta_j \) yields

\[ \frac{\partial p}{\partial s}(s_i(\boldsymbol{\eta}),\boldsymbol{\eta}) \frac{\partial s_i}{\partial \eta_j} + \frac{\partial p}{\partial \eta_j}(s_i(\boldsymbol{\eta}),\boldsymbol{\eta}) = 0. \]

Assuming \( \partial p / \partial s \neq 0 \) at the root, the pole sensitivity is

\[ \frac{\partial s_i}{\partial \eta_j} = - \frac{ \dfrac{\partial p}{\partial \eta_j}(s_i(\boldsymbol{\eta}),\boldsymbol{\eta}) }{ \dfrac{\partial p}{\partial s}(s_i(\boldsymbol{\eta}),\boldsymbol{\eta}) }. \]

In practice, this derivative is often approximated numerically via finite differences in software. However, the analytic expression is useful to reason about how sensitive stability is to a particular parameter and to guide the choice of sweep ranges and resolutions.

Python, together with the python-control library, NumPy, and Matplotlib, is well suited for parameter sweeps. For robotics, roboticstoolbox-python can provide realistic inertia and damping parameters from robot models, which can then be fed into the control analysis.

import numpy as np
import matplotlib.pyplot as plt
import control  # python-control: pip install control
# Optional robotics toolbox (Peter Corke)
# from roboticstoolbox import DHRobot, RevoluteDH

# Parameter ranges for robot joint PD control
J_vals = np.linspace(0.5, 2.0, 20)      # inertia range
B_vals = np.linspace(0.2, 1.0, 20)      # viscous friction range
Kp_vals = np.linspace(5.0, 40.0, 10)    # proportional gain range
Kd_vals = np.linspace(0.1, 5.0, 10)     # derivative gain range

zeta_min = 0.5
Mp_max = 0.15
ts_max = 2.0

def closed_loop_metrics(J, B, Kp, Kd):
    # Plant G(s) = 1 / (J s^2 + B s)
    numG = [1.0]
    denG = [J, B, 0.0]
    G = control.TransferFunction(numG, denG)

    # PD controller C(s) = Kp + Kd s
    numC = [Kd, Kp]
    denC = [1.0]
    C = control.TransferFunction(numC, denC)

    T = control.feedback(C*G, 1)  # unity feedback

    # Step response and performance metrics
    t, y = control.step_response(T)
    info = control.step_info(T)

    Mp = info["Overshoot"] / 100.0  # convert percent to fraction
    ts = info["SettlingTime"]

    # Approximate damping by dominant poles
    poles = control.pole(T)
    # take pole with largest real part (closest to imaginary axis)
    p_dom = poles[np.argmax(np.real(poles))]
    wn = abs(p_dom)
    zeta = -np.real(p_dom) / wn if wn > 0 else 1.0

    return zeta, Mp, ts, poles

stable_count = 0
robust_count = 0
total_count = 0

for J in J_vals:
    for B in B_vals:
        for Kp in Kp_vals:
            for Kd in Kd_vals:
                total_count += 1
                zeta, Mp, ts, poles = closed_loop_metrics(J, B, Kp, Kd)

                if np.all(np.real(poles) < 0.0):
                    stable_count += 1
                    if (zeta >= zeta_min) and (Mp <= Mp_max) and (ts <= ts_max):
                        robust_count += 1

print("Stable fraction:", stable_count / total_count)
print("Robust-performance fraction:", robust_count / total_count)

# Example visualization: map of spec satisfaction over (Kp, Kd) for nominal (J,B)
J_nom, B_nom = 1.0, 0.5
sat_map = np.zeros((len(Kp_vals), len(Kd_vals)), dtype=int)

for i, Kp in enumerate(Kp_vals):
    for j, Kd in enumerate(Kd_vals):
        zeta, Mp, ts, poles = closed_loop_metrics(J_nom, B_nom, Kp, Kd)
        if np.all(np.real(poles) < 0.0) and (zeta >= zeta_min) and (Mp <= Mp_max) and (ts <= ts_max):
            sat_map[i, j] = 1

Kd_grid, Kp_grid = np.meshgrid(Kd_vals, Kp_vals)
plt.figure()
plt.contourf(Kd_grid, Kp_grid, sat_map, levels=[-0.5, 0.5, 1.5])
plt.xlabel("Kd")
plt.ylabel("Kp")
plt.title("Spec-satisfying region (1 = ok, 0 = fail)")
plt.colorbar()
plt.show()
      

In a robotics setting, one may loop over configurations of a manipulator using roboticstoolbox-python to obtain configuration-dependent joint inertia J(q), then embed that inside the sweep to obtain configuration-robust PD gains.

In C++, parameter sweeps can be implemented efficiently using nested loops. Linear algebra can be handled with Eigen, and integration with robotics frameworks such as ROS or ROS 2 allows these sweeps to run as off-line analysis nodes that read robot parameters from URDF models or dynamic parameter servers.

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

// Example: analytic metrics for second-order PD-controlled joint
struct Metrics {
    double zeta;
    double wn;
    bool stable;
};

Metrics metrics_from_params(double J, double B, double Kp, double Kd) {
    // Stability conditions for J s^2 + (B + Kd) s + Kp
    bool stable = (J > 0.0) && (B + Kd > 0.0) && (Kp > 0.0);

    double wn = std::sqrt(Kp / J);
    double zeta = (B + Kd) / (2.0 * std::sqrt(J * Kp));

    return {zeta, wn, stable};
}

int main() {
    std::vector<double> J_vals  = {0.8, 1.0, 1.2};
    std::vector<double> B_vals  = {0.3, 0.5, 0.7};
    std::vector<double> Kp_vals = {10.0, 20.0, 30.0};
    std::vector<double> Kd_vals = {0.5, 1.0, 2.0};

    double zeta_min = 0.5;
    double Mp_max   = 0.15;  // overshoot limit (fraction)
    double ts_max   = 2.0;   // approximate 2% settling time

    int total_count  = 0;
    int stable_count = 0;
    int robust_count = 0;

    for (double J : J_vals) {
        for (double B : B_vals) {
            for (double Kp : Kp_vals) {
                for (double Kd : Kd_vals) {
                    ++total_count;
                    Metrics m = metrics_from_params(J, B, Kp, Kd);

                    if (!m.stable) continue;
                    ++stable_count;

                    // Time-domain specs using standard second-order formulas
                    if (m.zeta <= 0.0 || m.zeta >= 1.0) continue;
                    double Mp = std::exp(-M_PI * m.zeta / std::sqrt(1.0 - m.zeta * m.zeta));
                    double ts = 4.0 / (m.zeta * m.wn);

                    if (m.zeta >= zeta_min && Mp <= Mp_max && ts <= ts_max) {
                        ++robust_count;
                    }
                }
            }
        }
    }

    std::cout << "Stable fraction: "
              << static_cast<double>(stable_count) / total_count << std::endl;
    std::cout << "Robust fraction: "
              << static_cast<double>(robust_count) / total_count << std::endl;

    return 0;
}
      

In a ROS-based robotics application, the nested loops could be wrapped inside a node that retrieves J and B from a dynamics library (e.g., using KDL or pinocchio), while logging safe regions of Kp and Kd for later use in controller tuning.

Java-based control and robotics systems (e.g., industrial control frameworks, educational robotics via WPILib) can use simple numerical utilities or libraries such as Apache Commons Math to implement parameter sweeps. Below is a self-contained parametric study of the PD-controlled joint using analytic second-order formulas.

public class PDSweep {

    static class Metrics {
        final double zeta;
        final double wn;
        final boolean stable;
        Metrics(double zeta, double wn, boolean stable) {
            this.zeta = zeta;
            this.wn = wn;
            this.stable = stable;
        }
    }

    static Metrics metricsFromParams(double J, double B, double Kp, double Kd) {
        boolean stable = (J > 0.0) && (B + Kd > 0.0) && (Kp > 0.0);
        double wn = Math.sqrt(Kp / J);
        double zeta = (B + Kd) / (2.0 * Math.sqrt(J * Kp));
        return new Metrics(zeta, wn, stable);
    }

    static double overshoot(double zeta) {
        if (zeta <= 0.0 || zeta >= 1.0) return 0.0;
        return Math.exp(-Math.PI * zeta / Math.sqrt(1.0 - zeta * zeta));
    }

    static double settlingTime(double zeta, double wn) {
        return 4.0 / (zeta * wn);
    }

    public static void main(String[] args) {
        double[] Jvals  = {0.8, 1.0, 1.2};
        double[] Bvals  = {0.3, 0.5, 0.7};
        double[] Kpvals = {10.0, 20.0, 30.0, 40.0};
        double[] Kdvals = {0.5, 1.0, 1.5, 2.0};

        double zetaMin = 0.5;
        double MpMax   = 0.15;
        double tsMax   = 2.0;

        int total   = 0;
        int stableN = 0;
        int robustN = 0;

        for (double J : Jvals) {
            for (double B : Bvals) {
                for (double Kp : Kpvals) {
                    for (double Kd : Kdvals) {
                        total++;
                        Metrics m = metricsFromParams(J, B, Kp, Kd);
                        if (!m.stable) continue;
                        stableN++;

                        double Mp = overshoot(m.zeta);
                        double ts = settlingTime(m.zeta, m.wn);

                        if (m.zeta >= zetaMin && Mp <= MpMax && ts <= tsMax) {
                            robustN++;
                        }
                    }
                }
            }
        }

        System.out.println("Stable fraction: " + (double) stableN / total);
        System.out.println("Robust fraction: " + (double) robustN / total);
    }
}
      

In Java-based robotics (e.g., differential drive or arm control in WPILib), the same parameter sweep logic can be adapted to test candidate gains for motor controllers against ranges of estimated inertias and frictions derived from identification experiments.

MATLAB with Control System Toolbox offers high-level functions for closed-loop construction and performance analysis. Simulink enables graphical models where controller and plant blocks are parameterized and swept using scripts or tools such as Simulink Design Optimization. A typical script-level sweep for the PD-controlled robot joint is:

% Parameter ranges
J_vals  = linspace(0.5, 2.0, 10);
B_vals  = linspace(0.2, 1.0, 10);
Kp_vals = linspace(5, 40, 8);
Kd_vals = linspace(0.1, 5, 8);

zeta_min = 0.5;
Mp_max   = 0.15;
ts_max   = 2.0;

stable_count = 0;
robust_count = 0;
total_count  = 0;

for J = J_vals
    for B = B_vals
        % Plant
        G = tf(1, [J B 0]);
        for Kp = Kp_vals
            for Kd = Kd_vals
                total_count = total_count + 1;

                % PD controller
                C = tf([Kd Kp], 1);

                T = feedback(C*G, 1);

                p = pole(T);
                if any(real(p) >= 0)
                    continue; % unstable
                end
                stable_count = stable_count + 1;

                S = stepinfo(T);
                Mp = S.Overshoot / 100;    % fraction
                ts = S.SettlingTime;

                % Approximate zeta from dominant pole
                [~, idx] = max(real(p));
                p_dom = p(idx);
                wn = abs(p_dom);
                zeta = -real(p_dom) / wn;

                if zeta >= zeta_min && Mp <= Mp_max && ts <= ts_max
                    robust_count = robust_count + 1;
                end
            end
        end
    end
end

fprintf("Stable fraction: %g\n", stable_count / total_count);
fprintf("Robust fraction: %g\n", robust_count / total_count);

% Visualization for nominal (J,B) using mesh plot
J_nom = 1.0; B_nom = 0.5;
G_nom = tf(1, [J_nom B_nom 0]);

[Kg, Kd_grid] = meshgrid(Kp_vals, Kd_vals);
sat = zeros(size(Kg));

for i = 1:numel(Kg)
    Kp = Kg(i);
    Kd = Kd_grid(i);
    C = tf([Kd Kp], 1);
    T = feedback(C*G_nom, 1);
    p = pole(T);
    if any(real(p) >= 0), continue; end
    S = stepinfo(T);
    Mp = S.Overshoot / 100;
    ts = S.SettlingTime;
    [~, idx] = max(real(p));
    p_dom = p(idx);
    wn = abs(p_dom);
    zeta = -real(p_dom) / wn;
    if zeta >= zeta_min && Mp <= Mp_max && ts <= ts_max
        sat(i) = 1;
    end
end

figure;
surf(Kp_vals, Kd_vals, sat');
xlabel("Kp"); ylabel("Kd"); zlabel("Spec OK (1/0)");
title("Spec-satisfying region at nominal J,B");
      

In Simulink, one can create a model with mask parameters J, B, Kp, and Kd, and then use a script that calls sim in nested loops to generate robustness maps, or use tools like parameter sweeps in Simulink Design Optimization to automate this process.

Wolfram Mathematica has symbolic and numeric capabilities for control. Parameter sweeps can be expressed using Table, and performance metrics can be extracted from closed-loop models defined via TransferFunctionModel and FeedbackConnect.

(* Parameter ranges *)
Jvals  = Range[0.5, 2.0, 0.3];
Bvals  = Range[0.2, 1.0, 0.2];
Kpvals = Range[5.0, 40.0, 5.0];
Kdvals = Range[0.1, 5.0, 0.5];

zetaMin = 0.5;
MpMax   = 0.15;
tsMax   = 2.0;

clearMetrics[J_, B_, Kp_, Kd_] := Module[
  {G, C, T, poles, stable, dom, wn, zeta, stepData, Mp, ts},
  G = TransferFunctionModel[{1}, {J s^2 + B s}, s];
  C = TransferFunctionModel[{Kd s + Kp}, {1}, s];
  T = FeedbackConnect[C*G, 1];

  poles = Eigenvalues[StateSpaceModel[T][[2]]]; (* A-matrix eigenvalues *)
  stable = Max[Re[poles]] < 0;

  dom  = First@SortBy[poles, Re]; (* largest real part *)
  wn   = Abs[dom];
  zeta = -Re[dom]/wn;

  stepData = StepResponse[T];
  Mp = (Max[stepData[[2]]] - 1);
  ts = Last@Select[stepData[[1]], 
    (Max @ Select[stepData[[2]], # > 0.98 && # < 1.02 &) &] &];

  << returning zeta, Mp, ts, stable;
  {zeta, Mp, ts, stable}
  ];

samples = Flatten[
   Table[
    {J, B, Kp, Kd, clearMetrics[J, B, Kp, Kd]},
    {J, Jvals}, {B, Bvals}, {Kp, Kpvals}, {Kd, Kdvals}
   ],
   3
];

stableCount = Count[samples[[All, 5, 4]], True];
totalCount  = Length[samples];

robustCount = Count[
   samples,
   {___, {z_, mp_, ts_, True}} /; 
     (z >= zetaMin && mp <= MpMax && ts <= tsMax)
];

N[stableCount/totalCount]
N[robustCount/totalCount]
      

For symbolic robustness, Mathematica can also manipulate parameter-dependent polynomials directly, enabling analytic checks of conditions such as \( J > 0 \), \( B + K_d > 0 \), and \( K_p > 0 \), before numerical sweeps are performed.

Problem 1 (Closed-loop parameters and specs): For the PD-controlled robot joint with \( G(s) = 1/(J s^2 + B s) \) and \( C(s) = K_p + K_d s \), derive expressions for \( \omega_n \) and \( \zeta \) in terms of \( J, B, K_p, K_d \). For fixed \( J \) and \( B \), show how to construct level sets in the \( (K_p,K_d) \) plane corresponding to constant \( \zeta \).

Solution: As derived earlier, the characteristic polynomial is

\[ J s^2 + (B + K_d) s + K_p = 0. \]

Comparing with \( s^2 + 2\zeta\omega_n s + \omega_n^2 = 0 \), we obtain

\[ \omega_n = \sqrt{\frac{K_p}{J}}, \qquad \zeta = \frac{B + K_d}{2 \sqrt{J K_p}}. \]

For fixed \( J \) and \( B \), the level set \( \zeta = \zeta_0 \) in the \( (K_p,K_d) \) plane satisfies

\[ \frac{B + K_d}{2 \sqrt{J K_p}} = \zeta_0 \quad \Rightarrow \quad K_d = 2 \zeta_0 \sqrt{J K_p} - B. \]

Thus, for each chosen \( \zeta_0 \), the relation between \( K_p \) and \( K_d \) is a curve of the form \( K_d(K_p) = 2 \zeta_0 \sqrt{J K_p} - B \), which can be plotted in software to guide selection of PD gains.

Problem 2 (Parameter sweep design for specifications): Using the same system, suppose the design specifications are \( M_p \le 0.1 \) and \( t_s \le 1.5 \) s. Assume \( J = 1 \), \( B = 0.5 \). For the three candidate gain pairs: \( (K_p,K_d) = (10,1), (20,2), (40,3) \), determine which meet the specifications using the standard formulas for \( M_p \) and \( t_s \).

Solution: First compute \( \omega_n = \sqrt{K_p} \) and \( \zeta = (B + K_d)/(2 \sqrt{K_p}) \). For each candidate:

• Case 1: \( K_p = 10, K_d = 1 \). Then \( \omega_n = \sqrt{10} \approx 3.16 \), \( \zeta = (0.5 + 1)/(2 \sqrt{10}) \approx 1.5/(6.32) \approx 0.237 \). Overshoot: \( M_p \approx \exp(-\pi \zeta/\sqrt{1-\zeta^2}) \approx 0.44 \), which violates \( M_p \le 0.1 \). This pair fails.
• Case 2: \( K_p = 20, K_d = 2 \). Then \( \omega_n = \sqrt{20} \approx 4.47 \), \( \zeta = (0.5 + 2)/(2 \sqrt{20}) = 2.5/(2 \cdot 4.47) \approx 0.28 \). Overshoot is still large (\( M_p \approx 0.35 \)). Settling time \( t_s \approx 4 / (\zeta \omega_n) \approx 4/(0.28 \cdot 4.47) \approx 3.2 \) s. Both specs are violated.
• Case 3: \( K_p = 40, K_d = 3 \). Then \( \omega_n = \sqrt{40} \approx 6.32 \), \( \zeta = (0.5 + 3)/(2 \sqrt{40}) = 3.5/(2 \cdot 6.32) \approx 0.277 \). Overshoot \( M_p \approx 0.36 \), and \( t_s \approx 4/(0.277 \cdot 6.32) \approx 2.3 \) s. Again, both specs fail.

None of the three gain pairs satisfy the required overshoot and settling-time specifications. A parameter sweep in software over a finer grid in \( (K_p,K_d) \) would identify gains with higher damping ratios and shorter settling times.

Problem 3 (Pole sensitivity formula derivation): Let \( p(s,\eta) \) be a scalar polynomial depending smoothly on a scalar parameter \( \eta \), and let \( s(\eta) \) be a simple root such that \( p(s(\eta),\eta) = 0 \). Derive the expression

\[ \frac{d s}{d\eta} = - \frac{\dfrac{\partial p}{\partial \eta}(s(\eta),\eta)} {\dfrac{\partial p}{\partial s}(s(\eta),\eta)}. \]

Solution: Differentiating the identity \( p(s(\eta),\eta) = 0 \) with respect to \( \eta \) gives, by the chain rule,

\[ \frac{d}{d\eta} p(s(\eta),\eta) = \frac{\partial p}{\partial s}(s(\eta),\eta)\,\frac{d s}{d\eta} + \frac{\partial p}{\partial \eta}(s(\eta),\eta) = 0. \]

Solving for \( ds/d\eta \) yields the stated formula, provided \( \partial p / \partial s \neq 0 \) at the root so that the implicit -function theorem applies.

Problem 4 (Grid resolution in parameter sweeps): Consider a two-parameter sweep over \( (K_p,K_d) \) on a rectangle \( K_p \in [K_p^{\min},K_p^{\max}] \), \( K_d \in [K_d^{\min},K_d^{\max}] \). Discuss qualitatively how the choice of grid resolution affects:

• the probability of missing small pockets of instability or poor performance, and
• the computational cost of the sweep.

Solution: A coarse grid (few points in each parameter direction) yields low computational cost but may miss narrow regions where poles cross the imaginary axis or where performance metrics sharply deteriorate. In particular, if the instability boundary is highly curved, coarse sampling can misclassify the system as robustly stable. Increasing grid resolution, i.e. using smaller step sizes in \( K_p \) and \( K_d \), reduces the risk of missing small regions but increases the number of samples approximately quadratically with resolution. In practice, one often uses multi-resolution strategies: coarse sweeps to locate interesting regions, followed by refined sweeps around boundaries of stability or performance satisfaction.

Problem 5 (Algorithmic structure of a sweep): Sketch a high-level algorithm (in words or pseudo-code) that, given ranges for \( J, B, K_p, K_d \), computes the fraction of parameter combinations that satisfy stability and a given set of performance inequalities.

Solution: A typical algorithm is:

flowchart TD
  S["Start: define ranges for J, B, Kp, Kd"] --> L1["Initialize counters: total, stable, robust"]
  L1 --> L2["For each J in grid"]
  L2 --> L3["For each B in grid"]
  L3 --> L4["For each Kp in grid"]
  L4 --> L5["For each Kd in grid"]
  L5 --> M["Build closed-loop model and compute metrics"]
  M --> C1["Check stability, update stable counter"]
  C1 --> C2["Check performance inequalities, update robust counter"]
  C2 --> N["After loops, compute ratios stable/total, robust/total"]
        

This structure is common across Python, C++, Java, MATLAB, and Mathematica implementations; only the inner step of “build closed-loop model and compute metrics” is tool-specific.

This lesson formalized parameter sweeps as finite-sample explorations of parameterized closed-loop systems, focusing on PD control of a robot joint as an illustrative example. We expressed stability and performance requirements in terms of closed-loop poles, damping ratio, natural frequency, overshoot, and settling time, and showed how these quantities depend on plant and controller parameters. We derived a local sensitivity formula for poles, then implemented global parameter sweeps in Python, C++, Java, MATLAB/Simulink, and Mathematica, with comments on integration with robotics toolchains. These techniques provide a practical bridge between theoretical robustness concepts (e.g., margins and sensitivity functions) and concrete controller tuning workflows used in engineering practice.