Consider a standard unity-feedback loop with controller \( C(s) \) and plant \( G(s) \). The open-loop transfer function is \( L(s) = C(s)G(s) \). In the frequency domain, we evaluate \( L(j\omega) \) for real \( \omega \) to obtain the sinusoidal steady-state response.

The Nichols plot represents the frequency response by plotting:

• Phase angle \( \varphi(\omega) \) (in degrees) on the horizontal axis.
• Magnitude in decibels \( m_{\mathrm{dB}}(\omega) = 20\log_{10}\left|L(j\omega)\right| \) on the vertical axis.

As the frequency \( \omega \) increases, the point \( (\varphi(\omega), m_{\mathrm{dB}}(\omega)) \) traces a curve in the phase–gain plane. This curve is the Nichols plot. It captures, in a single graph, the gain and phase behavior of the open-loop system across frequency.

Conceptually, constructing a Nichols plot from a mathematical model follows this workflow:

flowchart TD
  P["LTI plant + controller model"] --> LTF["Form open-loop L(s) = C(s)G(s)"]
  LTF --> FREQ["Evaluate L(jw) on frequency grid"]
  FREQ --> MP["Compute magnitude (dB) and phase (deg)"]
  MP --> NCH["Plot (phase, magnitude) curve = Nichols plot"]
        

In robotics and mechatronics, the plant \( G(s) \) may represent a motor-driven joint, a linear axis, or an actuated link. Nichols plots allow the designer to assess closed-loop behavior (studied in later lessons) using only the open-loop model.

Let \( L(s) \) be a proper rational transfer function. For sinusoidal steady-state analysis we restrict to the imaginary axis \( s = j\omega \), \( \omega \in \mathbb{R} \), and define:

\[ L(j\omega) = \Re\!\big(L(j\omega)\big) + j\,\Im\!\big(L(j\omega)\big). \]

The magnitude and phase of the complex number \( L(j\omega) \) are

\[ \left|L(j\omega)\right| = \sqrt{\Re\!\big(L(j\omega)\big)^2 + \Im\!\big(L(j\omega)\big)^2}, \quad \varphi(\omega) = \operatorname{atan2}\!\Big(\Im\!\big(L(j\omega)\big),\, \Re\!\big(L(j\omega)\big)\Big), \]

where \( \operatorname{atan2} \) returns a phase in radians, typically in \( (-\pi,\pi] \). In control engineering we work in degrees, so we define

\[ \varphi_{\deg}(\omega) = \frac{180}{\pi}\,\varphi(\omega). \]

The magnitude in decibels is defined as

\[ m_{\mathrm{dB}}(\omega) = 20\log_{10}\left|L(j\omega)\right|. \]

The Nichols mapping \( \mathcal{N} \) assigns to each frequency \( \omega \) the point

\[ \mathcal{N}(\omega) = \Big(\,\varphi_{\deg}(\omega),\; m_{\mathrm{dB}}(\omega)\Big). \]

The Nichols curve of the open loop \( L(s) \) is the parametric curve \( \omega \mapsto \mathcal{N}(\omega) \), typically drawn for \( \omega > 0 \) over a wide range of frequencies relevant to the dynamics of interest (e.g., below and above the mechanical resonance of a robot joint).

The Nyquist plot represents the mapping \( \omega \mapsto L(j\omega) \in \mathbb{C} \) directly in the complex plane. Write \( L(j\omega) = x(\omega) + j y(\omega) \). Then the Nyquist curve is the locus of \( (x(\omega), y(\omega)) \).

The Nichols plot is obtained from the same frequency response by the transformation

\[ (x(\omega),y(\omega)) \mapsto \Big( \varphi_{\deg}(\omega), 20\log_{10}\sqrt{x(\omega)^2 + y(\omega)^2} \Big). \]

Thus, Nichols plots and Nyquist plots contain the same information but in different coordinates. The Nichols representation is particularly convenient for overlaying contours of closed-loop quantities (discussed in the next lesson).

