RouthStability

Note: Video lecture available for this section!

Authors: John D'Arcy, Matt Hagen, Adam Holewinski, and Alwin Ng

Stewards: Jeff Falta, Taylor Lebeis, Shawn Mayfield, Marc Stewart, Tom Welch

Date created: 10/26/2006; updated: 10/25/2007
 * First round reviews for this page
 * Rebuttal for this page

Introduction
The stability of a process control system is extremely important to the overall control process. System stability serves as a key safety issue in most engineering processes. If a control system becomes unstable, it can lead to unsafe conditions. For example, instability in reaction processes or reactors can lead to runaway reactions, resulting in negative economic and environmental consequences.

The absolute stability of a control process can be defined by its response to an external disturbance to the system. The system may be considered stable if it exists at a consistent state or setpoint and returns to this state immediately after a system disturbance. In order to determine the stability of a system, one often must determine the eigenvalues of the matrix representing the system’s governing set of differential equations. Unfortunately, sometimes the characteristic equation of the matrix (the polynomial representing its eigenvalues) can be difficult to solve; it may be too large or contain unknown variables. In this situation, a method developed by British mathematician Edward Routh can yield the desired stability information without explicitly solving the equation.

Recall that in order to determine system stability one need only know the signs of the real components of the eigenvalues. Because of this, a method that can reveal the signs without actual computation of the eigenvalues will often be adequate to determine system stability.

To quickly review, negative real eigenvalue components cause a system to return to a steady state point (where all partial derivatives equal zero) when it experiences disturbances. Positive real components cause the system to move away from a stable point, and a zero real component indicates the system will not adjust after a disturbance. Imaginary components simply indicate oscillation with a general trend in accordance with the real part. Using the method of Routh stability, one can determine the number of each type of root and thus see whether or not a system is stable. When unknown variables exist in the equation, Routh stability can reveal the boundaries on these variables that keep the roots negative.

The Routh Array
The Routh array is a shortcut to determine the stability of the system. The number of positive (unstable) roots can be determined without factoring out any complex polynomial.

Generating the Array
The system in question must have a characteristic equation of a polynomial nature. as shown below:

$$P(S) = a_{n}S^n + a_{n-1}S^{n-1} + \cdots + a_{1}S+ a_0 \qquad $$

In order to examine the roots, set P(S)=0, which will allow you to tell how many roots are in the left-hand plane, right hand plane, and on the j-omega axis. If the system involves trigonometric functions it needs to be fit to a polynomial via a Taylor series expansion. One necessary condition for stability is that $$ a_n > 0 \qquad $$. (If $$ a_n < 0 \qquad $$, all coefficients may be multiplied by -1 before checking). The other condition is that all values in column 1 of the Routh array must be positive for the system to be stable.



This flow diagram shows the generation of a Routh array for an idealized case with m,n representing the location in the matrix.

The coefficients of the polynomial are placed into an array as seen below. The number of rows is one more than the order of the equation. The number of sign changes in the first column indicate the number of positive roots for the equation.

In the array, the variables b1,b2,c1,c2,etc. are determined by calculating a determinant using elements from the previous two rows as shown below:

$$b_1 = \frac{a_{n-1}a_{n-2} - a_na_{n-3}}{a_{n-1}}, b_2 = \frac{a_{n-1}a_{n-4} - a_na_{n-5}}{a_{n-1}}, b_3 = \frac{a_{n-1}a_{n-6} - a_na_{n-7}}{a_{n-1}},\cdots$$

$$c_1 = \frac{b_1a_{n-3} - a_{n-1}b_2}{b_1}, \ c_2 = \frac{b_1a_{n-5} - a_{n-1}b_3}{b_1}, c_3 = \frac{b_1a_{n-7} - a_{n-1}b_4}{b_1},\cdots$$

The general expression for any element $$ x $$ after the first two rows with index (m,n) is as follows:

For $$A=\begin{bmatrix}x_{m-2,1} & x_{m-2,n+1}\\ x_{m-1,1} & x_{m-1,n+1}\end{bmatrix}$$

$$ x_{m,n} = \frac{-det(A)}{x_{m-1,1}} $$

Note, that if the Routh array starts with a zero, it may still be solved (assuming that all the other values in the row are not zero), by replacing the zero with a constant, and letting that constant equal a very small positive number. Subsequent rows within that column that have this constant will be calculated based on the constant choosen.

Once the array is complete, apply the following theorems to determine stability:

1) If all of the values in the fist column of the Routh array are >0, then P(S) has all negative real roots and the system is stable.

2) If some of the values in the first column of the Routh array are <0, then the number of times the sign changes down the first column will = the number of positive real roots in the P(S) equation.

3) If there is 1 pair of roots on the imaginary axis, both the same distance from the origin (meaning equidistant), then check to see if all the other roots are in the left hand plane. If so, then the imaginary roots location may be found using $$ AS^2 + B = 0 $$, where A and B are the elements in the Routh array for the 2nd to last row.

To clarify even further, an example with real numbers is analyzed.

