|
| 1 | +""" |
| 2 | +Simulation of a free-cart pendulum |
| 3 | +
|
| 4 | +Parameters: |
| 5 | + L - length of the rod |
| 6 | + m - mass of the load on the top of the rod |
| 7 | + M - mass of the cart |
| 8 | +
|
| 9 | +System of equations of motion: |
| 10 | + L * th'' = -g * sin(th) + x'' * cos(th), |
| 11 | + (m + M) * x'' + m * th'' * L * cos(th) - m * L * Th'^2 * sin(th) = 0 |
| 12 | +
|
| 13 | +System: |
| 14 | + th' = Y, |
| 15 | + Y' = (g * sin(th) + b * L * Y**2 * sin(th) * cos(th)) / (L * (1 + b * cos(th)**2)), |
| 16 | + x' = Z, |
| 17 | + Z' = b * (L * Y**2 * sin(th) - g * sin(th) * cos(th)) / (1 + b * cos(th)**2), |
| 18 | + where b = m / (M + m) |
| 19 | +
|
| 20 | +State: |
| 21 | + [th, Y, x, Z] |
| 22 | +
|
| 23 | +References: |
| 24 | + - (Original example)[https://matplotlib.org/gallery/animation/double_pendulum_sgskip.html] |
| 25 | +""" |
| 26 | + |
| 27 | +import numpy as np |
| 28 | +import matplotlib.pyplot as pp |
| 29 | +import scipy.integrate as integrate |
| 30 | +import matplotlib.animation as animation |
| 31 | +from matplotlib.patches import Rectangle |
| 32 | + |
| 33 | +from math import pi |
| 34 | +from numpy import sin, cos |
| 35 | + |
| 36 | +# physical constants |
| 37 | +g = 9.8 |
| 38 | +L = 1.0 |
| 39 | +m = 0.5 |
| 40 | +M = 1.0 |
| 41 | +b = m / (m + M) |
| 42 | + |
| 43 | +# simulation time |
| 44 | +dt = 0.05 |
| 45 | +Tmax = 20 |
| 46 | +t = np.arange(0.0, Tmax, dt) |
| 47 | + |
| 48 | +# initial conditions |
| 49 | +Y = .0 # pendulum angular velocity |
| 50 | +th = pi/3 # pendulum angle |
| 51 | +x = .0 # cart position |
| 52 | +Z = .0 # cart velocity |
| 53 | + |
| 54 | +state = np.array([th, Y, x, Z]) |
| 55 | + |
| 56 | +def derivatives(state, t): |
| 57 | + ds = np.zeros_like(state) |
| 58 | + |
| 59 | + ds[0] = state[1] |
| 60 | + ds[1] = (g * sin(state[0]) + b * L * state[1]**2 * sin(state[0]) * cos(state[0])) / (L * (1 + b * cos(state[0])**2)) |
| 61 | + ds[2] = state[3] |
| 62 | + ds[3] = b * (L * state[1]**2 * sin(state[0]) - g * sin(state[0]) * cos(state[0])) / (1 + b * cos(state[0])**2) |
| 63 | + |
| 64 | + return ds |
| 65 | + |
| 66 | +print("Integrating...") |
| 67 | +# integrate your ODE using scipy.integrate. |
| 68 | +solution = integrate.odeint(derivatives, state, t) |
| 69 | +print("Done") |
| 70 | + |
| 71 | +ths = solution[:, 0] |
| 72 | +xs = solution[:, 2] |
| 73 | + |
| 74 | +pxs = L * sin(ths) + xs |
| 75 | +pys = L * cos(ths) |
| 76 | + |
| 77 | +fig = pp.figure() |
| 78 | +ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2)) |
| 79 | +ax.set_aspect('equal') |
| 80 | +ax.grid() |
| 81 | + |
| 82 | +patch = ax.add_patch(Rectangle((0, 0), 0, 0, linewidth=1, edgecolor='k', facecolor='g')) |
| 83 | + |
| 84 | +line, = ax.plot([], [], 'o-', lw=2) |
| 85 | +time_template = 'time = %.1fs' |
| 86 | +time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes) |
| 87 | + |
| 88 | +cart_width = 0.3 |
| 89 | +cart_height = 0.2 |
| 90 | + |
| 91 | +def init(): |
| 92 | + line.set_data([], []) |
| 93 | + time_text.set_text('') |
| 94 | + patch.set_xy((-cart_width/2, -cart_height/2)) |
| 95 | + patch.set_width(cart_width) |
| 96 | + patch.set_height(cart_height) |
| 97 | + return line, time_text, patch |
| 98 | + |
| 99 | + |
| 100 | +def animate(i): |
| 101 | + thisx = [xs[i], pxs[i]] |
| 102 | + thisy = [0, pys[i]] |
| 103 | + |
| 104 | + line.set_data(thisx, thisy) |
| 105 | + time_text.set_text(time_template % (i*dt)) |
| 106 | + patch.set_x(xs[i] - cart_width/2) |
| 107 | + return line, time_text, patch |
| 108 | + |
| 109 | +ani = animation.FuncAnimation(fig, animate, np.arange(1, len(solution)), |
| 110 | + interval=25, blit=True, init_func=init) |
| 111 | + |
| 112 | +pp.show() |
| 113 | + |
| 114 | +# Set up formatting for the movie files |
| 115 | +# Writer = animation.writers['ffmpeg'] |
| 116 | +# writer = Writer(fps=15, metadata=dict(artist='Sergey Royz'), bitrate=1800) |
| 117 | +# ani.save('free-cart.mp4', writer=writer) |
| 118 | + |
| 119 | + |
| 120 | + |
0 commit comments