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]

**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}].

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 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.

Here’s the code and model:

good job

Hey, there might be an error in your calculation.. Im slightly unsure if b2 and b3 are correct

Hello! I’m fairly sure it is correct. If you go through the derivation of (I just did – admittedly not all the way as it gets pretty lengthy!) subbing in eq (6) for , you end up with an extra , which you can further reduce terms to provide the (in b2). I’m a little rusty though so please share if you see something different.

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.

I have the same doubt as Narad

Hi, good job

I have a question about the PID control system. How do you determine the coefficients of the controller?

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

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]

Hi John,

Can this simulink model be simulated as a 2-link arm segment model?

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.

Hi Fellipe,

It was a while ago so I was probably on 2011b. What errors do you get?

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?

SimMechanics Library*

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.

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