Example Array
The following polynomial was generated from a sample system.

$$P(S) = 5S^3 -10S^2 + 7S+ 20 \qquad $$

The preceding polynomial must be investigated in order to determine the stability of the system. This is done by generating a Routh array in the manner described above. The array as a result of this polynomial is,

In the array shown above, the value found in the third row is calculated as follows. $$b_1 = \frac{{-10*7} - {5*20}}{-10}$$ The array can now be analyzed. When looking down the first column, it can be seen that 5 is positive in magnitude, then the sign changes in the -10 entry, and the sign changes a second time to positive 17. This counts as two changes in sign, which corresponds to two positive roots, making the system unstable.

Finding Stable Control Parameter Values
Often, for a unit operation, a PID parameter such as controller gain (Kc), the integral time constant (Ti), or the derivative time constant (Td) creates an additional variable in the characteristic equation. This can be carried through the computations of the Routh array to indicate which values of the variable will provide stability to the system through by preventing positive roots from occuring in the equation. For example, if a controller output is governed by the function:

$$ 10s^3 + 5s^2 + 8s + (T_d + 2)  \qquad $$

The stable values of Td can be found via a Routh array:

We reveal $$ -2 < T_d < 2 \qquad $$ in order to keep the first column elements positive, so this is the stable range of values for this parameter. If multiple parameters were in the equation, they would simply be solved for as a group, yielding constraints along the lines of "Ti + Kc > 2" etc, so any value chosen for one parameter would give a different stable range for the other.

Special Cases
There are a few special cases that one should be aware of when using the Routh Test. These variances can arise during stability analysis of different control systems. When a special case is encountered, the traditonal Routh stability solution methods are altered as presented below.

One of the coefficients in the characteristic equation equals zero
If the power of the 0*Sn is $$ \ge 1$$, replace the zero with a quantity, ε, which would be positive and will approach zero. Then continue with your analysis as normal. Essentially this gives the limit of the roots as that coefficient approaches zero. (If the power = 0, see Case 3). For example,

Equation:	$$2S^3 - 24S + 32 = 0 \qquad $$

Working Equation:    $$2S^3 + \epsilon S^2 - 24S + 32 = 0 \qquad $$

Since ε is positive we know that in the first column row 2 will be positive, row 4 will be positive, and row 3 will be negative. This means we will have a sign change from 2 to 3 and again from 3 to 4. Because of this, we know that two roots will have positive real components. If you actually factor out the equation you see that $$(S-2)^2 (2S + 8) = 0\qquad$$, showing that we do have 2 positive roots. Both of these roots are equal to 2, so there is technically only one root, but in any case we know the system is unstable and must be redesigned.

One of the roots is zero
This case should be obvious simply from looking at the polynomial. The constant term will be missing, meaning the variable can be factored from every term. If you added an &epsilon; to the end as in case 1, the last row would be &epsilon; and falsely indicate another sign change. Carry out Routh analysis with the last zero in place. Equation:	$$S^3 - S^2 - 2S = 0 \qquad $$

As you can see in column one we have row 1 positive, row 2 and 3 negative, and row 4 zero. This is interpreted as one sign change, giving us one positive real root. Looking at this equation in factored form,

$$(S + 1)(S - 2) S = 0 \qquad $$

we can see that indeed we have only one positive root equal which equals 2. The zero in the last row indicates an additional unstable root of zero. Alternatively, you may find it easier to just factor out the variable and find the signs of the remaining eigenvalues. Just remember there is an extra root of zero.

A row full of zeros
When this happens you know you have either a pair of imaginary roots, or symmetric real roots. The row of zeros must be replaced. The following example illustrates this procedure.

Equation:	$$S^4 - 6S^3 + 10S^2 -6S +9 = 0 \qquad $$

Row 4 contains all zeros. To determine its replacement values, we first write an auxiliary polynomial A determined by the entries in Row 3 above.

$$ A(S) = 9S^2 + 9 \qquad $$

Notice that the order decreases by 1 as we go down the table, but decreases by 2 as we go across.

We then take the derivative of this auxiliary polynomial.

$$ A'(S) = 18S \qquad $$

The coefficients obtained after taking the derivative give us the values used to replace the zeros. From there, we can proceed the table calcuations normally. The new table is

In fact, the purely imaginary or symmetric real roots of the original polynomial are the same as the roots of the auxiliary polynomial. Thus, we can find these roots.

$$ 9S^2 + 9 = 0 \qquad $$

$$ S = \pm i\qquad $$

Because we have two sign changes, we know the other two roots of the original polynomial are positive.

In fact, after factoring this polynomial, we obtain

$$(S^2 + 1)(S-3)^2 = 0 \qquad $$

Therefore, the roots are $$S = \pm i, 3 \qquad $$, where in this case, the root 3 has multiplicity 2.

Limitations
Routh arrays are useful for classifying a system as stable or unstable based on the signs of its eigenvalues, and do not require complex computation. However, simply determing the stability is not usually sufficient for the design of process control systems. It is important to develop the extent of stability as well as how close the system is to instability. Further stability analysis not accounted for in the Routh analysis technique include finding the degree of stability, the steady state performance of the control system, and the transient response of the system to disturbances.