A Bode plot consists of two separate graphs:

• Magnitude \( 20\log_{10}\left|L(j\omega)\right| \) vs. \( \log_{10}\omega \).
• Phase \( \varphi_{\deg}(\omega) \) vs. \( \log_{10}\omega \).

Given a Bode magnitude curve \( M_{\mathrm{dB}}(\omega) \) and phase curve \( \Phi(\omega) \), the corresponding Nichols curve is obtained by eliminating the explicit dependence on \( \omega \):

\[ \mathcal{N}(\omega) = \big(\Phi(\omega),\, M_{\mathrm{dB}}(\omega)\big). \]

That is, Nichols plots can be seen as a recombination of the Bode magnitude and phase plots into a single two-dimensional curve in the phase–gain plane.

Consider a plant with one integrator and one real pole, controlled by a proportional gain \( K > 0 \):

\[ L(s) = \frac{K}{s(1 + sT)}, \quad T > 0. \]

Substituting \( s = j\omega \) yields

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

Using magnitude and phase of products and reciprocals, we obtain

\[ \left|L(j\omega)\right| = \frac{K}{\left|j\omega\right|\;\left|1 + j\omega T\right|} = \frac{K}{\omega\,\sqrt{1 + (\omega T)^2}}, \]

\[ \varphi_{\deg}(\omega) = \arg K - 90^{\circ} - \arctan(\omega T) = -90^{\circ} - \arctan(\omega T), \]

assuming \( K > 0 \), so that \( \arg K = 0^{\circ} \). The magnitude in decibels is

\[ m_{\mathrm{dB}}(\omega) = 20\log_{10} K - 20\log_{10}\omega - 10\log_{10}\big(1 + (\omega T)^2\big). \]

From these expressions we can deduce asymptotic behavior:

• As \( \omega \to 0^{+} \): \( \left|L(j\omega)\right| \approx K/\omega \) (slope \( -20 \) dB/decade), and \( \varphi_{\deg}(\omega) \approx -90^{\circ} \).
• As \( \omega \to \infty \): \( \left|L(j\omega)\right| \approx K/(\omega^2 T) \) (slope \( -40 \) dB/decade), and \( \varphi_{\deg}(\omega) \approx -180^{\circ} \).

Plotting \( \big(\varphi_{\deg}(\omega), m_{\mathrm{dB}}(\omega)\big) \) as \( \omega \) increases from very small to very large values yields a Nichols curve that starts near \( (\,-90^{\circ}, +\infty\,) \) and ends near \( (\,-180^{\circ}, -\infty\,) \), bending smoothly between these asymptotes.

Many robotic actuators (e.g., a motor plus load inertia with viscous friction) can be approximated by models of this form near a nominal operating point, so understanding this Nichols shape is directly relevant to servo design.

In practice, we compute Nichols plots numerically on a finite grid of frequencies. A typical algorithm for a SISO LTI system with transfer function \( L(s) \) proceeds as follows:

1. Choose minimum and maximum frequencies \( \omega_{\min} \), \( \omega_{\max} \) that cover the dynamic range of interest.
2. Construct a logarithmically spaced grid \( \{\omega_k\}_{k=0}^{N-1} \):

   \[ \omega_k = \omega_{\min} \left(\frac{\omega_{\max}}{\omega_{\min}}\right)^{\frac{k}{N-1}}, \quad k = 0,\dots,N-1. \]

3. For each \( \omega_k \), form \( s_k = j\omega_k \) and compute \( L(j\omega_k) \).
4. Compute magnitude and phase: \( m_{\mathrm{dB}}(\omega_k) \) and \( \varphi_{\deg}(\omega_k) \).
5. Phase unwrapping: adjust the sequence \( \varphi_{\deg}(\omega_k) \) to avoid artificial jumps of around \( \pm 360^{\circ} \) introduced by the principal value of the argument. This ensures a visually continuous curve.
6. Plot the discrete set of points \( (\varphi_{\deg}(\omega_k), m_{\mathrm{dB}}(\omega_k)) \) and join them with line segments.

The overall computational pipeline for numerical tools can be summarized as:

flowchart TD
  A["Specify L(s) (tf or state-space)"] --> B["Select [w_min, w_max] and N"]
  B --> C["Compute L(jw_k) for grid"]
  C --> D["Get mag_dB(w_k), phase_deg(w_k)"]
  D --> E["Unwrap phase sequence"]
  E --> F["Plot phase_deg vs mag_dB in phase–gain plane"]
        

Numerical robustness requires avoiding frequencies where \( L(j\omega) \) is extremely close to zero (to avoid underflow) or extremely large (overflow). Well-chosen frequency ranges and logarithmic spacing generally prevent such issues for typical control plants.

In Python, the python-control library is widely used for education and research in control and robotics. Consider a simplified DC motor model for a single robot joint:

• Inertia \( J \), viscous friction \( b \).
• Motor torque constant \( K_t \).
• P controller with gain \( K_p \).

The transfer function from motor command to joint angle (linearized) can often be approximated as \( G(s) = \dfrac{K_t}{Js^2 + bs} \), so the open-loop transfer function is \( L(s) = K_p G(s) \).

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

# DC motor / joint parameters (simple illustrative values)
J = 0.01   # kg*m^2
b = 0.1    # N*m*s/rad
Kt = 0.01  # N*m/A
Kp = 5.0   # proportional controller gain

# Transfer functions: G(s) = Kt / (J s^2 + b s), C(s) = Kp
s = ct.TransferFunction.s
G = Kt / (J * s**2 + b * s)
C = Kp
L = C * G  # open-loop transfer function

# Frequency grid (rad/s)
w = np.logspace(-1, 3, 400)

# Frequency response
mag, phase, omega = ct.freqresp(L, w)  # mag, phase are complex arrays

# Convert to Nichols coordinates
mag = np.abs(mag.squeeze())
phase_rad = np.angle(phase.squeeze())  # some versions return separate phase
mag_db = 20.0 * np.log10(mag)
phase_deg = 180.0 / np.pi * phase_rad

# Optional phase unwrapping for visual continuity
phase_rad_unwrapped = np.unwrap(phase_rad)
phase_deg_unwrapped = 180.0 / np.pi * phase_rad_unwrapped

# Plot Nichols curve
plt.figure()
plt.plot(phase_deg_unwrapped, mag_db)
plt.xlabel("Phase (deg)")
plt.ylabel("Magnitude (dB)")
plt.title("Nichols Plot of Robotic Joint Open Loop")
plt.grid(True)
plt.gca().invert_xaxis()  # conventional: phase decreases to the right
plt.show()
      

In a robotics workflow (e.g., using ROS), such a Nichols plot can be used to tune \( K_p \) (and later more sophisticated controllers) so that the closed-loop tracking and robustness requirements are met. Higher-level libraries for robot modeling (such as pinocchio or robotics-toolbox-python) can provide state-space models that are converted to transfer functions before computing Nichols plots.

In embedded robotic controllers (e.g., running under ROS 2 or a real-time OS), C++ is often used. We can compute Nichols data on the host to analyze controllers before deployment. The following example computes magnitude and phase for \( L(s) = \dfrac{K}{s(1 + sT)} \) on a logarithmic grid:

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

using cd = std::complex<double>;

