Notebook

CDS 110, Lecture 7

Frequency Domain Analysis using Bode/Nyquist plots

Richard M. Murray, Winter 2024

The purpose of this lecture is to introduce tools that can be used for frequency domain modeling and analysis of linear systems. It illustrates the use of a variety of frequency domain analysis and plotting tools.

In [ ]:

# Import standard packages needed for this exercise
import numpy as np
import matplotlib.pyplot as plt
import math
try:
  import control as ct
  print("python-control", ct.__version__)
except ImportError:
  !pip install control
  import control as ct

# Use ctrlplot defaults for matplotlib
plt.rcParams.update(ct.rcParams)

Stable system: servomechanism¶

We start with a simple example a stable system for which we wish to design a simple controller and analyze its performance, demonstrating along the way the basic frequency domain analysis functions in the Python control toolbox (python-control).

Consider a simple mechanism for positioning a mechanical arm whose equations of motion are given by

$$ J \ddot \theta = -b \dot\theta - k r\sin\theta + \tau_\text{m}, $$

which can be written in state space form as

$$ \frac{d}{dt} \begin{bmatrix} \theta \\ \theta \end{bmatrix} = \begin{bmatrix} \dot\theta \\ -k r \sin\theta / J - b\dot\theta / J \end{bmatrix} + \begin{bmatrix} 0 \\ 1/J \end{bmatrix} \tau_\text{m}. $$

The system consists of a spring loaded arm that is driven by a motor, as shown below.

The motor applies a torque that twists the arm against a linear spring and moves the end of the arm across a rotating platter. The input to the system is the motor torque $\tau_\text{m}$. The force exerted by the spring is a nonlinear function of the head position due to the way it is attached.

The system parameters are given by

$$ k = 1,\quad J = 100,\quad b = 10, \quad r = 1,\quad l = 2,\quad \epsilon = 0.01, $$

