Modelling an Inverted Pendulum in Matlab

Uncategorized 7 Minutes

Modelling an Inverted Pendulum

There are a number of ways to model and control an inverted Pendulum using MATLAB and Simulink. Each approach has its benefits. Four approaches will be outlined in the following posts, two in MATLAB and two in Simulink.

The modelling of an inverted Pendulum is a useful exercise to familiarise oneself with MATLAB and Simulink while applying your knowledge of dynamics. It is a particularly good exercise since the inverted pendulum problem is a practical practical problem which closely relates to a number of common place engineering problems.

\bf {MATLAB \ Command \ Approach}

We will model the inverted pendulum system directly using the MATLAB command window, by first deriving the necessary equations with which we can define the system and then inserting the equation/s into MATLAB, to allow the software to model the system for us. This is obviously far more math heavy than the Simulink approaches we will cover later on.

Two ways of representing the motion of the pendulum and cart in MATLAB are as transfer functions or in a state-space form. In either case we will need to look at the dynamics of the system and resolve the forces acting on the Pendulum and cart, for an impulsive external force F. A simple inverted Pendulum attached to a wheeled cart is shown below:

Figure 1: Illustration of a Simple inverted Pendulum loaded onto a cart
(The figure was take from Michigan University Control Tutorial page)

\bf {Transfer \ Functions}

From the Figure above we see that the sum of the horizontal forces acting on the cart give the following equation of motion:

m \ddot x + b \dot x + N = F

Where b is the coefficient of friction on the between the car and the ground, and N is the reaction force on the cart. For now we will only look at the sum of the horizontal forces.

Similarly, we then sum up the forces acting on the pendulum to show that the reaction force N for the pendulum can be defined by the angular velocity the equation:

N= m \ddot x + mL \ddot \theta cos \theta - mL \dot \theta ^2 sin \theta

By substituting the equation for N back into the original equation of motion of the cart, we obtain:

(M+m)\ddot x + b \dot x + mL \ddot \theta cos \theta - mL \dot \theta ^2 sin \theta = F

The sum of the forces perpendicular to the pendulum can be described by

Psin \theta + Ncos \theta - mg sin \theta = mL \ddot \theta + m \ddot x cos \theta

In order to substitute out the reaction force terms P and N, we can simply use the equation for the sum of moments about the Pendulum, which is a function of the reaction forces and the moment of inertia I:

I \ddot \theta = - PL Sin \theta - NL cos \theta

By substituting this into the equation describing the forces about the pendulum, above, we can rearrange to get:

(I + mL^2) \ddot \theta =mgLsin \theta = - mL \ddot x cos \theta

As we are, at least for the purpose of this exercise, only applying this control design to a linear system, we can linearize the equations, such that \theta = \pi + \phi , where the \phi defines the deviation of the pendulum position from the upright equilibrium position. This is a reasonably valid assumption provided that under the control law we implement, the pendulum does not deviate any more than 20 degrees.

On this basis, we can also then linearize the trigonometric terms in the equation such that:

cos \theta = cos (\pi + \phi ) \approx -1

sin \theta = sin(\pi + \phi ) \approx - \phi

The linearized forms of the equations of motion derived will give the following:

(I + mL^2) \ddot \phi- mgL \phi = mL \ddot x

and

(M+m)\ddot x + b \dot x - mL \ddot \phi = u

Writing the Transfer Function of the pendulum and cart:

In order to compute the transfer function describing the linearized system, we must take the Laplace transform of the equations, assuming zero initial conditions, such that u(0)=0, x(0)=0 and \Phi (0)=0.

The resulting Laplace transforms of the two equations are shown below:

(I+mL^2) \Phi (s)s^2 - mgL\Phi (s) = mLX(s) s^2

and

(M+m) X(s) s^2 + bX(s)s -mL \Phi (s)s^2= U(s)

Note: A transform function represents the relationship between a single input and single output at a time. Therefore, we will need to compute a separate transfer function for the position of the cart and the angular deviation of the pendulum.

We will then need to rearrange the derived equations in terms of it’s inputs and outputs such that the Transfer function is defined: T(s)= \frac { \Phi (s)}{ U(s)}

The transfer function of the pendulum will be defined by the input U(s) and output \phi (s). Therefore, we will need to substitute out the X(s) term of the cart position in order to compute the transfer function in terms of a single input, rotation angle \Phi (s). We can use the previous two equations derived to obtain an expression of X(s) in terms of \Phi (s).

