## Controller Stability Analysis

The purpose of controller stability analysis is to determine the range of controller gains between lower `K_{cL}` and upper `K_{cU}` limits that lead to a stable controller.

$$K_{cL} \le K_c \le K_{cU}$$

The principles of stability analysis presented here are general for any linear time-invariant system whether it is for controller design or for analysis of system dynamics. Several characteristics of a system in the Laplace domain can be deduced without transforming a system signal or transfer function back into the time domain. Some of the analysis relies on the roots of the transfer function denominator, also known as poles. The roots of the numerator, also known as zeros, do not affect the stability directly but can potentially cancel an unstable pole to create an overall stable system.

**Converge or Diverge**

A first point of analysis is whether the system converges or diverges. This is determined by analyzing the roots of the denominator of the transfer function. If any of the real parts of the roots of the denominator are positive then the system is unstable.

A simple rule to determine whether there are positive real roots is to examine the signs of the polynomial. If there are mixed signs (+ or -) then the system will be unstable because there is at least one positive real root.

Before modern computational methods, there were several methods devised to determine the stability of a system. One such approach is the Routh-Hurwitz stability criterion. The leading left edge of a table determines whether the system is stable for or any *n*th-degree polynomial

$$a_n s^n + a_{n-1} s^{n-1} + \cdots + a_1 s + a_0$$

The coefficients of the polynomial are placed into tabular form and additional coefficients *b* and *c* are computed from higher rows.

`a_n` | `a_{n-2}` | `a_{n-4}` | `\ldots` |

`a_{n-1}` | `a_{n-3}` | `a_{n-5}` | `\ldots` |

`b_{1}=\frac{a_{n-1}a_{n-2}-a_n a_{n-3}}{a_{n-1}` | `b_{2}=\frac{a_{n-1}a_{n-4}-a_n a_{n-5}}{a_{n-1}}` | `b_{3}` | `\ldots` |

`c_{1}=\frac{b_1 a_{n-3}-a_{n-1} b_{2}}{b_1}` | `c_{2}=\frac{b_1 a_{n-5}-a_{n-1} b_{3}}{b_1}` | `c_{3}` | `\ldots` |

`\vdots` | `\vdots` | `\vdots` | `\ddots` |

The terms *b* and *c* are:

$$b_i=\frac{a_{n-1}\times{a_{n-2i}}-a_n\times{a_{n-2i-1}}}{a_{n-1}}$$

$$c_i=\frac{b_1\times{a_{n-2i-1}}-a_{n-1}\times{b_{i+1}}}{b_1}$$

A changing sign (+ or -) of the leading left edge `a_n`, `a_{n-1}`, `b_{1}`, `c_{1}` indicates that the system is unstable.

Several additional methods can be used to determine stability as summarized below.

Method | Stability Condition |

Routh Array | Leading left edge terms are greater than or equal to zero |

Root Locus Plot | All poles (roots of denominator) have negative real parts |

Bode Plot | The response magnitude at -180 deg phase is less than one |

Simulation | Simulate a closed loop response with increasing or decreasing controller gains until instability occurs |

In addition to analysis in the Laplace domain, stability can be determined from a model in state space form.

$$\dot x = A x + B u$$

$$y = C x + D u$$

A state space model is stable when the eigenvalues of the *A* matrix have negative real parts.

**Oscillatory or Smooth**

A second point of analysis is whether the system exhibits oscillatory or smooth behavior. If any of the roots of the denominator have an imaginary component then the system has oscillations. Imaginary roots always come in pairs with the same positive and negative imaginary values.

**Final Value Theorem**

The Final Value Theorem (FVT) gives the steady state gain `K_p` of a transfer function `G(s)` by taking the limit as `s \to 0`

$$K_p = \lim_{s \to 0}G(s)$$

The FVT also determines the final signal value `y_\infty` for a stable system with output `Y(s)`. Note that the Laplace variable `s` is multiplied by the signal `Y(s)` before the limit is taken.

$$y_\infty = \lim_{s \to 0} s \, Y(s)$$

The FVT may give misleading results if applied to an unstable system. It is only applicable to stable systems or signals.

**Initial Value Theorem**

The Initial Value Theorem (IVT) gives an initial condition of a signal by taking the limit as `s \to \infty`. Like the FVT, the Laplace variable `s` is multiplied by the signal `Y(s)` before the limit is taken.

$$y_0 = \lim_{s \to \infty} s \, Y(s)$$

**Controller Stability**

Controller stability analysis is finding the range of controller gains that lead to a stabilizing controller. There are multiple methods to compute this range between a lower limit `K_{cL}` and an upper limit `K_{cL}`.

