2D Physics of the Inverted Pendulum on a Cart

Equations of motion and a Python simulation of a uniform rod on a cart.

Introduction

Balancing a broom on your finger takes some practice. Here, we model the same problem with an inverted pendulum on a cart.

When a broom balanced on your finger starts leaning right, you move your hand right to catch it. If you keep your hand still, gravity makes the lean grow until the broom falls. An inverted pendulum on a cart behaves in much the same way. The cart takes the place of your hand, and we can push it left or right to try to keep the pendulum upright.

Before we can simulate this, we need to work out how the cart and pendulum respond to a push. Instead of modelling a full 3D system, let’s make our life easier and consider a 2D system. The cart moves horizontally, and the pendulum swings in the same plane around a joint called the pivot. We’ll derive the equations of motion for this system, then use them to build a simulation we can try to balance.

Equations of motion

Defining the geometry and symbols

Figure 1: Geometry and forces on the cart and rod. The arrows marked N and −N act on different objects. Resistance is drawn for a cart moving right; the cart’s vertical forces are omitted.

In Figure 1, \(x\) represents the cart’s horizontal position, measured to the right from a fixed point on the track. \(\theta\) represents the angle between the rod and the upward vertical, measured clockwise in radians.

The cart has mass \(M\). We model the pendulum as a uniform rod of mass \(m\) and length \(L\), with its mass spread evenly along its length. Its center of mass is therefore at its midpoint, a distance \(l=L/2\) from the pivot.

From the dashed triangle in Figure 1, the center of mass is offset horizontally from the pivot by \(l\sin\theta\) and vertically by \(l\cos\theta\). Its horizontal position \(x_c\) and height \(y_c\) above the pivot are then

\[ x_c=x+l\sin\theta,\qquad y_c=l\cos\theta. \tag{1}\]

A dot denotes a rate of change with time. The cart’s velocity and acceleration are \(\dot{x}\) and \(\ddot{x}\), and the rod’s angular velocity and acceleration are \(\dot{\theta}\) and \(\ddot{\theta}\).

Table 1: Symbols used in the drawing, equations, and simulation.

Symbol	Unit	What it means
\(x\)	\(\mathrm{m}\)	Cart position, measured from a fixed zero; right is positive
\(\theta\)	\(\mathrm{rad}\)	Rod angle from upright; clockwise is positive
\(\dot{x}\), \(\ddot{x}\)	\(\mathrm{m/s}\), \(\mathrm{m/s^2}\)	Cart velocity and acceleration
\(\dot{\theta}\), \(\ddot{\theta}\)	\(\mathrm{rad/s}\), \(\mathrm{rad/s^2}\)	Rod angular velocity and angular acceleration
\(M\), \(m\)	\(\mathrm{kg}\)	Cart mass and rod mass, respectively
\(L\), \(l\)	\(\mathrm{m}\)	Full rod length and pivot-to-center distance, \(l=L/2\)
\(x_c\), \(y_c\)	\(\mathrm{m}\)	Horizontal position and height of the rod’s center of mass
\(F\)	\(\mathrm{N}\)	Force we apply to the cart; right is positive
\(N\), \(P\)	\(\mathrm{N}\)	Pivot forces on the rod: horizontal (right positive) and vertical (up positive)
\(b\)	\(\mathrm{N\,s/m}\)	Cart resistance coefficient; the resisting force is \(-b\dot{x}\)
\(g\)	\(\mathrm{m/s^2}\)	Gravitational acceleration, \(9.81\ \mathrm{m/s^2}\)
\(I\)	\(\mathrm{kg\,m^2}\)	Rod’s moment of inertia about its center of mass
\(\tau_{\mathrm{net}}\)	\(\mathrm{N\,m}\)	Net torque about the rod’s center of mass; clockwise is positive
\(t\), \(\Delta t\)	\(\mathrm{s}\)	Time and the duration of one simulation step

Forces on the cart

We want to find how the cart accelerates when we push it and the rod pushes back. The cart moves on a horizontal track, which supports it against the vertical forces. Since the cart cannot move up or down, we only need to model its horizontal motion.

The horizontal forces in Figure 1 are the applied force \(F\), friction \(-b\dot{x}\), and the force \(-N\) from the rod. We model friction as proportional to the cart’s velocity, with coefficient \(b\). The minus sign makes it oppose the motion.

At the pivot, the cart exerts a horizontal force \(N\) on the rod. The rod pushes back with the equal and opposite force \(-N\), which is the force we include in the cart’s equation. Newton’s second law gives

