A small step: Modelling pendulum

Erratum: I have now fixed the errors that were present in the calculations in the previous version of the post.

This term I am teaching a course on modelling dynamical systems. In this post I will present a simple example of modelling a cart-pendulum system (see Fig. 1).

Fig. 1. A cart-pendulum system.

The system consists of a pendulum mounted on a cart sliding horizontally. The pendulum is assumed to be a point mass hanging on a massless wire at an angle of θ The cart is pushed by some external force u=f(t) and its position is defined as z The following parameters define the properties of the system:mass of the cart m1, mass of the pendulum m2 and the length of the pendulum l.

The problem is formulated as follows. Find the equations of motion for the system using Lagrange equations. Linearize the non-linear model obtained. Observe and compare the behaviour of both the linear and the non-linear models.

We will solve it in the following steps:
1) How many degrees of freedom (DoF) do we have here?
2) What form of Lagrange equations should we use?
3) What’s the kinetic energy (K) and the potential energy (U) of the system?
4) Solve Lagrange equations.
5) Define state-space non-linear model.
6) Linearize the model.
7) Simulate the behaviour of the models using Matlab.

Let’s go!

1) First, it is quite obvious that the system has two degrees of freedom: z and θ.
It ought to be said that there are several ways these generalized coordinates can be chosen, but this is the most convenient choice in this problem.

2) The Lagrange equations for this system are going to be in the following form:

$\frac{d}{dt} \frac{\partial\mathcal{L}}{\partial\dot{q_i}}\,-\,\frac{\partial\mathcal{L}}{\partial{q_i}}\,=\,F_{q_i}$

3) Since the Lagrangian of the system is defined as its kinetic energy K minus its potential energy U, we have to find out what K and U are. There are two bodies in the system that have masses (m1 and m2). K will therefore be:

$K\,=\,\frac{1}{2}m_1 v_1^2\,+\,\frac{1}{2}m_2 v_2^2$

We have to be careful while calculating velocities v1 and v2. The velocity v1 is straightforward, but v2 consists of the carrying linear motion of the cart plus the motion of the pendulum relative to the cart. Moreover, the velocity of the pendulum around the cart can be split into two orthogonal vectors: $v_x$ and $v_y$ .

$v_1\,=\,\dot{z}_1$

$\bar v_2\,=\,\bar v_1 + \bar v_{2/1}\,=\,(v_1+v_{x})\,\bar i \,+\, (v_y)\,\bar j$

Since $v_x$ and $v_y$ are calculated as follows:

$v_x\,=\,\dot\theta l \cos\theta$
$v_y\,=\,\dot\theta l \sin\theta$

The $v_2^2$ will be:

$v_2^2\,=\,(\dot z+\dot\theta l \cos\theta)^2\,+\,(\dot\theta l \sin\theta)^2$
$v_2^2\,=\,\dot z^2\,+\,2\dot z \dot\theta l \cos\theta\,+\,\dot\theta^2l^2$

The kinetic energy K will therefore be:

$K\,=\,\frac{1}{2}m_1\dot{z}^2\,+\,\frac{1}{2}m_2(\dot z^2\,+\,2\dot z \dot\theta l \cos\theta\,+\,\dot\theta^2l^2)^2$

The potential energy U is associated with the mass of the pendulum m2 hanging in the gravitational field. Here, we are going to introduce temporary variable y which describes the vertical displacement of the pendulum:

$U\,=\,-m_2 g y \,=\, -m_2 g l\,\cos \theta$

4) Having calculated both K and U, we can formulate the Lagrangian:

$\mathcal{L} \,=\, K\,-\,U\,=\, \frac{m_1+m_2}{2}\,\dot z^2\,+\,m_2\dot z \dot\theta l\cos\theta\,+\,\frac{m_2}{2}\dot\theta^2l^2\,+\,m_2 g l \cos\theta$

In order to solve the equation, we have to calculate these pesky partial derivatives - first for z, and then for θ. I did this in WolframAlpha.

$\frac{d}{dt}\frac{\partial\mathcal{L}}{\partial\dot{z}}\,=\,(m_1+m_2)\ddot z\,+\,m_2 l (\ddot\theta\cos\theta - \dot\theta^2\sin\theta)$

$\frac{\partial\mathcal{L}}{\partial{z}}\,=\,0$

$\frac{d}{dt}\frac{\partial\mathcal{L}}{\partial\dot{\theta}}\,=\,m_2l(l \ddot\theta+\ddot z\cos\theta-\dot z \dot\theta\sin\theta)$

$\frac{\partial\mathcal{L}}{\partial{\theta}}\,=\,-m_2l(\dot z \dot\theta+g)\,\sin\theta$

Now, subsituting the derivatives and the generalized force u into the Lagrange equations, we get a system:

$(m_1+m_2)\ddot z\,+\,m_2l(\ddot\theta\cos\theta-\dot\theta^2\sin\theta)\,=\,u$
$\ddot z l \cos\theta+\ddot\theta l^2+gl\sin\theta\,=\,0$

We are not going to calculate this by hand. I have used the following WolframAlpha query (and then substituted the names of variables back):

solve for x, y: (m+n)x+nl(ycos(f)-w^2sin(f))=u, xlcos(f)+yl^2+glsin(f)=0

The equations of motion are:

$\ddot{z}\,=\, \frac{m_2 \sin \theta \cdot (g\cos \theta+l\dot\theta^2)\,+\,u}{-m_2\cos^2 \theta +m_1 +m_2}$

$\ddot\theta \,=\, - \frac{g\sin\theta \cdot (m_1+m_2)\,+\,\cos\theta \cdot (lm_2\dot\theta^2\sin\theta + u)}{l\cdot(-m_2\cos^2\theta+m_1+m_2)}$

5) We can now construct the (non-linear) state-space model. This is done by first naming the state-space variables:

$x_1\,=\,z$
$x_2\,=\,\dot{z}$
$x_3\,=\,\theta$
$x_4\,=\,\dot\theta$

and then substituting into the equations above:

$\dot{x}_1\,=\,x_2$
$\dot{x}_2\,=\,\frac{m_2 \sin x_3 \cdot (g\cos x_3\,+\,lx_4^2) \,+\,u}{-m_2 \cos^2 x_3\,+\,m_1\,+\,m_2}$
$\dot{x}_3\,=\,x_4$
$\dot{x}_4\,=\,-\frac{g\sin x_3 \cdot (m_1+m_2)\,+\,\cos x_3 \cdot (lm_2x_4^2\,\sin x_3\,+\,u)}{l \cdot (-m_2\,\cos^2 x_3 + m_1 + m_2)}$

This non-linear system is of form:

$\bf{\dot{x}}\,=\,\bf{f}(\bf{x},\,\bf{u},\,\bf{t})$
$\bf{y}\,=\,\bf{g}(\bf{x},\,\bf{u},\,\bf{t})$

where, for instance:

$f_1(\bf{x},\,\bf{u},\,\bf{t})\,=\,x_2$

6) Linearization of the system is done by replacing the non-linear functions with their first-order approximations at selected “operating point”. Mathematically this is one by assuming the following linear model:

$\Delta\bf{\dot{x}}\,=\,A\Delta\bf{x}\,+\,B\Delta\bf{u}$
$\Delta\bf{y}\,=\,C\Delta\bf{x}\,+\,D\Delta\bf{u}$

The Δx, Δu and Δy are defined as offsets from the operating point:

$\Delta x\,=\,x\,-\,x|_0$
$\Delta u\,=\,u\,-\,u|_0$
$\Delta y\,=\,y\,-\,y|_0$

Let’s define our operating point as:

$x|_0\,=\,[0,\,0,\,0,\,0]^T$
$u|_0\,=\,0$

The matrices A, B, C and D are calculated as:

