1. Stability in the frequency domain

The frequency domain allows us to find the stability of closed loop systems using only open loop transfer functions and simple operations.

This material is also covered in this video. The GeoGebra sheet is available here.

import numpy
from matplotlib import pyplot as plt
%matplotlib inline

1.1. Locating poles and zeros of a complex function

Let’s construct a complex transfer function by specifying the poles, zeros and gain separately.

zeros = [1]
poles = [-1 + 1j, -1 - 1j]
gain = 1

from numpy.polynomial.polynomial import polyvalfromroots

def G(s):
    return gain*polyvalfromroots(s, zeros)/polyvalfromroots(s, poles)

It will be useful for us to be able to plot a complex curve easily

def plotcomplex(curve, color='blue', marker=None):
    plt.plot(numpy.real(curve), numpy.imag(curve), color=color, marker=marker)

def plotpz():
    for p in poles:
        plotcomplex(p, color='red', marker='x')
    for z in zeros:
        plotcomplex(z, color='red', marker='o')

This function will change the axes to be a cross through the origin and have an equal aspect ratio (so that a circle appears as a circle)

from tbcontrol.plotting import cross_axis

Let’s construct a circular contour and see how the image of the contour moves around as the contour moves around. The image is \(G(s)\) as \(s\) goes through a countour

from ipywidgets import interact

def plotsituation(contour):
    plotcomplex(contour)
    plotcomplex(G(contour), color='red')
    plotpz()
    cross_axis()

theta = numpy.linspace(0, 2*numpy.pi, 1000)

def argumentprinciple(centerreal=(-2., 2.), centerimag=(-2., 2.), radius=(0.5, 3)):
    contour = radius*numpy.exp(1j*theta) + centerreal + 1j*centerimag
    plotsituation(contour)

interact(argumentprinciple)

You should be able to verify the Cauchy argument principle using the interaction above:

As \(s\) describes a simple contour enclosing \(N_p\) poles and \(N_z\) zeros, the image \(G(s)\) encircles the origin \(w = N_z - N_p\) times. \(w\) is the winding number.

1.2. Closed loop stability

Normally we will be looking at transfer functions of the form

\[\frac{GK}{1 + GK}\]

So we will want to check if the denominator of the above \((1 + GK)\) has roots in the RHP. To do this we can construct a special contour called the Nyquist D contour which encloses the whole of the RHP. It starts at the origin, then goes up to infinity, circles around at infinite distance from the origin in a clockwise direction, and then comes back up the imaginary axis. For most functions, the part at infinity just maps \(1 + GK\) to 1 as \(GK\) goes to zero as s goes to infinity.

omega = numpy.logspace(-2, 2, 1000)
Dcontour = numpy.concatenate([1j*omega, -1j*omega[::-1]]) # We're ignoring the infinite arc

Let’s assume that \(K=1\) and check if our system will be closed loop stable

K = 1

plotcomplex(K*G(Dcontour) + 1)
cross_axis(size=2)

Counting encirclements of the origin of \(1 + GK\) is the same as counting encirclements of \(-1\) by \(GK\):

def nyquistplot(K):
    plotcomplex(K*G(Dcontour))
    plotcomplex(-1, color='red', marker='o')
    cross_axis(size=2)

nyquistplot(K=1)

This enables us to reason easily about the effect of the controller gain on stability:

interact(nyquistplot, K=(0.5, 5.))

1.3. Nyquist stability criterion

Let \(N_P\) be the number of poles of KG(s) encircled by the D contour and \(N_Z\) be the number of zeros of \(1+KG(s)\) encircled by the D contour. \(N_Z\) is the number of poles of the closed loop system in the right half plane. The resultant image shall encircle (clock-wise) the point \((-1+j0)\) \(w\) times such that \(w = N_Z - N_P\).

For a stable \(G\) this boils down to spotting when the Nyquist plot encircles the -1 point.

1.4. Bode stability criterion

Nyquist plots are hard to draw by hand, though, so we often use the Bode stability criterion instead. This works by noticing that, in order for the Nyquist graph to encircle the -1 point, the phase angle must reach -180 ° and the magnitude must be bigger than 1. We can draw a Bode diagram and a Nyquist diagram next to each other to see the effect of changing gains.

def bodeplot(K):
    fig = plt.figure(figsize=(10,5))

    ax_gain = plt.subplot2grid((2, 2), (0, 0))
    ax_phase = plt.subplot2grid((2, 2), (1, 0))
    ax_complex = plt.subplot2grid((2, 2), (0, 1), rowspan=2)

    freqresp = K*G(1j*omega)

    ax_gain.loglog(omega, numpy.abs(freqresp))
    ax_gain.axhline(1, color='orange')
    ax_gain.set_ylim([0.1, 10])
    ax_gain.set_ylabel('|G|')

    ax_phase.semilogx(omega, numpy.unwrap(numpy.angle(freqresp)) - numpy.angle(freqresp[0])) # We know the angle should start at 0
    ax_phase.axhline(-numpy.pi, color='green')
    ax_phase.set_ylabel('∠G / rad')
    ax_phase.set_xlabel('ω / (rad/s)')

    plt.sca(ax_complex)
    nyquistplot(K)

    circle = numpy.exp(-1j*numpy.linspace(0, numpy.pi*2))
    ax_complex.plot(circle.real, circle.imag, color='orange')
    ax_complex.plot([-2, 0], [0, 0], color='green', linewidth=4, alpha=1, zorder=-1)

interact(bodeplot, K=(0.5, 5.))