Working with MATLAB

Uncovered Idiosyncrasies

Every engineering tool—from your favorite oscilloscope to a handy programming environment—has its own set of quirks that you find annoying. In this article, Shlomo discusses the idiosyncrasies of MathWorks’s MATLAB and provides tips you’ll find useful the next time you use it.

ROOT LOCUS DIAGRAM

Consider the system of Figure 1 but assume for the moment that the gain, K, is one. When looking at the system, you see what seems to be a very precise description of how the system is constructed and how its parts can be characterized. As an engineer, you know that despite the apparent precision, the description is only an approximation of the real system.

One way in which an engineer can try to make sure that the system will actually behave in a reasonable way is to vary the gain, K, and see how changes in K affect the system. The root locus diagram shows the location of the poles of the closed-loop system as a function of changes in K. (For more about the root locus diagram and how to plot it, you may wish to refer to my book, A Mathematical Introduction to Control Theory, 2^nd Ed., 2015.) Every point in the diagram is a pole of the closed-loop system for some value of K > 0 (but the values of K are not displayed in the plot). If all the points in the diagram have negative real parts, then the closed-loop system is stable for all positive gains, and you know that the system’s gain margin is infinite. Additionally, the root locus diagram can be used when one is designing a system to see how adding gain to a system affects the position of the poles of the closed-loop system.

PREDICT OVERSHOOT THE MATLAB WAY

To use MATLAB to plot the root locus diagram of the system in Figure 1, use the following command:

rlocus(tf([1 2],[1 1 0])

This causes MATLAB to plot the root locus of the system (see Figure 2). Using the MATLAB data cursor, we examine the point at which the branches of the root locus reenter the real axis and find that at that point the gain, K, is 5.83 and, according to MATLAB, there is no overshoot (see Figure 2). It is simple enough to have MATLAB plot the step response of the closed-loop system when K is 5.83. The step response of the closed-loop system as plotted by MATLAB is shown in Figure 3. The step response seems to have overshoot. How can that be?

• Figure 2
• Figure 3

PREDICT OVERSHOOT THE OLD-FASHIONED WAY

We now analyze the root locus diagram and the system itself to determine whether or not the system really suffers from overshoot. We know that:

Branches can enter or leave the real axis at real values of s for which:

The solutions of this equation are:

Considering Figure 2, it is clear that the point of reentry is s = –2 – √2. We find that:

Thus, the corresponding gain, K, is:

So far everything agrees with MATLAB. Let’s move on to examining the step response of the closed-loop system. The transfer function of the closed-loop system, H(s), is:

We find that the system seems to be critically damped. It has two identical real poles, and we recall that critically damped systems do not suffer from overshoot—or do they?

Actually, that rule about critically damped systems, which MATLAB seems to be applying, only applies when the numerator of the second order transfer is a constant. What really happens here?

The function H(s) is the Laplace transform of the impulse response which is also the derivative of the step response (as long as the degree of the numerator of H(s) is less than that of its denominator). Let’s find the derivative of the step response to see what it tells us.

We find that:

And, therefore:

where u(t) is the unit-step function. Clearly, this function is positive until t = 1/√2 and is negative afterwards. This shows that the step response increases initially and then decreases to its final value. Moreover, the maximum of the step response occurs when t = 1/√2 ≈ 0.71 s. Thus, the plot of the step response is correct, and the 0% overshoot that the data cursor reports on is wrong—and is seemingly a consequence of misapplying a rule of thumb.

SYSTEM ANALYSIS

Next we consider the root locus associated with the system of Figure 4, a unity-feedback system for which:

The root locus plot associated with the system is given in Figure 5, and a blowup of the plot for low gains—where we should see all six visible branches starting from –1—is given too. Unfortunately, the branches do not start from –1. They start near –1.

The elementary properties of the root locus—that require that the branches of the root locus begin at poles of the open-loop system—guarantee that MATLAB is making a mistake. To understand why MATLAB does not work properly here, consider what happens when MATLAB is asked to find the poles of G(s)—the zeros of the denominator of G(s). The zeros of the denominator are the solutions of the equation:

• Figure 4
• Figure 5

By inspection, we find that the solutions are –1, –1, –1, –1, –1, –1, and 0 and are the starting points of the seven branches of the root locus.

When MATLAB is asked to solve this equation (using its roots command), it responds with:

These numbers are close to the exact solutions but are not the exact solutions (except for the initial zero). On the other hand, comparing the solutions near –1 with the starting points of the branches as given in the right-hand figure in Figure 5, we see that the agreement is fairly good. It seems that the problem here is that MATLAB does not “see” the exact solutions of the polynomial equation even though they are obvious to us.

A SIMPLE PROBLEM

It is an easily proved fact, and one that is made use of in the development of Rader’s FFT, that when considering numbers modulo a prime number the sequence of the powers of any nonzero number, x, the sequence x¹, x², …, x^{p – 1}, must end in one and the sequence of powers of x must be periodic with period p – 1 or some factor of p – 1. Additionally, none of sequence’s elements can be zero.^[1]

It should be easy to check that this property holds by having MATLAB evaluate the command:

mod(x.^(1:p-1),p)

It should be easy, but this does not always work. Consider, for example, x = 11, and p = 17. In this case, MATLAB responds with:

This sequence ends with a zero and cannot be correct, but where is the problem?

Looking back at the command, the problem is clear: the 11 is being raised to each of the powers 1, …, 16. The number 11¹⁶ is very large—too large for MATLAB to store it as an integer. Once it is not stored as an integer, there is some imprecision in its value, and this is enough to lead what should have been a one to become a zero.

A simple method for causing MATLAB to calculate the values correctly is to use the sequence of commands in Listing 1. When calculating the sequence in this way, cur(k) is always between 0 and p – 1. Thus, the numbers here remain in the range that MATLAB stores as integers. Once again, it is necessary to understand the tool with which we are working and the problem we are solving to make sure that MATLAB really does what we need it to do.

Listing 1
A simple method for causing MATLAB to calculate the values correctly

cur(1) = x;
for k = 2:p-1
    cur(k) = mod(cur(k-1)*x,p);
end

OVERCOME SHORTCOMINGS

I love MATLAB. It makes it possible to perform calculations I would not want to perform by hand and that, in many cases, I could not calculate manually. I use it as much as possible in my work because it is a wonderful tool. On the other hand, like every other tool, you must learn its quirks, and you must be prepared to use your knowledge of the problem at hand to compensate for the shortcomings of the tool you are using. 

REFERENCE
[1] Wikipedia, “Multiplicative Group of Integers Modulo n,” https://en.wikipedia.org/wiki/Multiplicative_group_of_integers_modulo_n.

RESOURCE
S. Engelberg, A Mathematical Introduction to Control Theory, 2nd ed., Imperial College Press, London, England, 2015.

SOURCES
MATLAB
MathWorks | www.mathworks.com

PUBLISHED IN CIRCUIT CELLAR MAGAZINE • JULY 2016 #312 – Get a PDF of the issue

Shlomo Engelberg

+ posts

Shlomo Engelberg (shlomoe@jct.ac.il) earned Bachelor’s and Master’s degrees in Engineering from The Cooper Union. He holds a PhD in Mathematics from New York University’s Courant Institute. He is the Dean of the School of Engineering and Computer Science at the Jerusalem College of Technology. Shlomo is the author of many articles and several books, the latest of which are ADuC841 Microcontroller Design Manual: From Microcontroller Theory to Design Projects (Circuit Cellar, 2011) and A Mathematical Introduction to Control Theory, 2nd Ed. (Imperial College Press, 2015). Shlomo’s research interests include applied mathematics, instrumentation and measurement, signal processing, coding theory, and control theory.

This author does not have any more posts.