int main() {
    const double K = 10.0;
    const double T = 0.1;

    const double w_min = 0.1;
    const double w_max = 1000.0;
    const std::size_t N = 200;

    std::vector<double> w(N);
    std::vector<double> mag_dB(N);
    std::vector<double> phase_deg(N);

    const double log_w_min = std::log10(w_min);
    const double log_w_max = std::log10(w_max);

    for (std::size_t k = 0; k < N; ++k) {
        double alpha = static_cast<double>(k) / static_cast<double>(N - 1);
        double log_w = log_w_min + alpha * (log_w_max - log_w_min);
        w[k] = std::pow(10.0, log_w);

        cd s(0.0, w[k]);              // s = j w
        cd L = K / (s * (1.0 + s * T)); // L(jw)

        double mag = std::abs(L);
        double phase_rad = std::arg(L);

        mag_dB[k] = 20.0 * std::log10(mag);
        phase_deg[k] = 180.0 / M_PI * phase_rad;
    }

    // Export data (phase_deg, mag_dB) for plotting in Python, gnuplot, etc.
    for (std::size_t k = 0; k < N; ++k) {
        std::cout << phase_deg[k] << " " << mag_dB[k] << "\n";
    }

    return 0;
}
      

In a robotics environment, the plant model \( L(s) \) can be obtained from an Eigen-based state-space model or from identification data. The resulting Nichols data can be plotted offline to tune gains for joint position or velocity controllers before being implemented inside the real-time control loop (for example in a ROS 2 controller using ros2_control).

Java appears in certain robotic stacks (for example, higher-level control or simulation frameworks). Below we implement a minimal complex-number utility and compute Nichols data for the same \( L(s) \) as in the C++ example:

public class NicholsBasic {

    // Minimal complex class for frequency response computations
    static class Complex {
        final double re;
        final double im;

        Complex(double re, double im) {
            this.re = re;
            this.im = im;
        }

        Complex add(Complex other) {
            return new Complex(this.re + other.re, this.im + other.im);
        }

        Complex mul(Complex other) {
            return new Complex(
                this.re * other.re - this.im * other.im,
                this.re * other.im + this.im * other.re
            );
        }

        Complex div(Complex other) {
            double denom = other.re * other.re + other.im * other.im;
            return new Complex(
                (this.re * other.re + this.im * other.im) / denom,
                (this.im * other.re - this.re * other.im) / denom
            );
        }

        double abs() {
            return Math.hypot(re, im);
        }

        double arg() {
            return Math.atan2(im, re);
        }
    }

    public static void main(String[] args) {
        double K = 10.0;
        double T = 0.1;

        double wMin = 0.1;
        double wMax = 1000.0;
        int N = 200;

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

            Complex s = new Complex(0.0, w); // s = j w
            Complex denom = s.mul(new Complex(1.0, 0.0).add(s.mul(new Complex(T, 0.0))));
            Complex L = new Complex(K, 0.0).div(denom);

            double mag = L.abs();
            double magDb = 20.0 * Math.log10(mag);
            double phaseDeg = 180.0 / Math.PI * L.arg();

            System.out.println(phaseDeg + " " + magDb);
        }
    }
}
      

The printed (phase, magnitude) pairs can be imported into plotting tools or robotics-oriented visualization frameworks to inspect the Nichols curve while designing motion controllers in Java-based environments.

MATLAB and Simulink are standard tools in control and robotics. The Control System Toolbox provides built-in functions for Nichols plots. Using the same DC motor joint model:

% Parameters
J  = 0.01;
b  = 0.1;
Kt = 0.01;
Kp = 5.0;

s = tf('s');
G = Kt / (J * s^2 + b * s);  % plant
C = Kp;                      % proportional controller
L = C * G;                   % open-loop transfer function

% Simple Nichols plot using Control System Toolbox
figure;
nichols(L);
grid on;
title('Nichols Plot of Robotic Joint Open Loop');

% Manual Nichols construction from Bode data
w = logspace(-1, 3, 400);
[mag, phase] = bode(L, w);   % mag, phase: 3D arrays

mag = squeeze(mag);
phase = squeeze(phase);

magDb = 20 * log10(mag);
phaseDeg = phase;            % already in degrees