and we assume that time is measured in msec and distance in cm. (The constants here are made up and don't necessarily reflect a real disk drive, though the units and time constants are motivated by computer disk drives.)

The system dynamics can be modeled in python-control using a NonlinearIOSystem object, which we create with the nlsys function:

In [ ]:

# Parameter values
servomech_params = {
    'J': 100,             # Moment of inertia of the motor
    'b': 10,              # Angular damping of the arm
    'k': 1,               # Spring constant
    'r': 1,               # Location of spring contact on arm
    'l': 2,               # Distance to the read head
    'eps': 0.01,          # Magnitude of velocity-dependent perturbation
}

# State derivative
def servomech_update(t, x, u, params):
    # Extract the configuration and velocity variables from the state vector
    theta = x[0]                # Angular position of the disk drive arm
    thetadot = x[1]             # Angular velocity of the disk drive arm
    tau = u[0]                  # Torque applied at the base of the arm

    # Get the parameter values
    J, b, k, r = map(params.get, ['J', 'b', 'k', 'r'])

    # Compute the angular acceleration
    dthetadot = 1/J * (
        -b * thetadot - k * r * np.sin(theta) + tau)

    # Return the state update law
    return np.array([thetadot, dthetadot])

# System output (end of arm)
def servomech_output(t, x, u, params):
    l = params['l']
    return np.array([l * x[0]])

# System dynamics
servomech = ct.nlsys(
    servomech_update, servomech_output, name='servomech',
    params=servomech_params,
    states=['theta_', 'thdot_'],
    outputs=['y'], inputs=['tau'])

print(servomech)
print("\nParams:", servomech.params)

Linearization¶

To study the open loop dynamics of the system, we compute the linearization of the dynamics about the equilibrium point corresponding to $\theta_\text{e} = 15^\circ$.

In [ ]:

# Convert the equilibrium angle to radians
theta_e = (15 / 180) * np.pi

# Compute the input required to hold this position
u_e = servomech.params['k'] * servomech.params['r'] * np.sin(theta_e)
print("Equilibrium torque = %g" % u_e)

# Linearize the system about the equilibrium point
P = servomech.linearize([theta_e, 0], u_e, name='P_ss')
P.name = 'P_ss'  # TODO: fix in nlsys_improvements
print("Linearized dynamics:", P)
print("Zeros: ", P.zeros())
print("Poles: ", P.poles())
print("")

# Transfer function representation
P_tf = ct.tf(P, name='P_tf')
print(P_tf)

Open loop frequency response¶

A standard method for understanding the dynamics is to plot the output of the system in response to sinusoids with unit magnitude at different frequencies.

We use the frequency_response function to plot the step response of the linearized, open-loop system.

In [ ]:

# Reset the frequency response label to correspond to a time unit of ms
ct.set_defaults('freqplot', freq_label="Frequency [rad/ms]")

# Frequency response
freqresp = ct.frequency_response(P, np.logspace(-2, 0))
freqresp.plot()

# Equivalent command
ct.bode_plot(P_tf, np.logspace(-2, 0), '--')

Feedback control design¶

We next design a feedback controller for the system using a proportional integral controller, which has transfer function

$$ C(s) = \frac{k_\text{p} s + k_\text{i}}{s} $$

We will learn how to choose $k_\text{p}$ and $k_\text{i}$ more formally in W9. For now we just pick different values to see how the dynamics are impacted.

In [ ]:

kp = 1
ki = 1

# Create tf from numerator/denominator coefficients
C = ct.tf([kp, ki], [1, 0], name='C')
print(C)

# Alternative method: define "s" and use algebra
s = ct.tf('s')
C = ct.tf(kp + ki/s, name='C')
print(C)

In [ ]:

# Loop transfer function
L = P * C
cplt = ct.bode_plot([P, C, L], label=['P', 'C', 'L'])
cplt.set_plot_title("PI controller for servomechanism")

Note that L = P * C corresponds to addition in both the magnitude and the phase.

Nyquist analysis¶

To check stability (and eventually robustness), we use the Nyquist criterion.

In [ ]:

fig = plt.figure(figsize=[7, 4])
ax1 = plt.subplot(2, 2, 1)
ax2 = plt.subplot(2, 2, 3)
ct.bode_plot(L, ax=[ax1, ax2])

# Tidy up the figure a bit
fig.align_labels()
ax1.set_title("Bode plot for L")

ax2 = plt.subplot(1, 2, 2)
ct.nyquist_plot(L, ax=ax2, title="")
plt.title("Nyquist plot for L")

plt.suptitle("Loop analysis for (unstable) servomechanism")
plt.tight_layout()

We see from this plot that the loop transfer function encircles the -1 point => closed loop system should be unstable. We can check this by making use of additional features of Nyquist analysis.

In [ ]:

# Get the Nyquist *response*, so that we can get back encirclements
nyqresp = ct.nyquist_response(L)
print("N = encirclements: ", nyqresp.count)
print("P = RHP poles of L: ", np.sum(np.real(L.poles()) > 0))
print("Z = N + P = RHP zeros of 1 + L:", np.sum(np.real((1 + L).zeros()) > 0))
print("Zeros of (1 + L) = ", (1 + L).zeros())
print("")

T = ct.feedback(L)
ct.step_response(T).plot(
    title="Step response for (unstable) servomechanism",
    time_label="Time [ms]");

Poles on the $j\omega$ axis¶

Note that we have a pole at 0 (due to the integrator in the controller). How is this handled?

A: use a small loop to the right around poles on the $j\omega$ axis => not inside the contour.

To see this, we use the nyquist_response function, which returns the contour used to compute the Nyquist curve. If we zoom in on the contour near the origin, we see how the outer edge of the Nyquist curve is computed.

In [ ]:

fig = plt.figure(figsize=[7, 5.8])

# Plot the D contour
ax1 = plt.subplot(2, 2, 1)
plt.plot(np.real(nyqresp.contour), np.imag(nyqresp.contour))
plt.axis([-1e-4, 4e-4, 0, 4e-4])
plt.xlabel('Real axis')
plt.ylabel('Imaginary axis')
plt.title("Zoom on D-contour")

# Clean up the display of the units
from matplotlib import ticker
ax1.xaxis.set_major_formatter(ticker.StrMethodFormatter("{x:.0e}"))
ax1.yaxis.set_major_formatter(ticker.StrMethodFormatter("{x:.0e}"))

ax2 = plt.subplot(2, 2, 2)
ct.nyquist_plot(L, ax=ax2)
plt.title("Nyquist curve")

plt.suptitle("Nyquist contour for pole at the origin")
plt.tight_layout()

Second iteration feedback control design¶

We now redesign the control system to give something that is stable. We can do this by moving the zero for the controller to a lower frequency, so that the phase lag from the integrator does not overlap with the phase lag from the system dynamics.

In [ ]:

# Change the frequency response to avoid crossing over -180 with large gain
Cnew = ct.tf(kp + (ki/200)/s, name='C_new')
Lnew = ct.tf(P * Cnew, name='L_new')

plt.figure(figsize=[7, 4])
ax1 = plt.subplot(2, 2, 1)
ax2 = plt.subplot(2, 2, 3)
ct.bode_plot([Lnew, L], ax=[ax1, ax2], label=['L_new', 'L_old'])

# Clean up the figure a bit
ax1.loglog([1e-3, 1e1], [1, 1], 'k', linewidth=0.5)
ax1.set_title("Bode plot for L_new, L_old", size='medium')

ax3=plt.subplot(1, 2, 2)
ct.nyquist_plot(Lnew, max_curve_magnitude=5, ax=ax3)
ax3.set_title("Nyquist plot for Lnew", size='medium')

plt.suptitle("Loop analysis for (stable) servomechanism")
plt.tight_layout()

We see now that we have no encirclements, and so the system should be stable.

Note however that the Nyquist curve is close to the -1 point => not that stable.

In [ ]:

# Compute the transfer function from r to y
Tnew = ct.feedback(Lnew)
cplt = ct.step_response(Tnew).plot(time_label="Time [ms]")
cplt.set_plot_title("Step response for (stable) spring-mass system")

Third iteration feedback control design (via loop shaping)¶

To get a better design, we use a PID controller to shape the frequency response so that we get high gain at low frequency and low phase at crossover.

In [ ]:

# Design parameters
Td = 1                    # Set to gain crossover frequency
Ti = Td * 10              # Set to low frequency region
kp = 500                  # Tune to get desired bandwith

# Updated gains
kp = 150
Ti = Td * 5; kp = 150

# Compute controller parmeters
ki = kp/Ti
kd = kp * Td

# Controller transfer function
ctrl_shape = kp + ki / s + kd * s

# Frequency response (open loop) - use this to help tune your design
ltf_shape = ct.tf(P_tf * ctrl_shape, name='L_shape')

cplt = ct.frequency_response([P, ctrl_shape]).plot(label=['P', 'C_shape'])
cplt = ct.frequency_response(ltf_shape).plot(margins=True)

cplt.set_plot_title("Loop shaping design for servomechanism controller")

In [ ]:

# Compute the transfer function from r to y
T_shape = ct.feedback(ltf_shape)
cplt = ct.step_response(T_shape).plot(
    time_label="Time [ms]",
    title = "Step response for servomechanism with PID controller")

Closed loop frequency response¶

We can also look at the closed loop frequency response to understand how different inputs affect different outputs. The gangof4 function computes the standard transfer functions:

In [ ]:

cplt = ct.gangof4(P_tf, ctrl_shape)

Stability margins¶

Another standard set of analysis tools is to identify the gain, phase, and stability margins for the system:

Gain margin: the maximimum amount of additional gain that we can put into the loop and still maintain stability.
Phase margin: the maximum amount of additional phase (lag) that we can put into the loop and still maintain stability.
Stability margin: the maximum amount of combined gain and phase at the critical frequency that can be put into the loop and still maintain stability.

The first two of the items can be computed either by looking at the frequency response or by using the margin command.

The stabilty margin is the minimum distance between -1 and $L(jw)$, which is just the minimum value of $|1 - L(j\omega)|$.

In [ ]:

plt.figure(figsize=[7, 4])

# Gain and phase margin on Bode plot
ax1 = plt.subplot(2, 2, 1)
plt.title("Bode plot for Lnew, with margins")
ax2 = plt.subplot(2, 2, 3)
ct.bode_plot(Lnew, ax=[ax1, ax2], margins=True)

# Compute gain and phase margin
gm, pm, wpc, wgc = ct.margin(Lnew)
print(f"Gm = {gm:2.2g} (at {wpc:.2g} rad/ms)")
print(f"Pm = {pm:3.2g} deg (at {wgc:.2g} rad/ms)")

# Compute the stability margin
resp = ct.frequency_response(1 + Lnew)
sm = np.min(resp.magnitude)
wsm = resp.omega[np.argmin(resp.magnitude)]
print(f"Sm = {sm:2.2g} (at {wsm:.2g} rad/ms)")

# Plot the Nyquist curve
ax3 = plt.subplot(1, 2, 2)
ct.nyquist_plot(Lnew, ax=ax3)
plt.title("Nyquist plot for Lnew [zoomed]")
plt.axis([-2, 3, -2.6, 2.6])

#
# Annotate it to see the margins
#

# Gain margin (special case here, since infinite)
Lgm = 0
plt.plot([-1, Lgm], [0, 0], 'k-', linewidth=0.5)
plt.text(-0.9, 0.1, "1/gm")

# Phase margin
theta = np.linspace(0, 2 * math.pi)
plt.plot(np.cos(theta), np.sin(theta), 'k--', linewidth=0.5)
plt.text(-1.3, -0.8, "pm")

# Stability margin
Lsm = Lnew(wsm * 1j)
plt.plot([-1, Lsm.real], [0, Lsm.imag], 'k-', linewidth=0.5)
plt.text(-0.4, -0.5, "sm")

plt.suptitle("")
plt.tight_layout()

Unstable system: inverted pendulum¶

When we have a system that is open loop unstable, the Nyquist curve will need to have encirclements to be stable. In this case, the interpretation of the various characteristics can be more complicated.

To explore this, we consider a simple model for an inverted pendulum, which has (normalized) dynamics:

$$ \dot x = \begin{bmatrix} 0 & 1 & \\ -1 & 0.1 \end{bmatrix} x + \begin{bmatrix} 0 \\ 1 \end{bmatrix} u, \qquad y = \begin{bmatrix} 1 & 0 \end{bmatrix} x $$

Transfer function for the system can be shown to be

$$ P(s) = \frac{1}{s^2 + 0.1 s - 1}. $$

This system is unstable, with poles $\sim\pm 1$.

In [ ]:

P = ct.tf([1], [1, 0.1, -1])
P.poles()

PD controller¶

We construct a proportional-derivative (PD) controller for the system,

$$ u = k_\text{p} e + k_\text{d} \dot{e} $$

which is roughly the equivalent of using state feedback (since the system states are $\theta$ and $\dot\theta$).

In [ ]:

# Transfer function for a PD controller
kp = 10
kd = 2
C = ct.tf([kd, kp], [1])

# Loop transfer function
L = P * C
L.name = 'L'
print(L)
print("Zeros: ", L.zeros())
print("Poles: ", L.poles())

In [ ]:

# Bode and Nyquist plots
plt.figure(figsize=[7, 4])
ax1 = plt.subplot(2, 2, 1)
plt.title("Bode plot for L", size='medium')
ax2 = plt.subplot(2, 2, 3)
ct.bode_plot(L, ax=[ax1, ax2])

ax3 = plt.subplot(1, 2, 2)
ct.nyquist_plot(L, ax=ax3)
plt.title("Nyquist plot for L", size='medium')

plt.suptitle("Loop analysis for inverted pendulum")
plt.tight_layout()

In [ ]:

# Check the Nyquist criterion
nyqresp = ct.nyquist_response(L)
print("N = encirclements: ", nyqresp.count)
print("P = RHP poles of L: ", np.sum(np.real(L.poles()) > 0))
print("Z = N + P = RHP zeros of 1 + L:", np.sum(np.real((1 + L).zeros()) >= 0))
print("Poles of L = ", L.poles())
print("Zeros of 1 + L = ", (1 + L).zeros())
print("")

T = ct.feedback(L)
ct.initial_response(T, X0=[0.1, 0]).plot();

Note that we get a warning when we set the initial condition. This is because T is a transfer function and so it doesn't have a unique state space realization. If the initial state is zero this doesn't matter, but if the initial state is nonzero then the assignment of states is not well defined.

Gang of 4¶

Another useful thing to look at is the transfer functions from noise and disturbances to the system outputs and inputs:

In [ ]:

ct.gangof4(P, C);

We see that the response from the input $r$ (or equivalently noise $n$) to the process input is very large for large frequencies. This means that we are amplifying high frequency noise (and comes from the fact that we used derivative feedback).

High frequency rolloff¶

We can attempt to resolve this by "rolling off" the derivative action at high frequencies:

In [ ]:

Cnew = (kp + kd * s) / (s/20 + 1)**2
Cnew.name = 'Cnew'
print(Cnew)

Lnew = P * Cnew
Lnew.name = 'Lnew'

plt.figure(figsize=[7, 4])
ax1 = plt.subplot(2, 2, 1)
ax2 = plt.subplot(2, 2, 3)
ct.bode_plot([Lnew, L], ax=[ax1, ax2])
ax1.loglog([1e-1, 1e2], [1, 1], 'k', linewidth=0.5)
ax1.set_title("Bode plot for L, Lnew", size='medium')

ax3 = plt.subplot(1, 2, 2)
ct.nyquist_plot(Lnew, ax=ax3)
ax3.set_title("Nyquist plot for Lnew", size='medium')

plt.suptitle("Stability analysis for inverted pendulum")
plt.tight_layout()

While not (yet) a very high performing controller, this change does get rid of the issues with the high frequency noise:

In [ ]:

# Check the gang of 4
ct.gangof4(P, Cnew);

In [ ]:

# See what the step response looks like
Tnew = ct.feedback(Lnew)
ct.step_response(Tnew, 10).plot()

In [ ]: