Pendulum

System of differential equations for a pendulum is derived using Newtonian and Lagrangian mechanics. The system is then linearized about three different equilibium points. Transfer function and state space model are developed for the linearized systems. Simulation code is also provided.

System

The pendulum is made of a point mass attached to the end of a massless bar. The other end of the bar is attached to a fixed ground. The mass and the bar rotate about this fixed point. Gravity acts downwards.

The point mass has mass m and the bar is of length l. The angle the bar makes with the downward direction is θ. When the mass is at the bottom-most position, θ = 0 and θ increases as the pendulum rotates counter-clockwise. There is also air resistance which acts in the direction opposite to the velocity and whose magnitude scales with the speed by μ. An external torque τ is applied to this pendulum. τ is positive counter-clockwise.

Differential Equation (Newtonian Mechanics)

Differential equation describing the dynamics of the pendulum is derived using Newtonian mechanics. Since there is only one unknown, one equation is needed.

Torques around the axis of rotation can be visualized using the following free body diagram:

where I is the inertia and is ml² for a point mass rotating l away from the center of rotation. α is angular acceleration.

Solving for angular acceleration, the differential equation describing the motion of the pendulum can be written as:

Differential Equation (Lagrangian Mechanics)

Differential equation describing the dynamics of the pendlum is derived using Lagrangian mechanics.

where the height h is 0 at the axis of rotation and positive upwards so h = -lcosθ.

Velocity-proportional friction forces can be accounted in Lagrangian mechanics using Rayleigh dissipation function:

and the equation of motion with external force is solved with Lagrange's equation,

from sympy import *

t, m, g, l, mu, tau = symbols('t m g l mu tau')
theta = symbols(r'\theta', cls=Function)
theta = theta(t)
theta_d = diff(theta, t)
theta_dd = diff(theta_d, t)

v = l * theta_d
h = -l * cos(theta)

T = 1/2 * m * v ** 2
V = m * g * h
L = T - V
D = 1/2 * mu * theta_d ** 2

# Lagrangian
LE = diff(diff(L, theta_d), t) - diff(L, theta) + diff(D, theta_d) - tau

# Set equal to zero and solve for theta_dd
rhs = solve(LE, theta_dd)

pprint(rhs)

State Space Model

The equation of motion derived in the previous sections is a second order differential equation. This is converted to a system of first order differential equations by letting

Linearizing at Equilibrium Points

This refers to the state when the pendulum is not moving nor accelerating. Letting

Doing Taylor expansion about vecδx=0 and δu=0 and ignoring second and higher order terms, the nonlinear state space model can be linearized as

Equilibrium 1: when the pendulum is straight down

This occurs when θ = 0 and τ = 0 (when no external torque is applied). Letting

The eigenvalues of the system matrix A provides insights into the stability of the system. Analytical solution is given as:

If we let g = 9.81, m = 0.5, l = 0.2, μ = 0.15, the eigenvalues of the system matrix A are -3.75 ± 5.91502i. Since all real parts are negative, the system is stable. A phase plane also shows that this is a stable node:

Next, it is possible to compute the transfer function using the following formula:

which is a simple harmonic oscillator if there is no friction (μ = 0) and no external torque (τ = 0).

from sympy import *

m, g, l, mu, u, s, t = symbols('m g l mu u s t')
x = MatrixSymbol('x', 2, 1)

F = Matrix([
    [x[1]],
    [-g/l * sin(x[0]) - mu/(m * l**2) * x[1] + 1/(m * l**2) * u]
])

G = Matrix([
    [x[0]]
])

x_eq = Matrix([
    [0], # 0 or pi
    [0]
])

u_eq = 0

A = F.jacobian(x).subs(x, x_eq)
B = F.jacobian([u]).subs(u, u_eq)
C = G.jacobian(x).subs(x, x_eq)
D = G.jacobian([u]).subs(u, u_eq)

transfer_function = C * (s * eye(2) - A) ** -1 * B + D

A_subs = A.subs([(g, 9.81), (m, 0.5), (l, 0.2), (mu, 0.15)])

pprint(A)
pprint(B)
pprint(C)
pprint(D)
pprint(A.eigenvals())
pprint(A_subs)
pprint(A_subs.eigenvals())
pprint(transfer_function)

