Activity 3: Modeling of a Simple Pendulum

Key Topics: Modeling Rotational Mechanical Systems, Nonlinear Systems, Underdamped Second-Order Systems, System Identification

Equipment needed
Purpose
Modeling from first principles
Linearization
System identification experiment
Model validation
Extensions

Equipment needed

Arduino board (e.g. Uno, Mega 2560, etc.)
simple pendulum (slender metal bar with end weight) with clamp or stand
rotary potentiometer (e.g. 10K-Ohm linear taper potentiometer)

The orientation of the simple pendulum will be measured employing a rotary potentiometer. The Arduino board is simply employed for data acquisition (and to supply excitation for the potentiometer). Specifically, an Analog Input on the Arduino board is employed to read the potentiometer output which is then fed to Simulink for visualization and for comparison to our resulting simulation model output.

Purpose

The purpose of this activity with the simple pendulum system is to demonstrate how to model a rotational mechanical system. Specifically, the theory of modeling is discussed with an emphasis on which simplifying assumptions are appropriate in this case. The associated experiment is employed to demonstrate how to identify different aspects of a physical system, as well as to demonstrate the accuracy of the resulting model.

Modeling from first principles

First we will employ our understanding of the underlying physics of the simple pendulum system to derive the structure of the system model. We will term this process, "modeling from first principles." In this example we employ the following variables.

(m)       mass of the pendulum bar

(M)       mass of the pendulum end weight

(l)       length to end weight center of mass

(theta)   pendulum angle from vertical (down)

To begin, we first draw the free-body diagram where the forces acting on the pendulum are its weight and the reaction at the rotational joint. We also include a moment due to the friction in the joint (and the rotary potentiometer). The simplest approach to modeling assumes the mass of the bar is negligible and that the entire mass of the pendulum is concentrated at the center of the end weight.

The equation of motion of the pendulum can then be derived by summing the moments. We will choose to sum the moments about the attachment point $O$ since that point is the point being rotated about and since the reaction force does not impart a moment about that point.

(1) $\Sigma M = (M+m)gl\sin\theta - T_{fric} = I_O\ddot{\theta}$

Assuming that the mass of the pendulum is concentrated at its end mass, the mass moment of inertia is $I_O \approx (M+m)l^2$ . A more accurate approach would be to consider the rod and end mass explicitly. In that case, the weight of the system could be considered to be located at the system's mass center $l_G = (Ml+0.5ml))/(M+m)$ . In that case, the mass moment of inertia is $I_O = ml^2/3 + Ml^2$ . Depending on the parameters of your particular pendulum, you can assess if this added fidelity is necessary.

For the experiment we will perform shortly, the simple pendulum employed consists of a rod of length $l = 0.43\ \mbox{m}$ and mass $m = 0.095\ \mbox{kg}$ with an end mass of $0.380\ \mbox{kg}$ . Therefore, the difference between $l = 0.43\ \mbox{m}$ and $l_G = 0.39\ \mbox{m}$ is significant enough to include. The difference between $I_O \approx (M+m)l^2 = 0.088\ \mbox{kg-m}^2$ and $I_O = ml^2/3 + Ml^2 = 0.076\ \mbox{kg-m}^2$ is also significant enough to include.

We will also initially assume a viscous model of friction, that is, $T_{Fric} = -b\dot{\theta}$ where $b$ is a constant. Such a model is nice because it is linear. We will assess the appropriateness of this model later. Sometimes the frictional moment is not linearly proportional to the angular velocity. Sometimes, the stiction in the joint is significant enough that it must be modeled too.

Taking into account the above assumptions, our equation of motion becomes the following.

(2) $-(M+m)gl_G\sin\theta - b\dot{\theta} = I_O\ddot{\theta}$

Linearization

Some parameters of the given pendulum system are relatively easy to measure, for example, things such as the pendulum length or the pendulum mass. Other parameters, such as the viscous coefficient of friction $b$ , are not as easy to measure directly. Therefore, we will perform a simple experiment to help identify $b$ . This experiment will also help to validate some of the simplifying assumptions we made in the process of generating our model.

