Can you balance a ruler upright on the palm of your hand? If I concentrate, I can just barely manage it by constantly reacting to the small wobbles of the ruler. This challenge is analogous to a classic problem in the field of control systems design: stabilizing an upside-down (“inverted”) pendulum.

One of the best things about Mathematica is that it makes specialist areas like control systems accessible to non-specialists. This lets you freely combine and develop new ideas without needing to be an expert in everything. It also makes Mathematica a great platform for learning and exploring new areas.

Using the new control systems features (one of several new application areas integrated into Mathematica 8), I’ve been experimenting with models of stabilized inverted pendulums. I’m no expert in control theory, but you’ll see that one doesn’t need to be.

Inverted pendulums arise as models of many real-world systems: Segway vehicles (rocking forward and backward), bicycles (rocking side-to-side), unicycles (both), rockets (during liftoff), and even skyscrapers.

In this post we’ll model a single inverted pendulum on a moving cart, shown above, and use the new control systems functionality to stabilize it. In later posts in this series, we’ll generalize the problem, looking at double (and triple and quadruple) pendulums, pendulums in 3D, and more.

First we need to simulate the pendulum-and-cart system, and then we can try to determine a “control” force f(t) that will move the cart to keep the pendulum balanced upright.

Assuming the pendulum arm is uniform (so its center of mass is at d/2), ignoring friction, and using dimensionless units in which the gravitational acceleration g = 1, the model pendulum-and-cart system shown above satisfies this pair of differential equations:

(We will derive these equations programmatically for a more general case in the next post. To see an outline of how to derive this case directly, download the Computable Document Format [CDF] file.)

In this post, we’ll explore the concrete case m = c = d = 1:

Using NDSolve, we can immediately simulate this system from a slightly off-center initial condition θ(0) = π/2 – 0.1. We’ll initially use no control force: f(t) = 0; that is, we make no attempt to move the cart to keep the pendulum upright.

Here’s an animation of the solution:

(For the definition of the function AnimatePendulum, download the CDF file.)

We can combine the above two steps into a handy function that animates the system for any control force f(t):

You can use the function SimulatePendulum to explore ideas for keeping the pendulum balanced upright.

The simplest plausible idea I can think of is to push the cart in whichever direction the pendulum is currently leaning—move the cart under the pendulum, to stop it falling over. The pendulum is perfectly upright when θ(t) = π/2, so let’s try pushing the cart with the force f(t) = π/2 – θ(t):

That wasn’t very good. The cart barely moved, until it was too late. Let’s try a stronger force.

Not bad! It’s better than me balancing a ruler on my hand, but it’s still a bit wobbly.

Maybe we can make it smoother by also taking into account how fast the pendulum is falling, which is θ′(t):

That stopped the wobbling nicely, but now the whole system is coasting steadily off to the right. What if we want the cart to stay at the point x(t) = 0?

We could go on modifying the control force by hand, but it’s much easier to turn to Mathematica‘s new control systems design features.

We want to keep the state of the system near the unstable equilibrium point:

The equilibrium is “unstable” because, in the absence of any control force, small deviations from this point tend to get larger.

The hand-made control forces we tried above are special cases of “linear control”, in which the control force is some constant linear combination of the deviation of the states from the equilibrium point:

The constants k_i are called “feedback gains”, and we can use a couple of new control theory functions in Mathematica 8 to find good values for them.

We begin by constructing a linear model for the behavior of the system near the equilibrium point. The function StateSpaceModel does this automatically if we give it the nonlinear equations:

In the third argument to StateSpaceModel, we specified that there is one control variable: f(t).

The result of StateSpaceModel looks like an augmented matrix. Its columns describe the effect of small deviations of each of the state variables {x(t), x ′ (t), θ(t), θ′ (t)} from equilibrium. For example, the third column says that, near the equilibrium point, small deviations in θ(t) away from the upright position (θ(t) = π/2) lead to an increase in both x′(t) (the speed of the cart) and θ′(t) (the rate at which the pendulum is falling). The augmented fifth column similarly describes the immediate effect of the control force on the state variables.

Now that we have a model for the system, we can use the function LQRegulatorGains to find values for the feedback gains k_i:

In the second argument of LQRegulatorGains we specified the following quadratic “cost” function for the system:

LQRegulatorGains finds values for the constants k_i that minimize this cost function over time (in the linearized model).

The cost coefficients {1, 10, 10, 100} are a guess, based on an intuition that deviations in θ′(t) are relatively more important than deviations in x′(t) and θ(t), which are in turn more important than deviations in x(t). Different coefficients would give better, worse, or just qualitatively different behavior. In real life the coefficients might represent actual costs associated with deviations of a system from the desired state.

Using these gains, our control force is:

Let’s test it out:

The cart is brought smoothly back to the equilibrium point!

Let’s make it more challenging by bumping the cart around a little bit. We’ll give the cart a push to the right at t = 8 and a push to the left at t = 16:

The bumping force looks like this:

Let’s simulate the bumped pendulum with our control force:

Our stabilized pendulum recovers nicely from each bump. Compare this with the performance of our hand-made control force f(t) = 3(π/2 – θ(t)) – θ'(t), which was able to keep the pendulum upright, albeit with a steady drift to the right:

You can see that our hand-made control force isn’t very robust. The first bump knocks it right over and sends it spiraling out of control.

Our pendulum stabilized via LQRegulatorGains can survive some quite bumpy conditions:

If you bump it hard enough, the linear control system we have used here won’t be able to recover. You can try that out by downloading the CDF file for this blog post. Another natural thing to explore would be the effect of different cost functions in LQRegulatorGains. It’s also fun to try to design your own superior control force by hand for this relatively simple system.

In future posts in this series, we’ll look at some generalizations of this problem, including double (and triple and quadruple) pendulums and 3D pendulums. Stay tuned!

Download CDF file