The Inverted Pendulum - Part 1

An inverted pendulum is simply a pendulum which has its centre of mass above its pivet point. The dynamic stabilisation of the inverted pendulum is a classic problem in control theory and is used to benchmark control strategies (see PID Tutorial). The T-Bot is an example of an inverted pendulum, as is SpaceX's reusable launch system.

You might imagine, if you stabilise an inverted pendulum well enough it should stay balanced if undisturbed (maybe not in the case of a rocket!). You can test that theory using your T-Bot, a stop watch and the python code below.

Turn the T-Bot off and balance it as well as you can. Start the clock as soon as you let go. Stop it when you hear a thud. Alternatively, you can video the T-Bot falling, and use Camera Stopwatch Lite to measure the falling time. You can run the Python code below and use the fallingtime() function to work out what your starting displacement was.

Where α is the angular acceleration, ω is the angular velocity and θ is the angle of the inverted pendulum with respect to the vertical axis and 𝑙 is the length of the inverted pendulum.

Import numpy and pyplot.

import numpy as np
import matplotlib.pyplot as plt
plt.ion()

Setup the functions:

################     Functions    ####################

def fallingtime(l,g,theta,theta2,dt):
	'''
	Calculates angular acceleration, angular velocity, velocity,
	and falling time for an inverted pendulum between two angles(in degrees)
	for a given length. Usage fallingtime(l,g,theta,theta2,dt)
	Returns angular acceleration, angular velocity, velocity, time.'''
	theta = theta*np.pi/180
	theta2 = theta2*np.pi/180
	omega = 0
	t=0
	while theta < theta2:
        # angular acceleration
		alpha =  np.sin(theta)*g/l
        # integrate angular acceleration to get angular velocity
		omega += alpha*dt
        # integrate angular velocity to get angle
		theta += omega*dt 
		# integrate dt to get time
		t += dt

    # return angular acceleration, angular velocity, velocity, time
	return(alpha,omega,l*omega,t)      

def v(l,g,theta):
    if theta != 0:
        return np.sqrt(2*g*(l-l*np.cos(theta*np.pi/180)))*theta/np.abs(theta)
    else:
        return 0.0


def falling(l,g,theta,theta2,dt):
	'''
	Calculates angular acceleration, angular velocity, velocity, and falling
	time for an inverted pendulum between two angles(in degrees) for a given
	length. Usage fallingtime(l,g,theta,theta2,dt) Returns angular acceleration,
	angular velocity, velocity, theta, and time up to each step in an array'''
       
	alpha_omega_v = np.array([[]]*5).T
	theta = theta*np.pi/180
	theta2 = theta2*np.pi/180
	omega = 0
	t=0
	while theta < theta2:
        # angular acceleration
		alpha =  np.sin(theta)*g/l
        # integrate angular acceleration to get angular velocity
		omega += alpha*dt 
        # integrate angular velocity to get angle
		theta += omega*dt 
        # integrate dt to get time
		t += dt                         
		alpha_omega_v = np.vstack((alpha_omega_v,[alpha, omega, omega*l,theta*180/np.pi,t]))                     
	return(alpha_omega_v)       # angular acceleration, angular velocity, velocity, time


def v(l,g,theta):
    if theta != 0:
        return np.sqrt(2*g*(l-l*np.cos(theta*np.pi/180)))*theta/np.abs(theta)
    else:
        return 0.0

Create the Data:

   
xdata = list(np.linspace(1.E-11, 50.e-7, 500))+list(np.linspace(60.E-7, 15.e-5,500))
ydata_earth = [fallingtime(0.08, g_earth,np.arcsin(x/8e-2)*180/np.pi, 90, 0.001)[3] for x in xdata]
ydata_moon = [fallingtime(0.08, g_moon,np.arcsin(x/8e-2)*180/np.pi, 90, 0.001)[3] for x in xdata]

Now we can plot the data:

plt.figure(figsize=(10, 4))

plt.title('Inverted Pendulum')
plt.plot(xdata, ydata_moon, c=(50/255.,50/255.,250/255.),linewidth = 2, label = 'Moon g = '+str(g_moon)+' ms^2')
plt.plot(xdata, ydata_earth,  c=(255/255.,10/255.,10/255.),linewidth = 2, label = 'Earth g = '+str(g_earth)+' ms^2')

plt.fill_between(xdata, 0, ydata_earth,facecolor=(100/255.,0/255.,0/255.),edgecolor=(191/255.,211/255.,255/255.), alpha = 0.6)
plt.fill_between(xdata, 0, ydata_moon,facecolor=(0/255.,0/255.,100/255.),edgecolor=(191/255.,211/255.,255/255.), alpha = 0.6)

plt.annotate('Atomic vibrational amplitude (10 pm)\nFalling time 2.17 s',xy=(1.E-11,2.16),xytext = (1.e-5,5), arrowprops = dict(facecolor=(0/255.,255/255.,63/255.),width=1,headwidth=7,shrink=0.01),)

plt.annotate('Width of H atom (100 pm)\nFalling time 1.96 s',xy=(1.E-10,2),xytext = (1.e-5,2.5), arrowprops = dict(facecolor='white',width=1,headwidth=7,shrink=0.01),)

plt.annotate('Average width of human hair (100 microns)\nFalling time 0.712 s',xy=(100.E-6,0.712),xytext = (60.e-6,3), arrowprops = dict(facecolor='orange',width=1,headwidth=7,shrink=0.01),)

plt.legend(loc = 'best',prop={ 'size': 8})
plt.xlabel('Initial horizontal displacement (m)')
plt.ylabel('Falling Time (s)')
#plt.axis('tight')
plt.xlim((xdata[0]-xdata[100]),(xdata[-1]+xdata[100]))
plt.subplots_adjust(bottom=0.15)
plt.savefig('Falling.svg')

Finally, we can save the figure:

plt.savefig('Falling.svg')