X(s)=\left [ \frac{I+mL^2}{mL}- \frac {g}{s^2} \right] \Phi(s)

This expression can then be used to substitute out the X(s) term in Laplace transform of the equation of motion to obtain an equation containing only a one input term \Phi(s) and output term U(s):

(M+m) \left [ \frac{I+mL^2}{mL}- \frac {g}{s^2} \right] \Phi (s)^2 + b \left [ \frac{I+mL^2}{mL}- \frac {g}{s^2} \right] \Phi (s)s -mL \Phi (s)s^2 = U(s)

We can then rearrange this equation to obtain an expression of the transfer function of the form:

T(s) = \frac {\Phi(s)}{U(s)} = \left ( \frac {1}{mLs^2} \right ) ([(M+m)(I+mL^2) - (mL)^2]s^4+b(I+mL^2)s^3-[(M+m)mgL]s^2 - [bmgl]s )^{-1}

The transfer function equation contains a large number of constants (M, m, L and I) which substantially lengthen the transfer function expression and increase the possibility for error. For our own ease of use, we can shorten the transfer function by simply substituting out some of these terms with a constant q.

For example we can take q = [(M+m)(I+mL^2)-(mL)^2

By substituting out these terms for q, we get a less chaotic expression:

T_{pend}(s)= \frac {{mL}s^2}{qs^4+ {b(I+mL^2)}s^3- {(M+m)mgL}s^2-{bmgL}s}

We can see that in this transfer function, for the motion of the Pendulum, there are common real zeros and poles. This can be simply canceled out to give:

T_{pend}(s)= \frac {{mL}s^2}{qs^3+ {b(I+mL^2)}s^2- {(M+m)mgL}s-{bmgL}}

This is the final transfer function for modelling Pendulum orientation. When writing this in MATLAB, it is important to remember to define the constant ‘q’, use to simplify our derivation, in the workspace along with the other constants of the system before running the simulation.

In addition to the Transfer function of the pendulum, we can proceed to simply formulate the transfer function of the cart, to give cart position as the output, using the above expression by simply substituting back in the variable \Phi (s) in terms of X(s). this will give the following:

T_{cart}(s) = \frac {X(s)}{U(s)} = \frac {{mL}\left [ \frac{I+mL^2}{mL}- \frac {g}{s^2} \right]s^2}{qs^4+{b(I+mL^2)}s^3-{(M+m)mgL}s^2-{bmgL}s}

We can now start transferring our work to MATLAB. We will need to start by defining our known constants (m, M, L, g, I, b and q) in the workspace, before transferring across the two transfer function expressions.

The steps for modelling the system in MATLAB are outlined further in the ‘Modelling in MATLAB’ Section below.

The command for modelling a system based on a transfer functions is:

> sys \_ tf(Tpend, Tcart)”>

\bf {State-Space \ Form}

We can also represent the model in state space form. To do this we need to identify the states with which we can define the system and compute the representative matrix. This method may seem more direct at first as it involves fewer steps and derivations but can be prove to be more mathematically tedious than the transfer function method, as it may result in having to work out more complex matrices, increasing the possibility for error.

Definitions: State space for is the Euclidean space in which the variable on the axis are the state variables, the state of the system can be represented a a vector within the space. As a result, when expressing the dynamics of the system in state-space form, the system in the inputs, outputs and states are also expressed as vectors.

We assign a state to each of the elements of the system’s dynamics and look at their dynamics separately.

In this case, using the two equations of motion previously derived for the inverted pendulum, stated above, we can identify 4 principal states that we can use to define the behaviour of the system: \dot x, \ddot x, \dot \theta, \ddot \theta

We first need to define these states using the two equations of the system derived from the sum of forces acting on the cart and pendulum. It is important to note that in this instance we are still assuming a linearised system whereby the following is still true:

(I+mL^2)\ddot \phi - mgl \phi = mL \ddot \phi

(M+m) \ddot x +b \dot x- mL \ddot \phi = u

We can use the two above equations to formulate expressions for the four states. The states \dot x and \dot \phi can simply remain as they are. However, the other two states will need to be defined in terms of the former two states and the input u.

\ddot \theta = \frac {1}{(I+mL^2)} (mL \ddot x + mgL \phi)

\ddot x = \frac {1}{(M+m)}(u+mL\ddot \phi - b \dot x)

It is important to note that each of the state expressions can only be defined in terms of constants and lower order derivatives i.e. \ddot x cannot be defined in term of \ddot \phi and vice versa. Using the two equations above we can substitute out the unwanted variables from the respective expressions, and with a bit of rearranging, show that:

\ddot \phi = - \frac {mLb}{(I+mL^2)} \dot x + \frac {(M+m)mgL}{I(M+m)+MmL^2} \phi + \frac {mL}{I(M+m)+MmL^2}u

\ddot x = -\frac {I+mL^2}{I(M+m)+MmL^2}b+\frac {m^2L^2g}{I(M+m)+MmL^2}\phi+ \frac {I+mL^2}{I(M+m)+MmL^2}u

From this we can write the state space equation in matrix form:

\left [\begin {matrix} \dot x \\ \ddot x \\ \dot\phi \\ \ddot\phi \end {matrix} \right]= \left [ \begin {matrix} 0&1&0&0 \\ 0&\frac{-(I+mL^2)b}{I(M+m)+MmL^2}&\frac {m^2L^2}{I(M+m)+MmL^2}&0 \\ 0&0&0&1 \\ 0&\frac {-mLb}{I(M+m)+MmL^2}&\frac {(M+m)mgL}{I(M+m)+MmL^2}&0 \end {matrix} \right ] \left [\begin {matrix} x \\ \dot x \\ \phi \\ \dot\phi \end {matrix} \right ] + \left [ \begin {matrix} 0 \\ \frac {I+mL^2}{I(M+m)+MmL^2} \\ 0 \\ \frac {mL}{I(M+m)+MmL^2} \end {matrix} \right ] u

The output of the state-space can be expressed in the matrix form as follows:

y=\left[ \begin {matrix} 1&0&0&0 \\ 0&0&1&0 \end {matrix} \right ] \left [ \begin {matrix} x \\ \dot x \\ \phi \\ \dot\phi \end {matrix} \right ] + \left [ \begin {matrix} 0 \\ 0 \end {matrix} \right ] u

Where the first and second rows of the output y are the position of the cart ad angle of rotation of the pendulum.

To transfer this over to MATLAB, we first need to define all our known constants, variables and matrices. For simplicity, we can express the state space equation in the form:

\left [\begin {matrix} \dot x \\ \ddot x \\ \dot\phi \\ \ddot\phi \end {matrix} \right]=A\left [\begin {matrix} x \\ \dot x \\ \phi \\ \dot\phi \end {matrix} \right ]+Bu and y=C\left [\begin {matrix} x \\ \dot x \\ \phi \\ \dot\phi \end {matrix} \right ]+ \matrix Du

So as not to overburden ourselves when writing in the expression.

We can then define our states, inputs an outputs. To model the system using the state-space form we can use the command:

>”>sys_ss=ss (A,B,C,D, ‘statename’, states, ‘inputname’, inputs. ‘outputname’, outputs)

This state space form of the system can also be converted into a transfer function using the command: >>sys_tf=tf(sys_ss)

\bf {Modelling \ the \ System \ in \ MATLAB}

In order to begin modelling the system in MATLAB, no matter which approach is taken, after the known constants and variables of the system have been defined in the workspace.

We then need to define our inputs and outputs as shown below:

>input=\{‘u’\}”>

>output=\{‘x’,’phi’\}”>

We can then use the relevant system command, of the the respective type of expression, to model our system i.e. sys_tf or sys_ss for the transfer function and state space-form expression, respectively.

We can set the names for the input and outputs using the commands:

>set(sys_tf, ‘inputname’, inputs)”>

>set(sys\_tf, ‘outputname’, outputs)”>

In the specific case of the state space method we will also need to define our 4 states before inserting the system command:

>states=\{‘x’, ‘xdot’, ‘phi’, ‘phidot’\}”>

The Matlab script used for to model either cases are presented below. However, this could always be written in a different format. I have also added the script you can use to produce graphs showing the model’s response to an arbitrary impulsive response.

Script for modelling the inverted Pendulum using transfer functions in the Laplace domain:

Script for modelling the inverted Pendulum in the state space-form is shown below:

We ca use the same impulsive function and the script as shown in the ‘transfer function’ code to test the function for an impulsive force and produce graphical representation of the response of the cart’s movement and angular deviation of the Pendulum.

Unknown's avatar

Published by Majed Yassin

Published

Leave a comment