Stabilizing an Inverted Pendulum on a Cart

Project should be handed-in on the date announced in class. Grade will be based on content as well as on clarity and on neatness of presentation.

The figure shows a rod attached to a cart through a pivot. A force f(t) (produced by a motor which is controlled by a computer) can be applied to the cart, as shown in the figure.1

cart.gif

We are using the variable ``p'' to indicate the distance from some given point along the x-axis, and the variable ``q'' to indicate the angle that the rod makes with the vertical (clockwise is positive).

It is possible to describe the complete behavior by means of two second-order equations, one for p(t) and one for q(t). Each of these two equations is of the type ``force = mass×acceleration'' as in Newton's law. Just as done with the harmonic oscillator, each of these two equations can be rewritten as a system of two first-order equations, after introducing two additional variables representing the velocity of the cart and the angular velocity of the rod. We end up with a system of 4 first-order equations.

Now, the equations are a bit messy, and they involve trigonometric functions because the tangential component of the gravity vector along the q variable is mgsinq (actually, it is even messier, because the mass of the cart appears as well). To avoid the trigs, for most of the project we will assume that the angle q » 0, so that sinq » q and therefore all nonlinear terms disappear; that is, we linearize around the equilibrium corresponding to the upwards position and stationary cart. This means that what we do will only be valid as long as we keep deviations of the rod close enough to the vertical. (Towards the end of the project, we will meet the full nonlinear equations.)

After the simplification q » 0, we end up with a system of four linear differential equations on the four variables q,q¢,p,p¢. It is convenient to collect these variables into a vector X = (x1,x2,x3,x4), where we write x1,x2,x3,x4 for the variables q,q¢,p,p¢ respectively. With the help of a friendly physics major, we obtain the following model, written for convenience in vector form:

dX
dt
(t)   =    æ
ç
ç
ç
ç
ç
ç
ç
ç
ç
ç
è
0
1
0
0
(M+m)g
ML
0
0
0
0
0
0
1
- mg
M
0
0
0
ö
÷
÷
÷
÷
÷
÷
÷
÷
÷
÷
ø
 X(t)       +      f(t) æ
ç
ç
ç
ç
ç
ç
ç
ç
ç
ç
è
0
-1
ML
0
1
M
ö
÷
÷
÷
÷
÷
÷
÷
÷
÷
÷
ø
  ,
where g = 9.8m/sec2 is the acceleration due to gravity, M and m are the masses of the cart and rod respectively, and L is the length of the rod.

These equations are easy to interpret: the first and third just say that x1¢ = x2 and x3¢ = x4, which are obvious. The second equation x2¢(t) = [((M+m)g)/( ML)]x1(t)-[1/( ML)]f(t) tells us that Mq¢¢ = cq- df(t), for some constants c and d. That is, ``mass times acceleration'' equals the ``net force'' on the q direction (to be precise, since we are dealing with rotations, we should really say ``moment of inertia'' and ``torques''). The first force term ``cq'' represents the effect of gravity (remember that mgsinq is approximated by mgq, because q » 0) and the second term represents the force being applied by the motor. The negative sign of the coefficient ``[(-1)/( ML)]'' of f(t) can be understood as follows: if we push to the right (as in the arrow representing f in the picture), the rod wants to move counterclockwise, and if we instead apply a negative f (``pull'' to the left), then the rod will want to move clockwise.

Part 1: Interpret the last equation (the one for x4 = p¢¢). Write a one-paragraph informal explanation, like we did for the other variables. (If you have trouble with the first term, keep reading and then come back to this point.)

From now on, to be concrete, let us suppose that M = 5kg, m = 0.5kg, and L = 1m. So the equations become:

dX
dt
(t)   =    æ
ç
ç
ç
ç
ç
è
0
1
0
0
10.78
0
0
0
0
0
0
1
- 0.98
0
0
0
ö
÷
÷
÷
÷
÷
ø
 X(t)     +     f(t) æ
ç
ç
ç
ç
ç
è
0
- 0.2
0
0.2
ö
÷
÷
÷
÷
÷
ø
  .