The fact that the system model we generated is nonlinear makes the parameter identification process a little more challenging. However, we can use a linearized version of the model to help us in the identification process. Presuming that for our experiment the pendulum swings through small angles (about $\theta = 0$ ), we can use the approximation that $\sin \theta \approx \theta$ . Therefore, our linearized model becomes the following.

(3) $\ddot{\theta} + \frac{b}{I_O}\dot{\theta} + \frac{(M+m)gl_G}{I_O}\theta = 0$

Examining the above, the linearized model has the form of a standard, unforced, second-order differential equation. Matching this equation to the canonical form, $\ddot{\theta} + 2\zeta\omega_n\dot{\theta} + \omega_n^2\theta = 0$ , we can see how the various system parameters influence the free response of the pendulum system. More specifically:

(4) $2\zeta\omega_n = \frac{b}{I_O}$

(5) $\omega_n^2 = \frac{(M+m)gl_G}{I_O}$

If we had simplified the pendulum as having its mass concentrated at its end, then $\omega_n \approx \sqrt{g/l}$ .

System identification experiment

In this experiment, we will simply record the angular displacement of the pendulum following its release from rest from some initial displacement. The resulting period and the rate at which the amplitude of the oscillation decays will then be employed to identify $\zeta$ and $\omega_n$ , which in turn can be used to identify $b$ .

Hardware setup

The potentiometer, which is physically attached to the rod of the pendulum (not shown), is connected to the Arduino board as shown. The potentiometer is used for sensing the pendulum's angular position and the Arduino board is simply employed to power the sensor, to acquire the data, and to communicate the data to Simulink. Specifically, the power and ground from the Arduino board are connected to the outermost pins of the potentiometer. The middle pin is the potentiometer's "output" and is read on one of the Analog Inputs of the board.

The potentiometer is in essence a variable resistor whose construction is shown below. The slider internal to the potentiometer breaks the wound resistance into two halves, $R_1$ and $R_2$ . The slider is attached to a shaft that is attached to the pendulum. As the pendulum swings, the slider moves changing the distribution of the total resistance between $R_1$ and $R_2$ . By applying voltage across the total wound resistance, the potentiometer behaves like a voltage divider. Hence the output voltage $e_o$ changes as $R_2$ changes.

Applying Ohm's law, we can generate the following expression for $e_o$ .

(6) $e_o = \frac{R_2}{R_1+R_2}e_i$

Since $R_1 + R_2$ is constant, the output voltage changes linearly with the angular displacement of the pendulum.

Software setup

In this experiment, we will employ Simulink to read the data from the potentiometer and to plot the data in real time. In particular, we will employ the IO package from the MathWorks. For details on how to use the IO package, refer to the following link. The Simulink model we will use is shown below and can be downloaded here, where you may need to change the COM port in the IO Setup block to match location where your Arduino board is connected. This model simply reads the potentiometer data (pendulum angle) via an Analog Read on channel A0 and then displays the stored data on a scope and a display. The model also writes the data to the MATLAB workspace for further analysis. The Arduino Analog Read block, the Arduino IO Setup block, and the Real-Time Pacer block are all part of the IO package. The remaining blocks are part of the standard Simulink library, specifically, they can be found under the Sinks library.

Double-clicking on the Analog Read block, we can change the "Sample time." The fastest that the system can be sampled, while allowing real-time communication and plotting with Simulink, is about once every 0.01 seconds. If you try to sample too quickly, the Simulink model will run slower than real-time, that is, it won't be able to sample at the rate specified. For our pendulum length of $l = 0.43\ \mbox{m}$ , a sample time of 0.01 is sufficient. It is not, however, recommended to use a pendulum that is too much shorter since it will cause the dynamics of the pendulum to be too fast (too few samples per period of the pendulum). This can be seen by examining the expression for $\omega_n$ , which in the case of the simplified pendulum model is approximately $\sqrt{g/l}$ .

Once the Simulink model has been created, it can then be run to collect a set of data like that shown below. Specifically, the pendulum is held at rest at an angle just less than 30 degrees. The Simulink model is then started and once it is done initializing and begins recording data, the pendulum can then be released. In choosing the initial pendulum angle, it is desired to choose a sufficiently large angle such that the resulting amplitude is large compared to the level of quantization (and possibly noise). However, the initial angle can't be too large because we wish to employ our linearized model for the purposes of identifying the system's parameters. The linearized model was derived employing a small angle approximation that is accurate only for angles near 0 degrees.