$$K_{cL} \le K_c \le K_{cU}$$

This range is important to know the range of tuning values that will not lead to a destabilizing controller. With modern computational tools and powerful computers, the simulation based option is frequently used for complex systems.

**Exercise**

Consider a feedback control system that has the following open loop transfer function.

$$G(s) = \frac{4K_c}{(s+1)(s+2)(s+3)}$$

Determine the values of `K_c` that keep the closed loop system response stable.

**Solution**

*Routh Array*

`a_n=1` | `a_{n-2}=11` |

`a_{n-1}=6` | `a_{n-3}=6+4 K_c` |

`b_{1}=\frac{66 - 6 - 4 K_c}{6}` | `b_{2}=0` |

`c_{1}=6+4 K_c` | 0 |

The leading edge cannot change signs for the system to be stable. Therefore, the following conditions must be met:

$$a_n=1 > 0$$

$$a_{n-1}=6 > 0$$

$$b_{1}=\frac{66 - 6 - 4 K_c}{6} > 0$$

$$c_{1}=6+4 K_c > 0$$

The positive constraint on `b_1` leads to `K_c<15`. The positive constraint on `c_1` means that `K_c > -1.5`. Therefore the following ranges are acceptable for the controller stability.

$$-1.5 < K_c < 15$$

This is a more comprehensive solution than the other methods shown below because it also includes a lower bound on the controller stability limit (if a direct acting controller were inadvertently used).

*Root Locus Plot*

Determine where the real portion of the roots crosses to the right hand side of the plane. In this case, the real part of two roots becomes positive at `K_c=15`.

*Bode Plot*

Determine the gain margin at -180^{o} phase. The magnitude at -180^{o} phase is about -23 dB. With `-23 = 20 log_{10} (AR)`, the gain margin is `1/{AR}` and approximately equal to 15. This is the upper bound on the controller gain to keep the system stable. This answer agrees with the root locus plot solution.

from scipy import signal

import matplotlib.pyplot as plt

# open loop

num = [4.0]

den = [1.0,6.0,11.0,6.0]

sys = signal.TransferFunction(num, den)

t1,y1 = signal.step(sys)

# closed loop

Kc = 1.0

num = [4.0*Kc]

den = [1.0,6.0,11.0,4.0*Kc+6.0]

sys2 = signal.TransferFunction(num, den)

t2,y2 = signal.step(sys2)

plt.figure(1)

plt.subplot(2,1,1)

plt.plot(t1,y1,'k-')

plt.legend(['Open Loop'],loc='best')

plt.subplot(2,1,2)

plt.plot(t2,y2,'r--')

plt.legend(['Closed Loop'],loc='best')

plt.xlabel('Time')

# root locus plot

import numpy.polynomial.polynomial as poly

n = 1000 # number of points to plot

nr = 3 # number of roots

rs = np.zeros((n,2*nr)) # store results

Kc = np.logspace(-2,2,n) # Kc values

for i in range(n): # cycle through n times

den = [1.0,6.0,11.0,4.0*Kc[i]+6.0] # polynomial

roots = poly.polyroots(den) # find roots

for j in range(nr): # store roots

rs[i,j] = roots[j].real # store real

rs[i,j+nr] = roots[j].imag # store imaginary

plt.figure(2)

plt.subplot(2,1,1)

plt.xlabel('Root (real)')

plt.ylabel('Root (imag)')

plt.grid(b=True, which='major', color='b', linestyle='-')

plt.grid(b=True, which='minor', color='r', linestyle='--')

for i in range(nr):

plt.plot(rs[:,i],rs[:,i+nr],'.')

plt.subplot(2,1,2)

plt.plot([Kc[0],Kc[-1]],[0,0],'k-')

for i in range(3):

plt.plot(Kc,rs[:,i],'.')

plt.ylabel('Root (real part)')

plt.xlabel('Controller Gain (Kc)')

# bode plot

w,mag,phase = signal.bode(sys)

plt.figure(3)

plt.subplot(2,1,1)

plt.semilogx(w,mag,'k-',linewidth=3)

plt.grid(b=True, which='major', color='b', linestyle='-')

plt.grid(b=True, which='minor', color='r', linestyle='--')

plt.ylabel('Magnitude')

plt.subplot(2,1,2)

plt.semilogx(w,phase,'k-',linewidth=3)

plt.grid(b=True, which='major', color='b', linestyle='-')

plt.grid(b=True, which='minor', color='r', linestyle='--')

plt.ylabel('Phase')

plt.xlabel('Frequency')

plt.show()

**Assignment**