Modelling a Double Pendulum in Simulink

April 26, 2013

One module I took during the final year of my degree was ‘System Modelling and Simulation’. A well taught and great module, one of the tasks was to model a double pendulum.

The approach involved deriving the equations of from the highest order of motion for each mass then working backwards through Simulink blocks to generate each term, which could then be to solve the equation – a bit of a chicken and egg problem! It was a an excellent task as the idea seems a little backwards at first and gave me a fresh approach to problem solving a model.

Below is the method extract from my report as the YouTube demo has generated some interest in the solution.

The double pendulum to be modelled.
The double pendulum to be modelled.


Double pendula are an example of a simple physical system which can exhibit chaotic behavior [1]. That is to say it is sensitive to initial conditions and that its movements, although predictable, appear random. Equations of motion for the pendulum of Fig.?? can be derived by first considering the x and y components of displacement:

\[\begin{align} x_1 = l_{1}\sin(\theta_{1}) \\ y_2 = -l_{1}\cos(\theta_{1}) \\ x_2 = l_{1}\sin(\theta_{1}) + l_{2}\sin(\theta{2}) \\ y_2 = -l_{1}\cos(\theta_{1}) - l_{2}\cos(\theta{2}) \end{align}\]

Leading to potential and kinetic energies of the system:

\[\begin{align} PE = V = m_1gy_1 + m_2gy_2 \\ = -(m_1+m_2)gl_1\cos{\theta_1} - m_2gl_2\cos{\theta_2} \\ KE = T = \frac{1}{2}m_1v_1^2+\frac{1}{2}m_2v_2 \\ = \underbrace{\frac{1}{2}m_1l_1^2\dot{\theta_1^2}}{m_1} + \underbrace{\frac{1} {2}m_2[l_1^2\dot{\theta_1^2}+l_2^2\dot{\theta_2^2}+2l_1l_2\dot{\theta_1}\dot{\theta_2}\cos{(\theta_1 - \theta_2)}]}{m_2} \end{align}\]

And Langrangian:

\[L \equiv T-V \\ = \frac{1}{2}(m_1+m_2)l_1^2\dot{\theta_1^2} + \frac{1}{2}m_2l_2^2\dot{\theta_2^2} \\ + m_2l_1l_2\dot{\theta_1}\dot{\theta_2}\cos{(\theta_1-\theta_2)} + (m_1+m_2)gl_1\cos{\theta_1} + m_2gl_2cos{\theta_2}\]

For \(\theta_1\), the simplified Euler-Lagrange differential equation becomes:

\[(m_1+m_2)l_1\ddot{\theta_1} + m_2l_1\ddot{\theta_2}\cos{(\theta_1-\theta_2)} + m_2l_2\dot{\theta_2^2}\sin{(\theta_1-\theta_2)} + (m_1+m_2)g\sin{\theta_1} = 0\]

Or in this case \(LHS=T\). Similarly, \(\theta_2\) becomes:

\[m_2l_2\ddot{\theta_2} + m_2l_1\ddot{\theta_1}\cos{(\theta_1-\theta_2)} - m_2l_1\dot{\theta_1^2}\sin{(\theta_1-\theta_2)} + m_2g\sin{\theta_2} = 0\]

Simulink Model

The Simulink model was constructed by first considering the various terms which make up the equations of motion (5), (6); each of the terms can be derived from the highest order of \(\theta\): \(\ddot{\theta}\). Hence, development begin by creating a mathematical flow of blocks, with each output a term within the equation [Fig.??].

The starting point for the Simulink model.
The starting point for the Simulink model.

Outputs for \(\dot{\theta_1^2}, \dot{\theta_2^2}, \sin{\theta_1}\), etc. had been created but the inputs \(\ddot{\theta_1}\) and \(\ddot{\theta_2}\) where unknown. By rearranging the equations of motion to a set of functions in terms of \(\ddot{\theta_1}\) and \(\ddot{\theta_2}\), they can be applied to the outputs and result fed back to the inputs:

\[\begin{align} \ddot{\theta_1} = \frac{b_1 - b_2}{b_3} \\ \ddot{\theta_2} = \frac{b_2 - b_4}{b_5} \end{align}\]


\[\begin{align} b_1 = T - (m_1+m_2)gsin{\theta_1} - m_2l_2\dot{\theta_2^2}sin{(\theta_1-\theta_2)} \\ b_2 = m_2l_2\dot{\theta_1}^2sin{(\theta_1-\theta_2)}-m_2gsin{\theta_2} \\ b_3 = (m_1+m_2)l_1 - m_2l_1\cos^2(\theta_1-\theta_2) \\ b_4 = m_2l_1\ddot{\theta_1}cos(\theta_1-\theta_2) \\ b_5 = m_2l_2 \end{align}\]

The outputs are combined using a Mux block, a multiplexor which creates an array of its inputs. Similarly, the masked system parameters are defined using MATLAB variables as Constant blocks and combined with another Mux. The \(b\) terms are then found by creating Function blocks referencing these arrays, which are then equated accordingly to find \(\ddot{\theta_1}\) and \(\ddot{\theta_2}\). Kinetic energies of the masses are also found at this stage in a similar manor.. The angular accelerations are then fed back to the start of the block diagram. Initial angle conditions are set in the displacement Integrator to the MATLAB variables.

The complete un-masked sub-system
The complete un-masked sub-system

The pendulum model was ‘masked’ into a sub-system block for cleaner integration with the proceeding models [Appendix: ??]. Firstly, for Steps.2-5 of the experiment, the system of Fig.?? was used, which saves the outputs of the pendulum to file for manipulation within MATLAB. Secondly, for Step.5, a closed loop system with a PID block was used [Fig.??].

For the purpose of experimenting with the system for the report, I created two simple MATLAB scripts:

  • Pre-processing – to define the variables (masses, length, gravity etc.)
  • Post-processing – to create output plots and an animation.

This made the whole model very easy to use and play around with. I also various disturbances and feedback loops to the system, that can be turned on and off.

The final system with masked sub-system, PID control and disturbances.
The final system with masked sub-system, PID control and disturbances.

Here’s the code and model: