Derivation of Rod Pendulum (Python Animation)

This is the part 1 in a series of blog posts devoted to the classic control problem of an inverted pendulum on a driven cart. This series is still in production (Nov. 2020) but when finished they will be structured as follows.

  1. Derivation of Simple Pendulum (Python Simulation) - link
  2. Building a Physical Inverted Pendulum
  3. Controlling Inverted Pendulum with PID controller
  4. State Space Model of Inverted Pendulum - link
  5. Designing the Optimal Controller for an Inverted Pendulum (LQR control) - link
  6. Simulating an Inverted Pendulum with an LQR Controller (Python) - link
  7. Deploying an LQR Controller on a Physical Inverted Pendulum

The associated code can be found on GitHub and a video series on YouTube (not currently released as of Nov. 2020). While this project aims to be encompassing and self contained, it does assume familiarity with classical mechanics, linear algebra, differential equations, some control theory, and programming (Python).

Animation of simple pendulum

This blog is supplemental documentation pertaining to the derivation for the equation of motion used in the python program,

While the program uses the scipy.integrate.solve_ivp() function, I will also touch on the small angle approximation for solving the system analytically for small initial angles $\theta_0 << 1$. There are a few different ways to derive the equations of motion, I will be using the torque derivation, and mention in passing the Lagrangian approach.

Firstly, since this is a rod-pendulum with uniform mass (unlike the slightly more common, point mass attached to a massless rod pendulum), the moment of inertia $I$ must be found.

$$I_{rod,CoM} = \int_{0}^M r^2dm = \int_{-\frac{L}{2}}^{\frac{L}{2}}\frac{r^2m}{L}dr = \frac{1}{12}mL^2 $$

Using the parallel axis theorem, we can shift the moment of inertia from the center of mass (CoM) to one of the ends:

$$I_{rod,end} = I_{rod,CoM} + md^2 = \frac{1}{3}mL^2$$

We now have the information required to find the equation of motion via the torque:

$$ \tau = I\alpha $$
$$ -mg\sin(\theta)(\frac{L}{2}) = I_{r,e}\frac{d^2 \theta}{d t^2}$$
$$\frac{d^2 \theta}{d t^2} = - \omega^2 \sin(\theta)$$

where $\omega = \sqrt{\frac{3g}{2L}}$ is the natural frequency of the pendulum. At this point if using a small initial angle $\theta_0 << 1$, one could solve the system analytically using the small angle approximation $\sin(\theta) \approx \theta$, and then integrating twice to give $\theta(t) = \theta_0 \cos(\omega t)$. Refer to the first couple of commits in the repository to see the code for a small angle approx (of a point-mass pendulum). At this point the pendulum animation can be created using matplotlib. Note: this blog post focuses on the physics of the system, not the actual software. For a deeper understanding of how scipy.integrate.solve_ivp() is used please refer to the scipy documentation.

To create the phase space plot I will calculate the Hamiltonian and Lagrangian (optional) of the system. The kinetic and potential energy can be found as:

$$T = \frac{1}{2}mv^2 + \frac{1}{2}I\dot\theta^2 = \frac{1}{6}mL^2\dot\theta^2$$

$$V = \frac{mgL}{2}(1-\cos(\theta))$$

The Lagrangian, $L=T-V$, can be used with the Euler-Lagrange equation to find the equation of motion (as opposed to the torque derivation, above). This is left as an exercise for the reader. Using the generalised coordinate, $q=\theta$, and generalized momentum, $p = \frac{\partial L}{\partial \dot\theta}$, the Hamiltonian can be expressed as:

$$H = T + V = \frac{p^2}{6m L^2} + mgL(1-\cos(q))$$

The Hamiltonian represents the total energy of the system and the different energy level contours are expressed in the phase space plot.

Show Comments