\[ M\ddot{x}=F-b\dot{x}-N. \tag{2}\]

The force from the rod is still unknown. To determine it, we need to describe how the rod moves.

Forces on the rod

The rod can move and turn at the same time. We’ll first describe how its center of mass moves, then calculate how its angle changes. Its lower end stays attached to the cart, so the pivot moves horizontally with the cart while the center of mass can also rise and fall.

In Figure 1, the cart exerts the horizontal force \(N\) and the vertical force \(P\) on the rod at the pivot. Gravity \(mg\) acts downward at the center of mass. We neglect air resistance and friction in the pivot.

A push at one end can move the whole rod as well as turn it. Newton’s second law uses the sum of all external forces on the rod to determine the acceleration of its center of mass. The points where those forces act will also matter when we calculate the rotation. Adding the horizontal and vertical forces separately gives

\[ N=m\ddot{x}_c,\qquad P-mg=m\ddot{y}_c. \tag{3}\]

These equations tell us how to find the pivot forces once we know the center-of-mass acceleration. We already have its position in Equation 1. Differentiating the two coordinates with respect to time gives the velocities

\[ \dot{x}_c=\dot{x}+l\dot{\theta}\cos\theta,\qquad \dot{y}_c=-l\dot{\theta}\sin\theta. \tag{4}\]

Differentiating once more, using the product and chain rules, gives the accelerations

\[ \begin{aligned} \ddot{x}_c&=\ddot{x}+l\ddot{\theta}\cos\theta-l\dot{\theta}^2\sin\theta,\\ \ddot{y}_c&=-l\ddot{\theta}\sin\theta-l\dot{\theta}^2\cos\theta. \end{aligned} \tag{5}\]

We can now substitute \(N=m\ddot{x}_c\) from Equation 3 into the cart’s equation, Equation 2, using the expression for \(\ddot{x}_c\) above. Collecting the two unknown accelerations on the left gives

\[ (M+m)\ddot{x}+ml\cos\theta\,\ddot{\theta} =F-b\dot{x}+ml\dot{\theta}^2\sin\theta. \tag{6}\]

This is the first equation of our system. It accounts for the horizontal forces on both the cart and the rod, but contains two unknown accelerations, \(\ddot{x}\) and \(\ddot{\theta}\). The second equation will come from the rod’s rotation.

Rotation of the rod

Figure 2: The solid arrows show the forces on the rod. The dashed arrows resolve the pivot forces into components perpendicular to the rod; they are not additional forces. Gravity acts through the center of mass, so its lever arm is zero.

To complete the description of the rod’s motion, we need to find how its angle changes. A force’s turning effect is its torque, measured about a chosen reference point. We’ll calculate torques about the center of mass shown in Figure 2. This lets us describe the rod’s motion as the movement of its center of mass together with a change in orientation, even while the pivot moves with the cart.

The rotational counterpart of \(F=ma\) relates the net torque \(\tau_{\mathrm{net}}\) to the angular acceleration \(\ddot{\theta}\),

\[ \tau_{\mathrm{net}}=I\ddot{\theta}. \tag{7}\]

Here, the torque and the moment of inertia \(I\) are both measured about the center of mass. The moment of inertia plays the role of mass in rotation; a larger value means the same torque produces less angular acceleration. For a thin, uniform rod, its value about the center of mass is

\[ I=\frac{mL^2}{12}. \tag{8}\]

A force’s torque has magnitude equal to the force multiplied by the perpendicular distance from the reference point to its line of action. In Figure 2, the upward force \(P\) has the horizontal lever arm \(l\sin\theta\). The rightward force \(N\) has the vertical lever arm \(l\cos\theta\). Gravity acts through the reference point itself, so its lever arm and torque about that point are zero.

We can also calculate these torques using the dashed force components perpendicular to the rod. Their magnitudes are \(P\sin\theta\) and \(N\cos\theta\); multiplying each by the distance \(l\) from the pivot to the center of mass gives the same torques.

For the orientation shown, \(P\) contributes a clockwise torque and \(N\) a counterclockwise torque about the center of mass. Counting clockwise as positive, we substitute these torques into Equation 7,

\[ I\ddot{\theta}=lP\sin\theta-lN\cos\theta. \tag{9}\]

We can remove \(N\) and \(P\) using the force equations we wrote earlier. From Equation 3, \(N=m\ddot{x}_c\) and \(P=mg+m\ddot{y}_c\). Substituting these into the torque equation gives

\[ I\ddot{\theta}=mgl\sin\theta +ml\left(\ddot{y}_c\sin\theta-\ddot{x}_c\cos\theta\right). \tag{10}\]

All that remains is to substitute the accelerations from Equation 5. In \(\ddot{y}_c\sin\theta-\ddot{x}_c\cos\theta\), the terms involving \(\dot{\theta}^2\) cancel. One is \(-l\dot{\theta}^2\cos\theta\sin\theta\) and the other is its positive counterpart. The remaining terms simplify using the Pythagorean trigonometric identity, \(\sin^2\theta+\cos^2\theta=1\),

\[ \begin{aligned} \ddot{y}_c\sin\theta-\ddot{x}_c\cos\theta &=-\ddot{x}\cos\theta-l\ddot{\theta}(\sin^2\theta+\cos^2\theta)\\ &=-\ddot{x}\cos\theta-l\ddot{\theta}. \end{aligned} \tag{11}\]

Substituting this result into Equation 10 and gathering the accelerations on the left gives our second equation of motion,

\[ ml\cos\theta\,\ddot{x}+(I+ml^2)\ddot{\theta}=mgl\sin\theta. \tag{12}\]

To see what Equation 12 predicts for the right-leaning rod in Figure 1, we can put its angular acceleration on its own,

\[ \ddot{\theta}=\frac{mgl\sin\theta-ml\ddot{x}\cos\theta}{I+ml^2}. \tag{13}\]

For a small rightward tilt, \(\sin\theta\) and \(\cos\theta\) are both positive. With the cart held still, \(\ddot{x}=0\), so gravity gives a positive angular acceleration. A rod released from rest will therefore lean farther right. Accelerating the cart right makes \(\ddot{x}\) positive, so the second term in the numerator subtracts from gravity’s contribution. If it is large enough to make the numerator negative, \(\ddot{\theta}<0\) and a rod falling right begins to slow. This is the effect we need when we move our hand underneath a leaning broom.

Matrix form

We now have two equations, Equation 6 and Equation 12, relating the cart’s acceleration to the rod’s angular acceleration. Both have to hold at the same time. Writing them together in matrix form lets us solve for \(\ddot{x}\) and \(\ddot{\theta}\) in one calculation.

\[ \begin{bmatrix} M+m & ml\cos\theta \\ ml\cos\theta & I+ml^2 \end{bmatrix} \begin{bmatrix} \ddot{x} \\ \ddot{\theta} \end{bmatrix} = \begin{bmatrix} F-b\dot{x}+ml\dot{\theta}^2\sin\theta \\ mgl\sin\theta \end{bmatrix}. \tag{14}\]

The first row is the cart’s equation, and the second row is the rod’s rotational equation. The matrix on the left contains the coefficients multiplying the two accelerations. At a given instant, the masses and length are fixed, the angle and velocities come from the current state, and we choose the applied force \(F\). That gives us all the values needed to solve for the accelerations.

We can write this as \(A\mathbf{a}=\mathbf{r}\), where \(A\) is the coefficient matrix, \(\mathbf{a}\) contains the two unknown accelerations, and \(\mathbf{r}\) is the column on the right. The function numpy.linalg.solve(A, r) calculates those accelerations for us. They tell us how the motion is changing right now. To turn them into a simulation, we need a way to advance the position and angle through time.

Euler’s method

Figure 3: Euler’s method for \(\dot{y}=3y\), starting at \(y(0)=1\). The two step sizes give different approximations to \(y=e^{3t}\). The marked derivative determines the rise of the final coarse step.

The beauty of Euler’s method is its simplicity. We can approximate a function from its rate of change and a starting value, without knowing an explicit formula for the function itself. The derivative gives us the slope at the current point. Over a short interval, we can approximate the function with a straight line following that slope. At the end of the step, we calculate a new slope and repeat.

In Figure 3, the differential equation is \(\dot{y}=3y\), so the slope is three times the current value of \(y\). Starting at \(y=1\), the slope is \(3\). A step of \(0.25\) gives \(1+3\times0.25=1.75\). The next slope is \(3\times1.75=5.25\), so the next value is \(1.75+5.25\times0.25=3.0625\). Repeating this calculation gives the four straight segments in the figure.

More generally, if \(y_n\) is our current estimate and \(\dot{y}_n\) is the derivative calculated from that estimate, an Euler step is

\[ y_{n+1}=y_n+\dot{y}_n\,\Delta t. \tag{15}\]

