Chapter 28: Computer-Aided Analysis and Design for Linear Control
Lesson 3: Root Locus, Bode, and Nyquist Design via Software
This lesson shows how classical SISO design methods based on root locus, Bode plots, and Nyquist diagrams are implemented and automated in standard software environments. We connect the underlying mathematics of the loop transfer function, closed-loop poles, and stability margins to concrete workflows in Python, MATLAB/Simulink, C++, Java, and Wolfram Mathematica, with remarks on their use in robotics-oriented libraries.
1. Mathematical Foundations for Software-Based Classical Design
Consider a unity-feedback single-input single-output (SISO) loop with loop transfer function \( L(s) = K G(s) \), where \( G(s) \) is the plant transfer function and \( K \) is a scalar gain (or more generally, a controller parameter). The closed-loop transfer function is
\[ T(s) = \frac{K G(s)}{1 + K G(s)} = \frac{L(s)}{1 + L(s)}. \]
The characteristic equation is \( 1 + L(s) = 0 \). All three classical design views are different ways of analyzing this equation:
- Root locus: For real nonnegative gains \( K \), study how the roots of \( 1 + K G(s) = 0 \) move in the complex plane as \( K \) varies.
- Bode plots: For \( s = j\omega \), study \( L(j\omega) \) in terms of magnitude and phase across frequency.
- Nyquist plot: Study the contour \( L(j\omega) \) as \( \omega \) sweeps from zero to high frequency, using encirclements of \( -1 \) to infer closed-loop stability.
The root locus is defined implicitly by the conditions
\[ 1 + K G(s) = 0 \quad \Longleftrightarrow \quad K G(s) = -1. \]
Writing \( G(s) = |G(s)| e^{j \angle G(s)} \), the root locus consists of points \( s \) satisfying:
\[ |K G(s)| = 1, \quad \angle G(s) = (2\ell + 1)\pi, \ \ell \in \mathbb{Z}. \]
Software packages implement root locus by computing the roots of the polynomial \( 1 + K G(s) \) (typically expressed as a ratio of polynomials with real coefficients) over a dense grid of \( K \) values and plotting the resulting pole locations.
Frequency-domain methods consider \( L(j\omega) \). Its magnitude and phase are
\[ |L(j\omega)| = \sqrt{\Re\{L(j\omega)\}^2 + \Im\{L(j\omega)\}^2}, \quad \angle L(j\omega) = \operatorname{atan2}\!\big(\Im\{L(j\omega)\}, \Re\{L(j\omega)\}\big). \]
The Bode magnitude in decibels (dB) is
\[ M_{\text{dB}}(\omega) = 20 \log_{10} |L(j\omega)|. \]
Software tools numerically evaluate these expressions over a frequency grid and plot magnitude and phase versus \( \log_{10} \omega \). The Nyquist diagram then traces the complex curve \( L(j\omega) \) and its mirror (for real-rational stable plants) in the complex plane.
flowchart TD
P["Plant model G(s)"] --> L["Loop L(s) = K G(s)"]
L --> RL["Solve 1 + L(s) = 0 for s & \nK (root locus)"]
L --> BO["Evaluate L(j*omega) \non grid (Bode)"]
L --> NY["Trace L(j*omega) in \ncomplex plane (Nyquist)"]
RL --> DES["Select K and controller structure"]
BO --> DES
NY --> DES
DES --> ROB["Implement controller in embedded/robotics software"]
2. Stability Margins and Nyquist Count in Software
The Nyquist stability criterion relates the open-loop transfer function \( L(s) \) to closed-loop pole locations through the mapping of a contour enclosing the right-half plane (RHP). Let \( P \) be the number of open-loop poles of \( L(s) \) in the RHP, and let \( Z \) be the number of closed-loop poles of \( T(s) \) in the RHP. Let \( N \) be the net number of clockwise encirclements of \( -1 \) by the Nyquist plot of \( L(s) \). Then
\[ N = Z - P. \]
For a loop with no open-loop RHP poles (\( P = 0 \)), closed-loop stability requires \( Z = 0 \), thus \( N = 0 \). Software tools detect encirclements of \( -1 \) automatically by tracking the phase of \( L(j\omega) \) as \( \omega \) sweeps a sufficiently dense frequency grid.
From the Bode plot of \( L(j\omega) \), the classical gain margin and phase margin are defined as:
\[ \omega_{\text{gc}} : M_{\text{dB}}(\omega_{\text{gc}}) = 0, \quad \varphi_m = 180^\circ + \angle L(j\omega_{\text{gc}}), \]
\[ \omega_{\text{pc}} : \angle L(j\omega_{\text{pc}}) = -180^\circ, \quad G_m = - M_{\text{dB}}(\omega_{\text{pc}}). \]
Here \( \omega_{\text{gc}} \) is the gain crossover frequency and \( \omega_{\text{pc}} \) the phase crossover frequency. Software packages compute these quantities by interpolating magnitude and phase on the sampled grid, then provide them as scalar outputs that can be used in automated design loops.
In robotics applications (e.g., joint position control of a manipulator link), these margins are used to ensure robust behavior under modeling uncertainty and varying loads. Software-generated margins can be used as constraints in gain tuning or compensator parameter optimization.
3. Software-Based Design Workflow
The typical computer-aided workflow for classical design combines root locus, Bode, and Nyquist plots:
- Specify the plant transfer function \( G(s) \) from a previously developed model (mechanical, electrical, or robotic joint dynamics).
- Encode performance specifications: overshoot and settling-time bounds in the time domain, and gain/phase margins or bandwidth in the frequency domain.
- Use root locus to select a feasible range of gains or to determine how placing additional poles/zeros (lead or lag compensators) will move closed-loop poles into the desired region.
- Cross-check with Bode and Nyquist plots to ensure adequate robustness margins and noise attenuation.
- Iterate the controller parameters until both time-domain and frequency-domain constraints are met.
flowchart TD
S["Start: plant G(s) & specs"] --> RL["Root locus: choose K, add poles/zeros"]
RL --> TD["Check time response \n(step, disturbance)"]
RL --> BD["Bode: compute margins, bandwidth"]
BD --> NY["Nyquist: verify encirclements and robustness"]
TD --> IT["Adjust controller params (K, lead/lag, filters)"]
NY --> IT
IT -->|repeat| RL
IT --> DONE["Freeze design & export to embedded/robotics code"]
In the following sections we illustrate this workflow in several software environments, using a representative SISO plant such as a DC motor or robotic joint approximated by \( G(s) = \dfrac{K_m}{s(J s + b)} \), where \( J \) is inertia, \( b \) viscous friction, and \( K_m \) the motor torque constant.
4. Python Implementation of Root Locus, Bode, and Nyquist
Python has mature libraries for linear control analysis. A common stack
is
numpy,
scipy, and
control (often called
python-control). In robotics, these
tools are frequently used together with ROS-based libraries or robotics
toolboxes that provide plant models for actuators and joints.
The example below defines a simple plant \( G(s) = \dfrac{K_m}{s(J s + b)} \), plots the root locus, computes Bode margins, and draws a Nyquist diagram.
import numpy as np
import control as ct
import matplotlib.pyplot as plt
# Physical parameters for a DC motor / robotic joint approximation
J = 0.01 # inertia
b = 0.1 # viscous friction
Km = 1.0 # motor torque constant
# Transfer function G(s) = Km / (s (J s + b))
s = ct.TransferFunction.s
G = Km / (s * (J * s + b))
# 1) Root locus plot and gain selection
plt.figure()
rlocus_data = ct.root_locus(G, kvect=np.linspace(0, 200, 400), plot=True)
# rlocus_data is (poles, gains); can be post-processed numerically
# Example: pick a gain K that gives acceptable damping (inspect plot)
K_des = 50.0
L = K_des * G # loop transfer function with proportional gain
T = ct.feedback(L, 1) # closed-loop transfer function
# 2) Bode plot and stability margins
plt.figure()
mag, phase, omega = ct.bode(L, dB=True, Hz=False, omega_limits=(0.1, 100), omega_num=500)
gm, pm, wgm, wpm = ct.margin(L)
print("Gain margin (dB):", 20 * np.log10(gm) if gm is not None else None)
print("Phase margin (deg):", pm)
print("Gain crossover rad/s:", wgm)
print("Phase crossover rad/s:", wpm)
# 3) Nyquist plot
plt.figure()
ct.nyquist_plot(L)
plt.show()
For robotic applications,
python-control can be combined with
robotics libraries (e.g., toolboxes that model manipulators) to
linearize joint dynamics about operating points, then apply exactly the
same root locus, Bode, and Nyquist calls to the resulting SISO axis
models.
5. MATLAB/Simulink Implementation
MATLAB provides native functions such as tf,
rlocus, bode, margin, and
nyquist. Simulink adds block-diagram simulation, which is
widely used in robotics and mechatronics. Robotics-oriented toolboxes
provide manipulator and mobile robot models that can be linearized for
axis-level control design.
% Parameters
J = 0.01;
b = 0.1;
Km = 1.0;
s = tf('s');
G = Km / (s * (J * s + b)); % Plant
% 1) Root locus
figure;
rlocus(G);
grid on;
title('Root locus of G(s)');
% Choose gain K from root locus (e.g., by clicking or inspection)
K_des = 50;
L = K_des * G; % Loop transfer function
T = feedback(L, 1); % Closed-loop transfer function
% 2) Bode and margins
figure;
margin(L); % Bode plot with gain/phase margins
grid on;
% 3) Nyquist plot
figure;
nyquist(L);
grid on;
% Step response to verify time-domain specs
figure;
step(T);
grid on;
title('Closed-loop step response');
In a robotics workflow, the same commands are applied to transfer
functions obtained from manipulator joint models. Simulink models of
robot actuators may include nonlinearities; linearization blocks can
extract the local linear model for which rlocus,
bode, and nyquist are used to design joint
controllers.
Simulink block diagrams also allow you to implement the final controller structure (e.g., proportional plus lead compensation), simulate with realistic disturbances and sensor noise, and then generate embedded code for deployment on robotic platforms.
6. C++ Implementation and Robotics Libraries
In C++, there is no single standard control library, but robotics and
embedded projects often rely on
Eigen for linear algebra together with
robotics frameworks such as ROS and various controller packages. Root
locus and frequency-response computations can be implemented using
polynomial solvers and complex arithmetic.
The following sketch evaluates the open-loop frequency response \( L(j\omega) \) for Bode and Nyquist-like plots, for a plant \( G(s) = \dfrac{K_m}{s(J s + b)} \).
#include <iostream>
#include <complex>
#include <vector>
int main() {
double J = 0.01;
double b = 0.1;
double Km = 1.0;
double K = 50.0; // loop gain
using cdouble = std::complex<double>;
// Frequency grid (rad/s)
std::vector<double> omega;
for (double w = 0.1; w <= 100.0; w *= 1.1) {
omega.push_back(w);
}
std::vector<cdouble> L_vals;
for (double w : omega) {
cdouble s(0.0, w); // j * w
cdouble G = Km / (s * (J * s + b));
cdouble L = K * G;
L_vals.push_back(L);
}
// Print magnitude (dB) and phase (deg) for Bode-style output
for (std::size_t i = 0; i < omega.size(); ++i) {
cdouble L = L_vals[i];
double mag = std::abs(L);
double phase = std::arg(L) * 180.0 / M_PI;
double mag_db = 20.0 * std::log10(mag);
std::cout << omega[i] << " " << mag_db
<< " " << phase << std::endl;
}
// The resulting data can be plotted using a plotting library or exported
// to Python/MATLAB. A Nyquist diagram uses the complex values L_vals directly.
return 0;
}
Root locus can be approximated by forming the characteristic polynomial coefficients of \( 1 + K G(s) \) symbolically in terms of \( K \), then sweeping over \( K \) and using polynomial root solvers to compute closed-loop poles for each value. In robotics codebases, similar numerical kernels are often integrated into controllers that have been tuned offline using higher-level tools but require on-board checks of margins or bandwidth.
7. Java Implementation and Robotics Context
Java projects can use numerical libraries such as Apache Commons Math for complex arithmetic and polynomial root solving. In educational and robotics competitions, Java-based robotics libraries (e.g., those used in mobile platforms or robotic arms) often embed linear control logic that can be validated offline with Bode and Nyquist computations.
The following example evaluates the open-loop frequency response for a Bode-style analysis:
import java.util.ArrayList;
import java.util.List;
public class BodeExample {
public static void main(String[] args) {
double J = 0.01;
double b = 0.1;
double Km = 1.0;
double K = 50.0;
List<Double> omega = new ArrayList<>();
for (double w = 0.1; w <= 100.0; w *= 1.1) {
omega.add(w);
}
for (double w : omega) {
// s = j * w
double real = 0.0;
double imag = w;
// G(s) = Km / (s (J s + b))
// Complex arithmetic by hand:
// s = j w, J s + b = b + j J w
double denom1_real = 0.0;
double denom1_imag = w;
double denom2_real = b;
double denom2_imag = J * w;
// denom = denom1 * denom2
double denom_real = denom1_real * denom2_real - denom1_imag * denom2_imag;
double denom_imag = denom1_real * denom2_imag + denom1_imag * denom2_real;
double G_real = Km * denom_real / (denom_real * denom_real + denom_imag * denom_imag);
double G_imag = -Km * denom_imag / (denom_real * denom_real + denom_imag * denom_imag);
double L_real = K * G_real;
double L_imag = K * G_imag;
double mag = Math.hypot(L_real, L_imag);
double phase = Math.toDegrees(Math.atan2(L_imag, L_real));
double magDb = 20.0 * Math.log10(mag);
System.out.println(w + " " + magDb + " " + phase);
}
}
}
More sophisticated Java code would wrap these calculations in classes that represent transfer functions and controllers, similar in spirit to MATLAB or Python. Robotics-oriented Java frameworks often expose high-level APIs for setting gains and bandwidth, while Bode and Nyquist checks are performed in a desktop design environment.
8. Wolfram Mathematica Implementation
Wolfram Mathematica has rich symbolic and numeric tools for linear
systems, including direct functions for root locus, Bode, and Nyquist
plots. A typical workflow uses
TransferFunctionModel and
plotting commands.
(* Parameters *)
J = 0.01;
b = 0.1;
Km = 1.0;
Kdes = 50.0;
s = LaplaceTransform[t, t, s]; (* symbolic placeholder, not strictly required *)
(* Plant and loop transfer functions *)
G = TransferFunctionModel[Km/(s (J s + b)), s];
L = Kdes G;
T = SystemsModelFeedbackConnect[L, 1];
(* 1) Root locus (vary scalar gain k) *)
RootLocusPlot[G, {k, 0, 200},
PlotLegends -> "Expressions",
PlotLabel -> "Root locus of G(s)"
]
(* 2) Bode magnitude and phase of L(s) *)
BodePlot[L, {10^-1, 10^2},
PlotRange -> All,
GridLines -> Automatic
]
(* 3) Nyquist plot *)
NyquistPlot[L,
PlotRange -> All,
GridLines -> Automatic
]
(* Step response of closed-loop system *)
StepResponsePlot[T, {0, 5},
PlotLabel -> "Closed-loop step response"
]
In robotics problems, Mathematica can be used to symbolically derive transfer functions for simplified joint dynamics or linearized models, after which the same commands above generate root locus and frequency-response plots for controller design.
9. Algorithmic View of Root-Locus and Frequency-Domain Design
From a numerical analysis point of view, software implementations of root locus, Bode, and Nyquist share a common structure:
- Represent \( G(s) \) as a rational function with numerator and denominator polynomials in \( s \) having real coefficients.
- For root locus, generate a grid of gain values \( K_1, \dots, K_N \), form the characteristic polynomial coefficients of \( 1 + K_i G(s) \), and compute its roots numerically for each \( i \).
- For Bode and Nyquist, generate a frequency grid \( \omega_1, \dots, \omega_M \) and evaluate \( L(j\omega_k) \) by substituting \( s = j\omega_k \).
- Compute magnitude, phase, and stability margins using standard formulas, and detect Nyquist encirclements of \( -1 \) via phase tracking.
For example, if \( G(s) = \dfrac{b_m s^m + \cdots + b_0}{a_n s^n + \cdots + a_0} \), then the characteristic polynomial at gain \( K \) is
\[ a_n s^n + \cdots + a_0 + K (b_m s^m + \cdots + b_0), \]
whose roots are computed using standard polynomial root algorithms. This is the core numerical mechanism behind the graphical root-locus plotting routines in software.
Robotics control design uses exactly these algorithms behind the scenes to tune joint gains, verify bandwidth, and ensure that all closed-loop poles lie in the left-half plane with adequate damping and that gain/phase margins exceed specified thresholds.
10. Problems and Solutions
Problem 1 (Root Locus Computation in Software): Consider the plant \( G(s) = \dfrac{1}{s(s + 4)} \) in unity feedback with proportional gain \( K \). The characteristic polynomial is \( s^2 + 4s + K \). Show that for any fixed \( K > 0 \) the closed-loop poles are given by
\[ s_{1,2}(K) = -2 \pm \sqrt{4 - K}, \]
and explain how a root-locus routine in software effectively plots these locations as \( K \) varies.
Solution:
Starting from the characteristic polynomial \( s^2 + 4s + K = 0 \), apply the quadratic formula:
\[ s = \frac{-4 \pm \sqrt{4^2 - 4 K}}{2} = \frac{-4 \pm \sqrt{16 - 4K}}{2} = -2 \pm \sqrt{4 - K}. \]
For each \( K \), there are two closed-loop poles at \( s_{1,2}(K) \). A software root-locus routine represents the polynomial coefficients \( [1, 4, K] \) numerically and computes its roots using a standard polynomial root function for a grid of \( K \) values. Plotting these roots for each \( K \) in the complex plane yields the root locus. The explicit expression above provides a closed-form description of the locus for this second-order example.
Problem 2 (Bode Margins from Software): For the same plant \( G(s) = \dfrac{1}{s(s + 4)} \) and gain \( K = 4 \), let \( L(s) = K G(s) \). Explain the steps a Bode routine takes to approximate the phase margin, and interpret the phase margin in terms of expected damping of the closed-loop response.
Solution:
The loop transfer function is \( L(s) = \dfrac{4}{s(s + 4)} \). A Bode routine proceeds as follows:
- Generate a logarithmically spaced frequency grid \( \omega_1, \dots, \omega_M \).
- For each \( \omega_k \), evaluate \( L(j\omega_k) = \dfrac{4}{j\omega_k (j\omega_k + 4)} \) and compute its magnitude \( |L(j\omega_k)| \) and phase \( \angle L(j\omega_k) \).
- Find the gain crossover frequency \( \omega_{\text{gc}} \) such that \( |L(j\omega_{\text{gc}})| \approx 1 \) (or 0 dB). This usually involves interpolation between adjacent points on the grid.
- Compute the phase margin \( \varphi_m = 180^\circ + \angle L(j\omega_{\text{gc}}) \).
A larger phase margin is associated with better damping and reduced overshoot in the closed-loop step response, assuming the loop has no RHP poles. For moderate-order plants, phase margins around \( 45^\circ \) to \( 60^\circ \) typically correspond to well-damped responses with insignificant oscillations. The software routine makes this quantitative by computing \( \varphi_m \) numerically and reporting it alongside the Bode plots.
Problem 3 (Nyquist Encirclements in Software): Suppose \( L(s) \) has no open-loop RHP poles. A Nyquist routine samples \( L(j\omega) \) on a dense grid \( \omega_k \) and plots the resulting curve in the complex plane. Explain how the routine can approximate the net number of clockwise encirclements \( N \) of \( -1 \), and why \( N = 0 \) is the condition for closed-loop stability.
Solution:
The Nyquist plot is the image of a contour enclosing the RHP under \( L(s) \). In practice, software approximates this by evaluating \( L(j\omega_k) \) for a set of frequencies covering a sufficiently large range. To approximate \( N \), the routine:
- Tracks the locus of \( L(j\omega_k) \) in the complex plane and monitors its angle with respect to the point \( -1 \).
- Computes the total change in angle of the vector from \( -1 \) to \( L(j\omega_k) \) as \( \omega_k \) increases.
- Divides the net change in angle by \( 2\pi \) to obtain the approximate number of encirclements (sign indicates clockwise or counterclockwise).
The Nyquist theorem states \( N = Z - P \), where \( Z \) is the number of closed-loop RHP poles and \( P \) is the number of open-loop RHP poles. If \( P = 0 \), then closed-loop stability (no RHP poles) requires \( Z = 0 \), hence \( N = 0 \). Thus the software checks whether the approximated Nyquist plot encircles \( -1 \) and with what net count, thereby determining stability.
Problem 4 (Software-Aided Gain Selection): For the plant \( G(s) = \dfrac{10}{s(s + 1)(s + 5)} \), consider proportional control with gain \( K \). Describe a software-based procedure (without writing full code) to find a gain \( K \) such that: (i) all closed-loop poles have real part less than \( -1 \), and (ii) the phase margin is at least \( 40^\circ \).
Solution:
Define \( L(s) = K G(s) \). A software procedure could be:
- Use a root-locus function to generate closed-loop pole locations for a grid of \( K \) values. For each \( K_i \), compute the real parts of all closed-loop poles and check whether they are strictly less than \( -1 \). Record the subset of gains satisfying this constraint.
-
For each candidate gain from step 1, compute the Bode plot of
\( L(s) \) and extract the phase margin
\( \varphi_m(K_i) \) using a
margin-type function. - Select gains with \( \varphi_m(K_i) \geq 40^\circ \). If multiple gains satisfy both constraints, choose the one that gives the desired compromise between settling time and control effort, as observed from the closed-loop step response.
This is a systematic numerical search guided by root locus (for pole location constraints) and Bode margins (for robustness), exactly mirroring graphical design procedures but executed via software.
11. Summary
This lesson connected the mathematical foundations of root locus, Bode plots, and Nyquist diagrams to their software realizations. We emphasized how the characteristic equation \( 1 + L(s) = 0 \) underlies root-locus computations, and how frequency-response sampling of \( L(j\omega) \) yields stability margins and Nyquist encirclement counts.
Implementations in Python, MATLAB/Simulink, C++, Java, and Wolfram Mathematica all follow the same numerical structure: polynomial root finding for closed-loop poles and complex evaluation of \( L(j\omega) \) on frequency grids. In robotics and mechatronic systems, these tools are routinely coupled with plant models of actuators and joints to design controllers that meet both time-domain performance specifications and frequency-domain robustness margins.
The next lessons in this chapter will build on these ideas to automate parameter sweeps, robustness studies, and reproducible reporting of linear control designs.
12. References
- Nyquist, H. (1932). Regeneration theory. Bell System Technical Journal, 11(1), 126–147.
- Evans, W. R. (1948). Graphical analysis of control systems. Transactions of the American Institute of Electrical Engineers, 67(1), 547–551.
- Bode, H. W. (1945). Network Analysis and Feedback Amplifier Design. D. Van Nostrand Company.
- MacColl, L. A. (1945). Fundamental theory of servomechanisms. Bell System Technical Journal, 24(3), 555–592.
- Zames, G. (1981). Feedback and optimal sensitivity: Model reference transformations, multiplicative seminorms, and approximate inverses. IEEE Transactions on Automatic Control, 26(2), 301–320.
- Doyle, J. C., Francis, B. A., & Tannenbaum, A. R. (1992). Feedback Control Theory. Macmillan (see chapters on classical loop shaping).
- Franklin, G. F., Powell, J. D., & Emami-Naeini, A. (2014). Feedback Control of Dynamic Systems, 7th ed. Pearson (classical design and software tools).