In the above figure, the resulting angles are expressed in numbers of bits. The Arduino Board employs a 10-bit analog-to-digital converter. This means (for the default) an Analog Input channel reads a voltage between 0 and 5 V and slices that range into $2^{10} = 1024$ pieces. Therefore, in the above figure 0 corresponds to 0 V and 1023 corresponds to 5 V. As you can see, we are only using a small percentage of the channel's range. In order to achieve increased sensitivity, we could apply more than 5 V to the potentiometer. This would require an external power source and we would need to orient the potentiometer so that the range through which the pendulum swung through kept the output voltage below 5 V. Another alternative, is to employ a different (smaller) Analog Reference than the default 5 V in order to reduce the range of the input channel. For example, it is possible to configure the board to use the internal 3.3 V source as the Analog Reference.

Aliasing

In order to illustrate what happens when we sample a periodic signal too slowly, let us change our Simulink model sample time to 1 second. If we then extend the model run length to 30 seconds and give the pendulum a larger initial angle, we can generate data like that shown below. Examination of the following seems to indicate we have a response that oscillates with decaying amplitude as we had previously, however, the period appears to be much larger. This phenomenon is referred to as aliasing and demonstrates that sampling a periodic signal too slowly will result in the sampled signal appearing to have a larger period than it actually has. In order to avoid this phenomenon, one needs to sample a periodic signal at a frequency at least twice the signal's highest frequency (referred to as the Nyquist criterion), though it is recommended to sample at even a higher rate than this in order to more accurately recreate the signal being sampled.

Below are a couple more examples to help illustrate the phenomenon of aliasing. For example, if a periodic signal was sampled once per cycle, then the input signal would appear stationary as in the following figure. In this figure, the blue line is the signal being sampled, the dots represent samples, and the red line is the signal that would be reconstructed based solely on the sampled data.

Another example is shown below, where the signal is a sinusoid that is being sampled at a frequency of about 5/3 the signal's frequency. Again, the reconstructed signal has a lower frequency than the actual signal.

Calibration

Let's now return to the case where we sampled our signal at a sufficiently high rate. In order to identify our system parameters, we don't need to convert the measured angles into engineering units that have physical meaning (such as degrees or radians). This is because we can identify the pendulum's parameters based on its period of oscillation and based on the ratio of its amplitudes at different times. Neither of these quantities require a conversion to engineering units. However, we want to be able to validate the accuracy of our resulting model, hence, we need to perform the conversion.

In order to perform the conversion into engineering units, we will calibrate our potentiometer to determine the mapping from bits to actual angles. The way we do this is to simply to hold the pendulum at a defined angle (as measured by a protractor, etc.) while running the Simulink model and recording the corresponding output in bits. Once this has been done at several different angular positions, one can then fit a curve to the data. For our example system, we recorded the following data where straight down is 0 degrees and the clockwise direction is positive.

Angular Position (degrees)	-45	-30	-15	0	15	-30	45
Output (bits)	19	97	177	261	324	396	475

Once the necessary data has been recorded, the MATLAB command polyfit can be employed to fit a curve to the data (in a least squares sense). Specifically, we will fit a polynomial of order n = 1 because we anticipate a linear relationship between the output voltage from the potentiometer and the pendulum's angular position (as described above). Entering the following commands at the MATLAB command line will generate the output shown below.

x = [19 97 177 261 324 396 475]; % output in bits
y = [-45:15:45];                 % angular position in degrees
n = 1;                           % order of the polynomial to be fit
p = polyfit(x,y,n)               % coefficients of the fitted polynomial

p =
    0.1985  -49.6077

In the above, the first element of p is the slope of the fitted line and the second element is the y-intercept. We can further get a visual sense for how well the resulting curve fits our data by plotting them against one another. Running the following additional commands in MATLAB generates the figure below which shows pretty good agreement between the calibration curve and the data.