import numpy as np
import matplotlib.pyplot as plt

v = 3
w = np.pi/4
Y, X = np.mgrid[-v:v:100j, -w:w:100j]

A = np.array([
    [0, 1],
    [-49.05, -7.5]
])

s, v = np.linalg.eig(A)

U = A[0,0]*X + A[0,1]*Y
V = A[1,0]*X + A[1,1]*Y

if s.imag[0] == 0 and s.imag[1] == 0:
    t = np.arange(-w, w, 0.01)
    plt.plot(t, (v[1,0]/v[0,0])*t)
    plt.plot(t, (v[1,1]/v[0,1])*t)
plt.streamplot(X, Y, U, V)
plt.xlabel(r'$\delta \theta$')
plt.ylabel(r'$\delta \dot{\theta}$')
plt.show()

Equilibrium 2: when the pendulum is straight up

This occurs when θ = π and τ = 0 (when no external torque is applied). Letting

If we let g = 9.81, m = 0.5, l = 0.2, μ = 0.15, the eigenvalues of A are -11.69434 and 4.19434. Since there exists an eigenvalue with positive real part, it is unstable. There is one positive and one negative eigenvalue so this is a saddle point and this can also be seen from the phase portrait:

Equilibrium 3: when the pendulum is tilted 45 degrees

The pendulum can be in equilibium in a tilted position if there is an external torque holding it in place. This happens at:

The external torque that just balances the torque produced by the gravity can be obtained using a free body diagram:

Another way to obtain the torque at this equilibium is by solving it algebraically:

It is also possible to find this equilibium point by using numerical method. If we let g = 9.81, m = 0.5, l = 0.2, μ = 0.15, an optimization solver is used to find the equilibium torque:

from scipy import optimize
import numpy as np

# Constants
g = 9.81
l = 0.2
m = 0.5
mu = 0.15

# Desired angle at equilibrium
theta_eq = np.pi/4

# Nonlinear pendulum differential equation: xdot = f(x, u)
def f(x, u):
    return [
        x[1],
        -g / l * np.sin(x[0]) - mu/(m * l**2) * x[1] + 1/(m * l**2) * u
    ]

# Optimization function
def fun(z):
    x = z[:2]
    u = z[2]
    xdot = f(x, u)
    theta_err = x[0] - theta_eq
    return [
        xdot[0],
        xdot[1],
        theta_err
    ]

# Initial guess
x0 = np.array([theta_eq, 0])
u0 = np.array([0])

# Find equilibrium
z0 = np.concatenate((x0, u0))
result = optimize.root(fun, z0)

print(result.x[:2]) # x at equilibium
print(result.x[2]) # u at equilibrium

Above code returns the equilibium torque to be 0.69367 which is identical to the analytical solution:

Simulation

The linearized dynamics is only valid for small angles around the equilibium point. The full nonlinear dynamics can be simulated using numerical integration. Python's ODE solver is used:

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Constants
g = 9.81
l = 0.2
m = 0.5
mu = 0.02

def xdot(t, x):
    u = 0 # No input torque
    return [
        x[1],
        -g / l * np.sin(x[0]) - mu/(m * l**2) * x[1] + 1/(m * l**2) * u
    ]

x0 = [0, 17]  # Initial state [theta0, theta_dot0]
t0 = 0
tf = 6
dt = 0.01

# Simulate
sol = solve_ivp(xdot, [t0, tf], x0, t_eval=np.arange(t0, tf, dt))

# Plot as function of time
fig, axs = plt.subplots(2)
axs[0].plot(sol.t, sol.y[0])
axs[0].set_ylabel(r'$\theta$')
axs[0].grid()
axs[1].plot(sol.t, sol.y[1])
axs[1].set_ylabel(r'$\dot{\theta}$')
axs[1].grid()
plt.show()

In this simulation, the pendulum is initially positioned at the bottom with a counter-clockwise angular speed. The pendulum passes through the top and comes back down. After oscillating a little, it slowly comes to rest due to air resistance. There is no input torque. The animation at the top of this page was created with this method.

(input)	τ(t)	=	external torque applied to pendulum
(output)	θ(t)	=	angle of pendulum	unknown
	m	=	mass
	l	=	arm length
	μ	=	air friction coefficient