figure;
plot(phaseDeg, magDb);
xlabel('Phase (deg)');
ylabel('Magnitude (dB)');
title('Nichols Plot (Constructed from Bode Data)');
grid on;
set(gca, 'XDir', 'reverse'); % conventional Nichols orientation
      

In Simulink, a robotics engineer might have a detailed multi-body model of a manipulator built with the Robotics System Toolbox (e.g., using a rigidBodyTree). By linearizing the model around a nominal operating point with the Linear Analysis Tool, one can directly generate Nichols plots for the open-loop transfer function between motor torques and joint positions and use them to tune classical controllers.

Wolfram Mathematica can also be used to construct Nichols plots symbolically or numerically. For the servo model \( L(s) = \dfrac{K_p K_t}{Js^2 + bs} \), we may write:

(* Parameters *)
J  = 0.01;
b  = 0.1;
Kt = 0.01;
Kp = 5.0;

(* Frequency response L(j w) *)
L[w_] := (Kp*Kt) / (J*(I*w)^2 + b*(I*w));

magDb[w_]  := 20*Log10[Abs[L[w]]];
phaseDeg[w_] := 180./Pi*Arg[L[w]];

wMin = 0.1;
wMax = 1000.0;

data = Table[
   {phaseDeg[w], magDb[w]},
   {w, wMin, wMax, (wMax/wMin)^(1./199) - 1.0 + wMin} (* approximate log spacing *)
];

ListLinePlot[
  data,
  AxesLabel -> {"Phase (deg)", "Magnitude (dB)"},
  PlotRange -> All,
  PlotTheme -> "Detailed",
  Prolog -> {AbsoluteThickness[1.5]},
  ScalingFunctions -> {None, None},
  ImageSize -> Large
]
      

Symbolic manipulation features also allow derivation of exact magnitude and phase expressions for simple transfer functions, which is useful when teaching the theoretical underpinnings of Nichols plots.

Problem 1 (Nichols Plot of a Pure Gain): Let \( L(s) = K \) with \( K > 0 \). Compute the Nichols plot and describe the curve.

Solution: For all \( \omega \), we have

\[ L(j\omega) = K, \quad \left|L(j\omega)\right| = K, \quad \varphi_{\deg}(\omega) = 0^{\circ}. \]

Thus \( m_{\mathrm{dB}}(\omega) = 20\log_{10} K \) is constant and the Nichols curve collapses to a single point \( (0^{\circ}, 20\log_{10} K) \). A positive pure gain does not introduce any phase shift; it only shifts the magnitude vertically in the Nichols plane.

Problem 2 (First-Order Factor): For \( L(s) = \dfrac{K}{1 + sT} \) with \( K > 0 \) and \( T > 0 \), derive expressions for magnitude and phase as functions of \( \omega \) and describe qualitatively the Nichols curve.

Solution: Substituting \( s = j\omega \):

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

The magnitude and phase are

\[ \left|L(j\omega)\right| = \frac{K}{\sqrt{1 + (\omega T)^2}}, \quad \varphi_{\deg}(\omega) = -\arctan(\omega T)\cdot\frac{180}{\pi}. \]

Hence

\[ m_{\mathrm{dB}}(\omega) = 20\log_{10} K - 10\log_{10}\big(1 + (\omega T)^2\big). \]

As \( \omega \to 0^{+} \), we have \( \left|L(j\omega)\right| \to K \) and \( \varphi_{\deg}(\omega) \to 0^{\circ} \). As \( \omega \to \infty \), we have \( \left|L(j\omega)\right| \to 0 \) and \( \varphi_{\deg}(\omega) \to -90^{\circ} \). The Nichols curve is therefore a smooth trajectory from \( (0^{\circ}, 20\log_{10} K) \) toward \( (-90^{\circ}, -\infty) \).

Problem 3 (First-Order Factor with Integrator): For \( L(s) = \dfrac{K}{s(1 + sT)} \), verify the formulas for magnitude and phase given in Section 4 and sketch the qualitative shape of the Nichols plot.

Solution: From Section 4 we already have

\[ \left|L(j\omega)\right| = \frac{K}{\omega\sqrt{1 + (\omega T)^2}}, \quad \varphi_{\deg}(\omega) = -90^{\circ} - \arctan(\omega T). \]

As \( \omega \to 0^{+} \), \( \left|L(j\omega)\right| \to \infty \) and \( \varphi_{\deg}(\omega) \to -90^{\circ} \); as \( \omega \to \infty \), \( \left|L(j\omega)\right| \to 0 \) and \( \varphi_{\deg}(\omega) \to -180^{\circ} \). Thus the Nichols curve begins near very high magnitude around \( -90^{\circ} \) and ends near very low magnitude around \( -180^{\circ} \). In between, it bends smoothly downward and rightward in the phase–gain plane.

Problem 4 (From Bode Data to Nichols Points): Suppose the Bode plot of an open-loop transfer function gives the following data:

• At \( \omega_1 = 1 \) rad/s: \( M_{\mathrm{dB}}(\omega_1) = -6 \) dB, \( \Phi(\omega_1) = -60^{\circ} \).
• At \( \omega_2 = 10 \) rad/s: \( M_{\mathrm{dB}}(\omega_2) = -20 \) dB, \( \Phi(\omega_2) = -120^{\circ} \).

What are the corresponding Nichols coordinates for these two frequencies?

Solution: By definition, the Nichols coordinates are \( (\Phi(\omega), M_{\mathrm{dB}}(\omega)) \). Hence:

• At \( \omega_1 \): Nichols point \( (-60^{\circ}, -6\ \mathrm{dB}) \).
• At \( \omega_2 \): Nichols point \( (-120^{\circ}, -20\ \mathrm{dB}) \).

These points lie on the Nichols curve; interpolating between many such points reconstructs the entire Nichols plot.

Problem 5 (Effect of Time Delay): Let \( L_0(s) \) be an open-loop transfer function with Nichols curve \( (\varphi_{\deg}(\omega), m_{\mathrm{dB}}(\omega)) \). Consider adding a pure time delay \( e^{-sT_d} \) to obtain \( L(s) = L_0(s) e^{-sT_d} \). Show that the Nichols curve of \( L(s) \) is a horizontal shift of the Nichols curve of \( L_0(s) \).

Solution: For \( s = j\omega \):

\[ L(j\omega) = L_0(j\omega) e^{-j\omega T_d}. \]

The magnitude of the delay term is \( \left|e^{-j\omega T_d}\right| = 1 \), so \( \left|L(j\omega)\right| = \left|L_0(j\omega)\right| \) and \( m_{\mathrm{dB}}(\omega) \) is unchanged. The phase of the delay term is \( \arg(e^{-j\omega T_d}) = -\omega T_d \) (radians), so

\[ \varphi_{\deg}(\omega) = \varphi_{0,\deg}(\omega) - \frac{180}{\pi}\,\omega T_d. \]

Therefore, for each frequency the Nichols point is shifted horizontally by \( -\dfrac{180}{\pi}\omega T_d \) degrees, while the magnitude in dB is unchanged. The Nichols curve of \( L(s) \) is a frequency-dependent horizontal shift of the Nichols curve of \( L_0(s) \), which tends to reduce phase margin in feedback design.

In this lesson we defined the Nichols plot as a parametric curve in the phase–gain plane, derived directly from the open-loop frequency response \( L(j\omega) \). We related Nichols plots to Nyquist and Bode representations, derived analytical expressions for simple first-order and integrator-based plants that approximate many robotic actuators, and formulated a practical numerical algorithm for constructing Nichols plots. Implementations in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica show how Nichols plots are generated in modern control and robotics workflows. In the next lesson we will use Nichols charts with overlaid closed-loop contours to interpret stability and performance directly in the Nichols plane.