$A\,=\,\nabla_x\,\bf{f}|_0$
$B\,=\,\nabla_u\,\bf{f}|_0$
$C\,=\,\nabla_x\,\bf{g}|_0$
$D\,=\,\nabla_u\,\bf{g}|_0$

It isn’t as scary as it looks. $\nabla_x\,\bf{f}$ just says that we have to calculate the derivatives of functions $f_1$ , $f_2$ and so on partially over $x_1$ , $x_2$ etc. and put them into the matrix A. The $|_0$ says: after calculating the derivative, substitute $x_0$ and $u_0$ (assumed previously) into the equation.

I have used the following WolframAlpha queries to compute the (non-trivial) derivatives A23, A24 and A43, A44:

d/dx (nsin(x)(gcos(x)+ld^2)+u)/(-ncos(x)^2+m+n)

d/dy (nsin(x)(gcos(x)+ly^2)+u)/(-ncos(x)^2+m+n)

d/dx -(gsin(x)(m+n)+cos(x)(lnd^2sin(x)+u))/(l(-ncos(x)^2+m+n))

d/dy -(gsin(x)(m+n)+cos(x)(lny^2sin(x)+u))/(l(-ncos(x)^2+m+n))

$A_{11}\,=\,\frac{\partial f_1}{\partial x_1}|_0\,=\,0$
$A_{12}\,=\,\frac{\partial f_1}{\partial x_2}|_0\,=\,1$
$A_{13}\,=\,\frac{\partial f_1}{\partial x_3}|_0\,=\,0$
$A_{14}\,=\,\frac{\partial f_1}{\partial x_4}|_0\,=\,0$

$A_{21}\,=\,\frac{\partial f_2}{\partial x_1}|_0\,=\,0$
$A_{22}\,=\,\frac{\partial f_2}{\partial x_2}|_0\,=\,0$
$A_{23}\,=\,\frac{\partial f_2}{\partial x_3}|_0\,=\,m_2 g$
$A_{24}\,=\,\frac{\partial f_2}{\partial x_4}|_0\,=\,0$

$A_{31}\,=\,\frac{\partial f_3}{\partial x_1}|_0\,=\,0$
$A_{32}\,=\,\frac{\partial f_3}{\partial x_2}|_0\,=\,0$
$A_{33}\,=\,\frac{\partial f_3}{\partial x_3}|_0\,=\,0$
$A_{34}\,=\,\frac{\partial f_3}{\partial x_4}|_0\,=\,1$

$A_{41}\,=\,\frac{\partial f_4}{\partial x_1}|_0\,=\,0$
$A_{42}\,=\,\frac{\partial f_4}{\partial x_2}|_0\,=\,0$
$A_{43}\,=\,\frac{\partial f_4}{\partial x_3}|_0\,=\,\frac{-g(1_1+m_2)}{m_1 l}$
$A_{44}\,=\,\frac{\partial f_4}{\partial x_4}|_0\,=\,0$

Similarly, B would be:

$B_{11}\,=\,\frac{\partial f_1}{\partial u_1}|_0\,=\,0$
$B_{21}\,=\,\frac{\partial f_2}{\partial u_1}|_0\,=\,\frac{1}{m_1}$
$B_{31}\,=\,\frac{\partial f_3}{\partial u_1}|_0\,=\,0$
$B_{41}\,=\,\frac{\partial f_4}{\partial u_1}|_0\,=\,-\frac{1}{m_1 l}$

In order to observe all state variables, define C = I. D has to be filled with zeros, with dimensions that agree with B and C.

7) We can now see how these work!
The scripts for this problem can be found at: https://github.com/dagothar/mos03

By default, they generate the object response with u=0, but with a starting condition in which the pendulum starts at a positive angle to the vertical. And here are the results:

Fig. 2. Reponse of the non-linear system (to the left) and of the linearized system (to the right).

The responses seem very similar when a small initial angle displacement is given. Things change when the pendulum is placed at a larger starting angle.
Can you guess what happens?

2017-03-11

Modelling pendulum

No comments:

Post a Comment