xfit = [0:1:500];
yfit = p(1)*xfit + p(2);
plot(xfit,yfit,x,y,'ro')
xlabel('Output Voltage (bits)')
ylabel('Angular Position (degrees)')
legend('fitted curve','recorded data','Location','NorthWest')
text(200,-20,'y = 0.1985x - 49.6077')

This calibration data can then be added to our Simulink model, or could be applied to the data we already took.

Parameter identification

When we previously ran our Simulink model with the pendulum swinging, we wrote our data to the workspace under the variable angle. Entering the code: plot(p(1)*angle+p(2)) at the MATLAB command window will convert this data from bits to degrees and will plot the result as shown below. The additional labels and annotation were added separately.

In the above figure, datatips have been added at each of the peaks of the oscillating signal. We can use these specific points of data to estimate our system's parameters. Identifying the first peak by the point $(t_1,x_1)$ , the second peak by the point $(t_2,x_2)$ , and so on, we can perform the following calculations. The damped natural frequency $\omega_d$ of the system can be estimated from the period $T$ of the response:

(7) $\omega_d = \frac{2\pi}{T} \approx \frac{2\pi}{t_{i+1}-t_i}$

From the graph above, the period is consistently very near 1.3 seconds which translates to $\omega_d \approx 4.8\ \mbox{rad/sec}$ . The damping ratio $\zeta$ of the system can be estimated by employing the logarithmic decrement as shown in the equation below.

(8) $\zeta = \frac{\ln\frac{x_i}{x_{i+1}}}{\sqrt{4\pi^2+\left(\ln\frac{x_i}{x_{i+1}}\right)^2}}$

Employing each successive pair of peaks, the calculated $\zeta$ jumps around quite a bit. One possible reason for this is the quantization that arises from the limitations of the analog-to-digital conversion. The potentiometer employed in this experiment has a range of approximately 200 degrees. Dividing this number by 1024 means that angular data is rounded to the nearest 0.20 degree increment. This amount of quantization is significant when compared to the differences in amplitude observed in the above graph. One means for addressing this limitation is to calculate the $\zeta$ across several periods thereby making the same level of quantization a smaller percentage of the difference in amplitudes. This can be done employing the following equation.

(9) $\zeta = \frac{\frac{1}{n-1}\ln\frac{x_1}{x_{n}}}{\sqrt{4\pi^2+\left(\frac{1}{n-1}\ln\frac{x_1}{x_{n}}\right)^2}}$

In the case of our experiment, $n = 6$ and the $\zeta$ is averaged over $n - 1 = 5$ periods. Based on the above, we then estimate that $\zeta \approx 0.004$ .

Another reason that the calculated $\zeta$ values may be jumping around a bit is because the structure of the assumed friction model may not be correct. For example, the pendulum system likely experiences some stiction that causes the pendulum to "stick" a bit when it changes direction. A friction model that captures this stiction phenomenon that you may be familiar with is the Coulomb friction model (the friction has constant magnitude and opposes the direction of motion). Whereas viscous friction causes the amplitude to decay exponentially, Coulomb friction causes the amplitude to decay linearly (you can prove this to yourself).

We will later attempt to further validate our model by building a simulation of the pendulum to compare to our experimental data.

Based on our estimates of $\omega_d$ and $\zeta$ , we can also estimate $\omega_n$ by employing the following relationship.

(10) $\omega_n = \frac{\omega_d}{\sqrt{1 - \zeta^2}}$

Since the estimated damping ratio is so small ( $\zeta \approx 0.004$ ), the system's undamped natural frequency is approximately equal to its damped natural frequency, $\omega_n \approx \omega_d$ . Recall the following relationships we derived by matching the linearized pendulum model to the standard form of a second-order system.

(11) $2\zeta\omega_n = \frac{b}{I_O}$

(12) $\omega_n^2 = \frac{(M+m)gl_G}{I_O}$

Employing these expressions, we can estimate the damping coefficient $b$ (and verify our assumed $I_O$ ). Since we have fair confidence in our measurements of $M$ , $m$ , and $l$ , we will assume that they are accurate. Therefore, estimating the pendulum's mass moment of inertia using the above equation, we have the following.

(13) $I_O = \frac{(M+m)gl_G}{\omega_n^2} \approx \frac{(0.475\ \mbox{kg})(9.81\ \mbox{m/s}^2)(0.39\ \mbox{m})}{(4.8\ \mbox{rad/sec})^2} \approx 0.079 \mbox{kg-m}^2$

This result is pretty close to our theoretical estimate, $I_O = 0.076\ \mbox{kg-m}^2$ . Using the experimentally derived $I_O$ , we can then estimate the damping coefficient $b$ as follows.

(14) $b = 2\zeta\omega_n I_O \approx 2(0.004)(4.8\ \mbox{rad/sec})(0.079 \mbox{kg-m}^2) \approx 0.003 \mbox{N-m-s}$

Model validation

Recalling the model we derived from first principles earlier.

(15) $I_O\ddot{\theta} + b\dot{\theta} + (M+m)gl_G\sin\theta = 0$

And then substituting our estimates of the various parameters, our model becomes the following.

(16) $0.079\ddot{\theta} + 0.003\dot{\theta} + 1.82\sin\theta = 0$

We will validate this model by building a corresponding Simulink model. The results of the simulation model can then be compared to actual experimental data. The Simulink model that includes the nonlinear model and the linearized model is shown below and can be downloaded here. The initial conditions can be set by double-clicking the Integrator blocks. Based on our experimental data, the pendulum was released from an angle of approximately 24.23 degrees. Converting to radians, the initial angle is 0.423 radians. The initial angular velocity (theta_dot) is zero since the pendulum was released from rest.

Defining the following variables in the MATLAB workspace will allow us to run the simulation above. Also, the To Workspace blocks are set to have Array outputs.

M = 0.380;                  % mass of the pendulum bob (kg)
m = 0.095;                  % mass of the pendulum rod (kg)
l = 0.43;                   % length of the pendulum rod (m)
lG = (M*l+0.5*m*l)/(M+m);   % location of pendulum mass center (m)
IO = 0.079;                 % estimate of pendulum mass moment of inertia (kg-m^2)
b = 0.003;                  % estimate of viscous friction coefficient (N-m-s)
theta_ic = 0.423;           % initial pendulum angular position (rad)
theta_dot_ic = 0;           % initial pendulum angular velocity (rad/s)
g = 9.81;                   % acceleration due to gravity (m/s^2)

Running the above simulation, we can then take the outputs theta_lin and theta_nl and compare them to the actual experimental data. Since in the experiment the pendulum was released at approximately 1.9 seconds, we will shift the simulated outputs by 1.9 seconds to align them with the experimental data.

         plot(p(1)*angle+p(2));
         xlabel('time (seconds)')
         ylabel('angle (degrees)')
         title('Pendulum Free Response')
         hold
         plot(time+1.9,theta_lin*180/pi,'g:');
         plot(time+1.9,theta_nl*180/pi,'r--');
         legend('experiment','linear sim','nonlinear sim')

Examining the above, one can see there is pretty good, though not perfect, agreement between the models and the experimental data. We've discussed some sources of error previously (quantization, limitations of the friction model, etc.). Another source is the fact that we used a linearized model for estimating some of the system parameters. Evidence of this effect can be seen somewhat by the fact that the linearized model and nonlinear model don't agree perfectly. In fact, the approximate linearized model agrees with the experimental data better than the nonlinear model. One way to improve this agreement is to estimate the parameters from an experiment where the pendulum is given a smaller initial displacement. Since the linearization was performed about the angle 0, a smaller amplitude would keep the data closer to the point of the linearization, and hence would make the approximate linearized model more accurate.

Extensions

In this activity we modeled and analyzed a pendulum system. A logical extension would be to design and implement a control system to maintain the pendulum angle at some commanded level. This could be done by attaching a motor directly to the pendulum via a gear or pulley system. Alternatively, a motor could be outfitted with propellers and then attached to one end of the pendulum system (where the end mass is used as a counterweight to lessen the demand on the motor). The speed of the motor is then controlled to generate lift to affect the attitude of the pendulum assembly. One such instantiation of this setup is shown below.