The straight segments in Figure 3 approximate the curved solution. They drift below it because the slope is increasing within each step, while we keep using the slope from the start. The dashed approximation uses \(\Delta t=0.05\) instead of \(0.25\). Updating the slope five times as often brings it closer to the exact curve. The exact curve is available for this example, but we don’t need it to perform the calculation.

Our cart and pendulum work the same way. Acceleration is the rate of change of velocity, and velocity is the rate of change of position. Starting from the current position, angle, and velocities, we solve Equation 14 for the accelerations. We then use those accelerations to update the velocities, and use the new velocities to update the position and angle,

\[ \begin{aligned} \dot{x}_{\mathrm{new}}&=\dot{x}+\ddot{x}\,\Delta t, & x_{\mathrm{new}}&=x+\dot{x}_{\mathrm{new}}\,\Delta t,\\ \dot{\theta}_{\mathrm{new}}&=\dot{\theta}+\ddot{\theta}\,\Delta t, & \theta_{\mathrm{new}}&=\theta+\dot{\theta}_{\mathrm{new}}\,\Delta t. \end{aligned} \tag{16}\]

Updating the velocities before the positions is called the Euler–Cromer method. The linked explanation compares this update order with using the old velocities.

We repeat the calculation with the updated state at each timestep. This lets us find where the pendulum is after, say, 12 seconds, even though we never derived a formula for its position as a function of time. For the simulation below, we’ll use \(\Delta t=0.001\ \mathrm{s}\). Halving the timestep and comparing the motion over the same interval shows how much the result depends on the step size.

Simulation

Figure 4: Unforced motion from a 0.1 radian initial tilt, using the parameters and 1 ms timestep in the code below.

With no applied force, the initial tilt grows and the rod falls, as shown in Figure 4. The cart’s position in the lower plot initially decreases, even though we aren’t applying a force to it. It then reverses direction as the rod swings around the pivot.

The angle in the upper plot passes \(180^\circ\), where the rod points straight down. By 1.2 seconds, the rod has passed the bottom and is rising on the other side. The cart moves throughout this fall because of the force transmitted through the pivot.

Starting exactly upright and motionless, with no applied force, would give zero acceleration for both the cart and rod. The ideal simulation would stay there indefinitely. Instability means that a small disturbance grows, so we give the rod the \(0.1\) radian initial tilt used above.

The code follows the same order as the calculation. A and r come from Equation 14; v is \(\dot{x}\) and omega is \(\dot{\theta}\). Both accelerations are calculated before we update either velocity.

import numpy as np

M, m, L = 1.0, 0.5, 1.0
b, g = 0.1, 9.81
l = L / 2
I = m * L**2 / 12


def step(state, force, dt):
    x, v, theta, omega = state
    coupling = m * l * np.cos(theta)

    A = np.array([
        [M + m, coupling],
        [coupling, I + m * l**2],
    ])
    r = np.array([
        force - b * v + m * l * omega**2 * np.sin(theta),
        m * g * l * np.sin(theta),
    ])
    cart_acc, angular_acc = np.linalg.solve(A, r)

    v += cart_acc * dt
    omega += angular_acc * dt
    x += v * dt
    theta += omega * dt
    return x, v, theta, omega


if __name__ == "__main__":
    state = (0.0, 0.0, 0.1, 0.0)
    for _ in range(1200):
        state = step(state, force=0.0, dt=0.001)

    x, v, theta, omega = state
    print(f"x = {x:.3f} m, theta = {np.degrees(theta):.1f} degrees")

After 1.2 seconds, the simulation prints

x = 0.176 m, theta = 251.4 degrees

Now we can try balancing it ourselves.

If you bring the rod back to upright and release the force while it is still rotating, it falls to the other side. Reaching the right angle is only part of balancing it; we also have to slow its rotation. In the next post, we’ll use a fuzzy controller to decide how to push based on both the angle and angular velocity.

Citation

Cited as:

Englert, Brunó B. (July 26, 2023). 2D Physics of the Inverted Pendulum on a Cart. englert.ai.

Or

@misc{englert2023physicscartpole,
  title   = {2D Physics of the Inverted Pendulum on a Cart},
  author  = {Englert, Brunó B.},
  journal = {englert.ai},
  year    = {2023},
  month   = {July},
  url     = {https://englert.ai/posts/001_inverted_pendulum/englert_ai_1_cart_physics.html}
}

Source code

To run the Python examples, set up an environment with venv or Miniconda. Install the required packages with

python3 -m pip install numpy arcade

• Browser simulation — the interactive version above.
• Python/Arcade simulation — the interactive desktop version.