Pendulum

In this example we have a look at the mathematical pendulum.

You can also find this example as a single file in the GitHub repository.

Above we have both the mechanical mode, which is a weight dangling from a string with length l and diverted by some initial angle, as well as the block diagram representation of the equation of motion to the right.

The equation of motion of the system is a nonlinear second order ODE:

\[\ddot{\phi} = - \frac{g}{l} \sin(\phi)\]

Lets transition to Python by importing the components we need to model the system from PathSim

import numpy as np
import matplotlib.pyplot as plt

# Apply PathSim docs matplotlib style for consistent, theme-friendly figures
plt.style.use('../pathsim_docs.mplstyle')

from pathsim import Simulation, Connection

from pathsim.blocks import (Integrator,
    Amplifier, Function, Adder, Scope)

# Using an adaptive runge-kutta method
from pathsim.solvers import RKCK54

and define the system parameters. Here we choose an initial angle that is close to the pendulum beeing at the top:

# Initial angle and angular velocity
phi0, omega0 = 0.9*np.pi, 0

# Parameters (gravity, length)
g, l = 9.81, 1

We can directly translate the block diagram above to PathSim blocks and connections:

# Blocks that define the system
In1 = Integrator(omega0)
In2 = Integrator(phi0)
Amp = Amplifier(-g/l)
Fnc = Function(np.sin)  # <- function blocks just take callables
Sco = Scope(labels=["angular velocity", "angle"])

blocks = [In1, In2, Amp, Fnc, Sco]

# Connections between the blocks
connections = [
    Connection(In1, In2, Sco[0]),
    Connection(In2, Fnc, Sco[1]),
    Connection(Fnc, Amp),
    Connection(Amp, In1)
]

The simulation is initialized with the blocks and connections. In this case we dont use the default solver but an adaptive integrator RKCK54 to ensure accuracy. Its an adaptive runge-kutta method from Cash and Karp, similar to Matlabs ode45, which is from Dormand and Prince and RKDP54 in PathSim. The tolerances we set here, are also for the integrator. The adaptive method controls the timestep such that the local truncation error (lte) stays below the set tolerances.

# Simulation instance from the blocks and connections
Sim = Simulation(
    blocks,
    connections,
    dt=0.1,
    Solver=RKCK54,
    tolerance_lte_rel=1e-6,
    tolerance_lte_abs=1e-8
)

10:59:13 - INFO - LOGGING (log: True)
10:59:13 - INFO - BLOCKS (total: 5, dynamic: 2, static: 3, eventful: 0)
10:59:13 - INFO - GRAPH (nodes: 5, edges: 6, alg. depth: 3, loop depth: 0, runtime: 0.075ms)

Finally we can run the simulation for some number of seconds and see what happens.

# Run the simulation for 15 seconds
Sim.run(duration=15)

# Plot the results directly from the scope
fig, ax = Sco.plot(".-")
plt.show()

10:59:13 - INFO - STARTING -> TRANSIENT (Duration: 15.00s)
10:59:13 - INFO - --------------------   1% | 0.0s<0.2s | 1334.3 it/s
10:59:13 - INFO - ####----------------  20% | 0.0s<0.1s | 2599.2 it/s
10:59:13 - INFO - ########------------  40% | 0.1s<0.0s | 2589.4 it/s
10:59:13 - INFO - ############--------  60% | 0.1s<0.0s | 2635.4 it/s
10:59:13 - INFO - ################----  80% | 0.1s<0.0s | 1168.9 it/s
10:59:13 - INFO - #################### 100% | 0.1s<--:-- | 2638.0 it/s
10:59:13 - INFO - FINISHED -> TRANSIENT (total steps: 207, successful: 174, runtime: 120.92 ms)

Since we are using an adaptive integrator, it might be interesting to also look at the timesteps the simulation takes. To do this, we just get the times from the scope and compute the differences:

# Read the recordings from the scope
time, *_ = Sco.read()

fig, ax = plt.subplots(figsize=(8,4), tight_layout=True, dpi=120)
ax.plot(time[:-1], np.diff(time), lw=2)
ax.set_ylabel("dt [s]")
ax.set_xlabel("time [s]")
ax.grid(True)
plt.show()

We can clearly see that the integrator takes smaller steps when the pendulum gets closer to regions where the solution trajectory is more dynamic to keep the local truncation error below the tolerances.