More involved techniques, such as those discussed in Eigenvalues and Eigenvectors, must be used to further characterize the system stability (with the exception of system polynomials resulting in pure imaginary roots). Another limitation of the Routh method occurs when the polynomial in question becomes so large that Routh stability is too computationally time consuming (a personal judgment). For this situation another method, such as a root locus plotmust be used.

Note that for defining stability, we will always start out with a polynomial. This polynomial arises from finding the eigenvalues of the linearized model. Thus we will never encounter other functions, say exponenential functions or sin or cos functions in general for stability analysis in control theory.

Advantages Over Root Locus Plots
Routh stability evaluates the signs of the real parts of the roots of a polynomial without solving for the roots themselves. The system is stable if all real parts are negative. Therefore unlike root locus plots, the actual eigenvalues do not need to be calculated for a Routh stability analysis. Furthermore, sometimes the system has too many unknowns to easily construct and interpret a root locus plot (e.g. with two PID controllers there are the variables Kc1, Kc2, τi1, τi2, τd1, and τd2).

Example 1
Assume the following polynomial describes the eigenvalues of a linearized model of your process. For this polynomial, complete a Routh array and determine the system's stability?

$$ P(x) = x^4 + 10x^3 + 35x^2 + 50x + 24 \qquad $$

Answer Since P(X) is a fourth-order polynomial, the Routh array contains five rows.

Rows 1 and 2 correspond to the coefficients of the polynomial terms:

Rows 3, 4, and 5 contain the determinants using elements from the previous two rows.

Therefore,

The complete Routh array:

Since all the values in the first column are positive, the equation P(x) has all negative roots and the system is stable.

Example 2
Consider a system with the following characteristic equation:

$$ 20s^3+59s^2+46s+(4+K_c)=0 \qquad$$

Using a P-only controller, find the range of controller gain that will yield a stable system.

Answer

Since the equation is a third-order polynomial, the Routh array has four rows. Routh Array:

For the system to be stable, all the elements of the first column of the Routh array have to be positive.

The first column will contain only positive elements if:
 * $$46 - \frac{20}{59}(4+K_c)>0$$
 * $$K_c<\frac{2634}{20}$$


 * $$ 4+K_c>0 \qquad$$
 * $$K_c>-4 \qquad$$

Acceptable stable range $$-4<K_c<\frac{2634}{20}$$

Example 3
Consider a system with the following characteristic equation: $$ s^5-3s^4+s^3+s^2+4=0 \qquad$$ Determine the stability of this system.

Answer One of the coefficients in the characteristic equation equals 0. We replace the zero with a quantity which would be positive (approach zero from the right-hand side) and continue with the analysis as normal. Working equation: $$ s^5-3s^4+s^3+s^2+ \epsilon s+4=0 \qquad$$

Since ε is positive, in the first column, there are two sign changes, from row 1 to row 2 and from row 2 to row 3. Thus, we know that the roots will have two positive real components. If you actually factor out the equation you will see that, $$(s^2+1)(s-2)^2(s+1)=0 \qquad$$, showing that we do have two positive roots, both equal to 2. Additional complication exists because at row 5, as ε goes to zero, the term also goes to zero, which means that for row 5, we are getting a row full of zeros. This means that we have a pair of imaginary roots, and this situation can be solved using the equation, $$ a_2S^2 + a_0 = 0 \qquad $$.

In this case, the working equation is, $$S^2+4=0 \qquad $$ The imaginary roots are, $$ S = \pm 2 i \qquad $$

Example 4
You are an engineer at an icecream factory. There is a storage vat with a cooling system that has a PI controller scheme. It has the following characteristic equation: $$ 10s^4+3s^3-K_cs^2+3s+(T_i+6)=0 \qquad$$ Your job is to determine the constraints on the values $$K_c$$ and $$T_i$$ such that the system is stable.

Answer The goal is to make the matrix such that the first column has no sign changes. Since the first two entries in the first column are numbers and positive, therefore all other values in this column must be positive. Working equation: $$ 10s^4+3s^3-K_cs^2+3s+(T_i+6)=0 \qquad$$

Since $$10-K_c\,$$ is in the first column, it gives the constraint $$K_c<10\,$$. Likewise, $$T_i+6\,$$ must be positive, giving the constraint $$T_i>-6\,$$. Which brings the final entry of $$3- \frac{3(T_i+6)}{(10-K_c)}$$. By dividing by three and rearrangement of the terms, it is seen that $$ \frac{(T_i+6)}{(10-K_c)} < 1$$. The terms are then rearranged to arrive at the inequality: $$T_i+K_c<4\,$$.

Thus the stable conditions for the constants are: $$K_c<10\,$$ $$T_i>-6\,$$ $$T_i+K_c\,<4$$

Sage's Corner
Routh Stability Anaylsis Example



Routh Stability Slides