Now let us ask what happens if we start the cart with the rod slightly displaced to the right, say with q(0) = 0.1 radians (how many degrees is that?). We first assume that no force f is applied (the motor is turned off). We also suppose, to make things easier, that we start with zero velocity (we are holding the rod, and drop it when we hear ``start the experiment now''), and that at time t = 0 the car is stationary at p = 0.

Part 2: Solve (numerically, with Maple) the differential equation on the interval from t = 0 to t = 1, plot q(t) and p(t), and interpret your result.

Help: This is one set of Maple commands that lets you do what is required:

with(DEtools):with(plots):
sys:=diff(x1(t),t)=x2(t),
     diff(x2(t),t)=10.78*x1(t),
     diff(x3(t),t)=x4(t),
     diff(x4(t),t)=-0.98*x1(t);
init:=x1(0)=0.1,x2(0)=0,x3(0)=0,x4(0)=0;
sol:=dsolve({sys,init},{x1(t),x2(t),x3(t),x4(t)},type=numeric):
odeplot(sol,[[t,x1(t)],[t,x3(t)]],0..1,numpoints=100);

More help: To interpret the graph obtained for p(t), you may want to perform the following experiment. Get a sharp pencil, or even better a ball-point pen, and hold it at an angle on a tabletop. Now drop it. What happens with the tip? Do you see how the tip moves in the opposite direction?

pencil.gif

(You may need to start with the pencil or pen very slanted, so that there is not too much friction.)

Recall that, for two dimensional systems, if some eigenvalue has a positive real part, then most trajectories become unbounded. (Of course, positive real numbers are a special case of complex numbers with positive real part, so this includes all three cases, namely real sources, saddles, and spiral sources.) Actually, the same is true in any dimension, not just for systems of two equations: if some eigenvalue has a positive real part, then most trajectories become unbounded. From the plot, then, the answer that you get next should not be surprising.

Part 3 Calculate the eigenvalues of the matrix which defines the system (when there is no force). Discuss.

Help: This is one set of Maple commands that lets you do what is required:

with(linalg):
M:=array(1..4,1..4,[[0,1,0,0],[10.78,0,0,0],[0,0,0,1],[-0.98,0,0,0]]);
eigenvals(M);

Now we start with something that is much more interesting. Our objective will be to program the computer (that is, to decide what forces must be applied) so as to bring the cart to the position p » 0 and also make the rod be vertical (q » 0). We also want that once we get to p » 0 and q » 0, we stay at that position, with a vertical rod, that is, we want to also achieve p¢ = [(dp)/( dt)] » 0 and q¢ = [(dq)/( dt)] » 0. In summary, we want to program things so that all four variables making up the vector X(t) converge to zero as t becomes large, no matter where we start.

Let us be a bit more mathematical. We wish to find coefficients c1,c2,c3,c4 with this property: if we substitute

f(t) = c1x1(t)+c2x2(t)+c3x3(t)+c4x4(t)
into the equation for the system (we get x2¢ = (10.78-0.2c1)x1 -0.2(c2x2+c3x3+c4x4), and so on), then all the solutions of the new system have the property that they converge to zero.2

Practically, this means that our computer is going to measure q(t) etc, it will calculate an appropriate force f(t), evaluating some expression which involves these measured values, and it will then then tell the motor to apply this force f(t). As we want f(t) to be applied at time t (not seconds later), this the calculation must be performed very quickly. (This is no problem in practice.3)

This problem is not easy! (As a challenge, try finding a good set of coefficients without peeking at what is written next.)

Let's understand a bit of theory. For two-dimensional systems, the origin is a sink precisely if all the eigenvalues have negative real part (of course, negative real numbers are a special case of complex numbers with negative real part, so this includes both real and spiral sinks). In other words, for all solutions of the differential equation to have the property that their coordinates approach zero as t®¥, it is required that all eigenvalues have negative real part. It turns out that this fact is true for systems of any dimension, not just two:

all solutions go to zero if and only if all the eigenvalues have negative real part

so our objective can be expressed in terms of eigenvalues.

