Example: Solution to the inverted pendulum problem using Root Locus method

The transfer function of the plant for this problem is given below:

where,

The design criteria (with the pendulum receiving a 1N impulse force from the cart) are:

• Settling time of less than 5 seconds.
• Pendulum should not move more than 0.05 radians away from the vertical.

To see how this problem was originally set up, consult the inverted pendulum modeling page.

## Transfer functions

The rlocus command in MATLAB can find the root locus for a system described by state-space equations or by a transfer function. For this problem it will be easier in the long run to use a transfer function (the reason for this will become clear later). The transfer function found from the Laplace transforms for the output Phi (the pendulum's angle) can be set up using MATLAB by entering the numerator and denominator as vectors. Create an m-file (or a '.m' file located in the same directory as MATLAB) and copy the following text to model the transfer function:

```
M = .5;
m = 0.2;
b = 0.1;
i = 0.006;
g = 9.8;
l = 0.3;

q = (M+m)*(i+m*l^2)-(m*l)^2;   %simplifies input

num = [m*l/q  0];
den = [1  b*(i+m*l^2)/q  -(M+m)*m*g*l/q  -b*m*g*l/q];
pend=tf(num,den);
```

### Block Diagram

The control of this problem is a little different than the standard control problems you may be used to. Since we are trying to control the pendulum's position, which should return to the vertical after the initial disturbance, the reference signal we are tracking should be zero. The force applied to the cart can be added as an impulse disturbance. The schematic for this problem should look like the following.

It will be easier to determine the appropriate transfer function to enter into MATLAB if we first rearrange the schematic as follows:

## Root locus design

This closed-loop transfer function can be modeled in MATLAB. The first thing to do is to look at the root locus for the plant by itself, with no compensator. To pick the gain for the proportional control, remember that the settling time has to be less than 5 seconds. This implies that sigma should be more than 4.6/5=0.92. We can put this criteria right on the root locus using the sigrid function, which has to be copied to a m-file. To do this, copy the following code to the end of your m-file :

```rlocus(pend)
sigrid(0.92)
axis([-6 6 -6 6])
```
You should see the following rlocus plot:

As you can see, one of the roots of the closed-loop transfer function is in the right-half-plane. This means that the system will be unstable. Look back at the root locus to see why. Part of the root locus lies between the origin and the pole in the right-half-plane. No matter what gain you chose, you will always have a closed-loop pole in this region, making your impulse response unstable. To solve this problem, we need to add another pole at the origin so all the zeros and poles at the origin will cancel each other out and multiple roots will be created in this right-half-plane region. The multiple roots can then be drawn into the left-half-plane to complete the design. Add the following to your m-file:

```contr=tf(1,[1 0]);
rlocus(contr*pend)
sigrid(0.92)
axis([-10 10 -10 10])
```
You should get the following root locus plot with multiple roots in the right-half-plane:

Now we can begin trying to draw the branches of the root locus into the left half plane. Enter the following two commands into the command window:

```[num,den]=tfdata(contr*pend,'v');
roots(num)
roots(den)
```
You should get the following output in the MATLAB command window,
```ans =
0

ans =
0
-5.6041
5.5651
-0.1428```
As you can see there are four poles and only one zero. This means there will be three asymptotes: one along the real axis in the negative direction, and the other two at 120 degree angles from this.

This configuration will never have the multiple roots in the left-half-plane. We must reduce the number of asymptotes from three to two by adding one more zero than pole the controller. If just a zero is added, then the intersection of the asymptotes (alpha) would be [(-5.6041+5.5651-0.1428+0+0)-(0+0+z2)]/2. This means the two asymptotes will leave the real axis at roughly -0.1-(1/2)z2. Making z2 as small as possible (assume near the origin) will not pull the multiple roots far enough into the left-half-plane to meet the design requirements (asymptotes will leave at about -0.1).

The solution to this problem is to add another pole far to the left of the other poles and zeros. To keep the right number of asymptotes, another zero should be added as well. The placement of the added poles and zeros is not important except that the pole should be relatively large and the zeros should be relatively small.

Try the m-file below to see what effect the poles and zeros have on the root locus.

```M = .5;
m = 0.2;
b = 0.1;
i = 0.006;
g = 9.8;
l = 0.3;

q = (M+m)*(i+m*l^2)-(m*l)^2;   %simplifies input

num = [m*l/q  0];
den = [1  b*(i+m*l^2)/q  -(M+m)*m*g*l/q  -b*m*g*l/q];
pend=tf(num,den);

z1 = 3;
p1 = 0;
z2 = 4;
p2 = 50;
numlag = [1 z1];
denlag = [1 p1];
contr=tf(numc,denc);

rlocus(contr*pend)
sigrid(0.92)
axis([-50 50 -50 50])
figure
rlocus(contr*pend)
sigrid(0.92)
axis([-10 10 -10 10])
[k,poles]=rlocfind(contr*pend)
sys_cl=feedback(pend,k*contr);

figure
impulse(sys_cl)
axis([0 2 -0.05 0.05])
```

Note: With some versions of MATLAB, this root locus does not work correctly due to the way that the vector of gains is automatically chosen. If this happens with your version, you should go to the MathWorks home page and download the bug fix. Or, you can work around it for now and tell MATLAB which gains you would like it to use. In this case, the following code can be used to replace the rlocus command.
```gain=0:10:2000;
rlocus(contr*pend,gain)
```
The root locus will plot with "x" instead of a connected line, but the rlocfind command still works appropriately.

The poles and zeros were found by trial and error. The only things to keep in mind was that one pole had to be at the origin, the other had to be far to the left, and the two zeros had to be small. Furthermore, if the two zeros are close together and to the right of the farthest left plant pole, the response is better. You should see first the following root locus plot with the zeros and poles listed above:

The second plot should be of the same root locus magnified a little so that the root locus around the origin can be seen.

When prompted to pick a location on the root locus, chose a spot on the multiple roots just before they return to the real axis. Your velocity response to the impulse disturbance should look similar to the following:

The response now meets all of the requirements, so no further iteration is needed.

## What happens to the cart's position?

At the beginning on this solution page, the block diagram for this problem was given. The diagram was not entirely complete. The block representing the the position was left out because that part was not being controlled. It would be interesting though, to see what is happening to the cart's position when the controller for the pendulum's angle is in place. To see this, we need to consider the actual system block diagram:

Rearranging a little bit, you get the following block diagram:

The feedback loop represents the controller we have designed for the pendulum. The transfer function from the cart's position to the impulse force, with the feedback controller which we designed, is given as follows:

Now that we have the transfer function for the entire system, let's take a look at the response. First we need the transfer function for the cart's position. To get this we need to go back to the Laplace transforms of the system equations and find the transfer function from U(s) to X(s). Below is this transfer function:

where,

For more about the Laplace transform please refer to the inverted pendulum modeling page.

Now, add the following to your m-file and run it in the command window:

```num2 = [(i+m*l^2)/q  0  -m*g*l/q];
den2 = [1  b*(i+m*l^2)/q  -(M+m)*m*g*l/q  -b*m*g*l/q  0];
G2=tf(num2,den2);

figure;
t=0:0.01:6;
subplot(2,1,1);
impulse(sys_cl,t)
axis([0 6 -0.05 0.05])
subplot(2,1,2);
xpos=feedback(1,k*contr*pend)*G2;
impulse(xpos,t)
axis([0 6 -0.1 0.1])
```
You should see the following plot:

The top curve represents the pendulum's angle, and the bottom curve represents the cart's position. As you can see, the cart moves at first, then is stabilized at near zero for almost five seconds, and then goes unstable. It is possible that friction (which was neglected in the modeling of this problem) will actually cause the cart's position to be stabilized. Keep in mind that if this is in fact true, it is due more to luck than to anything else, since the cart's position was not included in the control design.

Root Locus Examples
Cruise Control | Motor Speed | Motor Position | Bus Suspension | Inverted Pendulum | Pitch Controller | Ball & Beam

Inverted Pendulum Examples
Modeling | PID | Root Locus | Frequency Response | State Space | Digital Control | Simulink

Tutorials
MATLAB Basics | MATLAB Modeling | PID Control | Root Locus | Frequency Response | State Space | Digital Control | Simulink Basics | Simulink Modeling | Examples