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

Modelling a Double Pendulum in Simulink

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.

[latexpage]

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

Theory

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.\ref{fig:pendulum} can be derived by first considering the x and y components of displacement:
\begin{eqnarray}
x_1 &=& l_1\sin{\theta_1} \cr
y_1 &=& -l_1\cos{\theta_1} \cr
x_2 &=& l_1\sin{\theta_1}+l_2\sin{\theta_2} \cr
y_2 &=& -l_1\cos{\theta_1}-l_2\cos{\theta_2}
\end{eqnarray}

Leading to potential and kinetic energies of the system:
\begin{eqnarray}
\textrm{PE} &=& V = m_1gy_1 + m_2gy_2 \nonumber \cr
&=& -(m_1+m_2)gl_1\cos{\theta_1} – m_2gl_2\cos{\theta_2} \cr
\textrm{KE} &=& T = \frac{1}{2}m_1v_1^2+\frac{1}{2}m_2v_2 \nonumber \cr
&=& \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} \label{eq:kenetic}
\end{eqnarray}

And Langrangian:
\begin{equation}
L \equiv T-V \nonumber
\end{equation}
\begin{eqnarray}
= \lefteqn{\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}} \nonumber \cr && + 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}
\end{eqnarray}

For $\theta_1$, the simplified Euler-Lagrange differential equation becomes:
\begin{eqnarray} \label{eq:theta1}
\lefteqn{(m_1+m_2)l_1\ddot{\theta_1} + m_2l_1\ddot{\theta_2}\cos{(\theta_1-\theta_2)}} \nonumber \cr && + m_2l_2\dot{\theta_2^2}\sin{(\theta_1-\theta_2)} + (m_1+m_2)g\sin{\theta_1} = 0
\end{eqnarray}

Or in this case $LHS=T$. Similarly, $\theta_2$ becomes:
\begin{eqnarray} \label{eq:theta2}
\lefteqn{m_2l_2\ddot{\theta_2} + m_2l_1\ddot{\theta_1}\cos{(\theta_1-\theta_2)}} \nonumber \cr && – m_2l_1\dot{\theta_1^2}\sin{(\theta_1-\theta_2)} + m_2g\sin{\theta_2} = 0
\end{eqnarray}

Simulink Model

The Simulink model was constructed by first considering the various terms which make up the equations of motion (\ref{eq:theta1}), (\ref{eq:theta2}); 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.\ref{fig:mathblock}].

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{eqnarray}
\ddot{\theta_1} = \frac{b_1 – b_2}{b_3} \cr
\ddot{\theta_2} = \frac{b_2 – b_4}{b_5}
\end{eqnarray}

Where:

\begin{eqnarray}
b_1 &=& T – (m_1+m_2)gsin{\theta_1} – m_2l_2\dot{\theta_2^2}sin{(\theta_1-\theta_2)} \cr
b_2 &=& m_2l_2\dot{\theta_1}^2sin{(\theta_1-\theta_2)}-m_2gsin{\theta_2} \cr
b_3 &=& (m_1+m_2)l_1 – m_2l_1\cos^2(\theta_1-\theta_2) \cr
b_4 &=& m_2l_1\ddot{\theta_1}cos(\theta_1-\theta_2) \footnotemark \cr
b_5 &=& m_2l_2
\end{eqnarray}

[Fig.\ref{fig:subsystem}] 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, referring to (\ref{eq:kenetic}). 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: \ref{sec:simuapp}]. Firstly, for Steps.2-5 of the experiment, the system of Fig.\ref{fig:dbpmain} 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.\ref{fig:dbpfbk}].

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:

Pendulum

cropped-JBRLogo.png

15 thoughts on “Modelling a Double Pendulum in Simulink”

    1. Hello! I’m fairly sure it is correct. If you go through the derivation of [\ddot{\theta_1}] (I just did – admittedly not all the way as it gets pretty lengthy!) subbing in eq (6) for [\ddot{\theta_2}], you end up with an extra [m_2l_1\cos^2(\theta_1-\theta_2)], which you can further reduce terms to provide the [m_2g\sin(\theta_2)] (in b2). I’m a little rusty though so please share if you see something different.

  1. Hello 🙂 Im not too sure about the exacts, I reckon that when you substitute for thetaddot2 in equation 6 you have two other terms sin terms the one with thetadot squared and the other with just mgsin(theta2). Haha hope im not making a fool of myself by missing something very obvious.

  2. Please elaborate how u entered the b1, b2, b3….b5 in fcn function block. Really need your help. I m doing the same work.

    1. Each variable required in the b1,b2,b3..b5 functions have be multiplexed onto on input line (the black junction boxes). To enter the function into the function block, you have to reference them using ‘u(n)’, where n is the number of the input (count down from the top at 1).

      For example. rotational velocity (dtheta^2) would be u(1). Torque (T) is u(7).

      The final block function (Eq. b5) also gives this away as it’s name is explicit:

      u(9)*u(10) [m2*l2]

  3. Hi John, i tried simulate in my matlab but i found some errors. What version you use to simulate this?
    And thx for share this work with us.

      1. I fixed it, was a bug on my pc. I have a question, i need use a torsional spring in the base of pendulum, i tried using anothers blocks in simmetechs library, but is failed. The equation of this spring is given by: M=k.theta=k(theta1-theta2) and U=0.5*(k*theta^2) (http://goo.gl/YVVPYQ) page 87. I tried modeling with the first equation but i failed, can you help me?

  4. Hi,

    I’m just worried about how the simulation gets 2xn matrix for t1 ? Then you take the first row as t and second row as t1 but why is that happening ? What is “t” there ? If time, then the values of t are separated according to what ?

    To be more clear, I’m trying to do a sim by using not that T(Torque) but using a sine wave. And the thing is as I increase the amplitude of the sine wave the number of values for t t1 and t2 increasing so the simulation becomes very hard to solve. Becoming hard is not a problem but I’m confused about how does it determine the values of time there ?

    Thanks.

  5. Could you please share the equations of theta1 and theta2? I tried to do modeling the same way what you did? But mine was little bit wrong! Please help

Leave a Reply