After consultation with our friendly control theory neighbor4 we are told that the force should be computed as follows:

f(t) = 150x1(t)+40x2(t)+8x3(t)+12x4(t)  .

Hmmm, is she right?

Part 4: Calculate the eigenvalues of the new matrix which is obtained after substituting this formula for f. Show that all the eigenvalues now have negative real part, so all solutions should indeed converge to zero.

The matrix obtained for the ``closed loop'' system after substituting the formula for f is, by the way:

æ
ç
ç
ç
ç
ç
è
0
1
0
0
-19.22
-8.0
-1.6
-2.4
0
0
0
1
29.02
8.0
1.6
2.4
ö
÷
÷
÷
÷
÷
ø

Just to convince ourselves that this works, let us simulate one solution:

Part 5: Solve (numerically, with Maple) the differential equations (now with the force substituted-in), with x1(0) = 0.1, x2(0) = 0, x3(0) = 0.5, and x4(0) = 0, on the time interval t = 0 to t = 10. Interpret your results.

(This time, I will not help you with the Maple code. Figure it out yourself. It is a simple modification of the commands given above; just change the coefficients as appropriate.)

Finally, let us go back to the first assumption that we made. We studied the linearized problem. The full nonlinear model, as our physics-major friend tells us after sweating a bit, is like this:5

q¢¢   =    mL(q¢)2sinqcosq-(M+m)gsinq+f(t)cosq
mLcos2q-(M+m)L
p¢¢   =    mgsinqcosq-mL(q¢)2sinq-f(t)
mcos2q-(M+m)
in terms of accelerations. (I did tell you it was messy, remember?)

Part 6: Rewrite these as a system of four first-order equations (in the variables xi as before). Remember to use M = 5, m = 0.5, and L = 1.

The ``feedback law'' f(t) = 150x1(t)+40x2(t)+8x3(t)+12x4(t) works well, when the initial conditions are small, for the true system:

Part 7: Plot the solution of the nonlinear differential equation which you obtained in Part 6, and substituting the formula given for f, with x1(0) = 0.1, x2(0) = 0, x3(0) = 0.5, and x4(0) = 0, on the time interval t = 0 to t = 10. Interpret your results.

However, it fails when we are farther from the linearization point.

Part 8: Plot the solution of the nonlinear differential equation which you obtained in Part 6, after inserting the formula given above for f, with x1(0) = 1, x2(0) = 0, x3(0) = 0.5, and x4(0) = 0, on the time interval t = 0 to t = 10. Interpret your results.

(Warning: it took my computer almost a minute to come up with the solution. Be patient.)

As a challenge, you may want to see if you can find different values of the coefficients of the control law (instead of c1 = 150, etc). Or even come up with a non linear feedback law, that stabilizes the system for a larger set of initial conditions.

The control of the (full nonlinear model) inverted pendulum is actually a challenging problem even in advanced control theory, and there are lots of questions that are still very poorly understood!


Footnotes:

1 This is the same problem as balancing a stick on the end of one's finger, and, more interestingly, the problem that must be solved when balancing a rocket during liftoff. Our Controls Lab has one of these cart/rod systems - you may want to ask for a demo sometime.

2 Of course, we could also look for a nonlinear f(t) as a function of the system variables, but things are already interesting enough with linear functions.

3 At our lab, a Pentium-133 processor can easily handle doing all the required operations about 100 times per second, i.e. ``at sampling frequency 100Hz''.

4 Or, perhaps, after buying and reading a control book, for example ``Mathematical Control Theory'' by E. Sontag.

5 Actually, everything we've done assumed that friction effects were negligible. This is actually not true in practice.


File translated from TEX by TTH, version 2.00.
On 20 Apr 1999, 16:00.
This page was last updated on September 05, 2006 at 10:28 am and is maintained by webmaster@math.rutgers.edu.
For questions regarding courses and/or special permission, please contact mclausen@math.rutgers.edu.
For questions or comments about this site, please contact help@math.rutgers.edu.
© 2012 Rutgers, The State University of New Jersey. All rights reserved.