ECE 602 Fall 2010

Dynamic System Modeling

by

Stanislaw H. Żak

September 07, 2010

To catch a mouse alive, try this trick: When you see the beast running around a room, open a black umbrella and hold it on the floor at an angle that allows the mouse to run in (it helps if someone else can chase the mouse). Your prey, drawn to dark places, will probably head straight for the umbrella. Then quickly close it, go outside, and set your prisoner free. [3, page 135]

1 Mathematical Modeling

Mathematical models of systems that we wish to control are key elements in their analysis and controller design. We devote this section 1 to mathematical modeling.

The goal of mathematical modeling is to provide a mathematical description of the interrelations between the system quantities as well as the relations between these quantities and external inputs. A mathematical model can be viewed as “a mathematical representation of the significant, relevant aspects of an existing or proposed physical system. (Significance and relevance being in relation to an application where the model is to be used).” [2, page xi]. Another definition of mathematical model that we found in [1, page 1] is: “A mathematical model is a system of equations that purports to represent some phenomenon in a way that gives insight into the origins and consequences of its behavior.”

At this point, we wish to point out that it is important to distinguish between a physical system and its model. A physical system may be modeled using different models. In particular, it may be advantageous to model the physical system using different models for its different operating regions. For example, it is common to model an electric circuit using different models for different ranges of frequencies in which the circuit operates. In what follows, we refer to the physical system’s model as the system.

In general, the more accurate the mathematical model is, the more complex are the modeling equations. This means that during the mathematical modeling process we may be forced to find a middle course between accuracy of the model and its complexity. Our goal is to generate the most accurate possible description of the given system while still being able to Figure 1: Diagram of mathematical modeling process and its use in the controller design process.

carry out the model analysis and computations in a reasonable amount of time. At first, we tend to make oversimplifying assumptions in order to obtain a very simple model. However, we should follow the principle of maximum simplicity, which states that ”everything should be as simple as possible—but no simpler.” After extracting some information about the system behavior based on this first, possibly crude model, it may be desirable to construct a more accurate model. We try to create a mathematical description of a given physical process that incorporates realistic assumptions and physical constraints. We then extract information from this model and compare it with the experimental data. We perform this comparative study to determine the value of the obtained model. If the model appears to be inadequate for our needs, then it must be adjusted. The process of mathematical modeling is depicted in the diagram in Figure 1. In the process of mathematical modeling, one can use different types of equations and their combinations like algebraic, ordinary and partial differential equations, difference, integral, or functional equations. In this book, only difference and ordinary differential equations are used to model dynamic systems.

For the controller design purposes, we first construct a truth model of the given plant and then its design model.

2 Review of Work and Energy Concepts

The purpose of this section is to review and reacquaint ourselves with some basic notions from physics like work, energy, etc. After a brief review of work and energy concepts, we derive the Lagrange equations of motion and use them to devise mathematical models of dynamic systems.

Suppose we are given a particle of a constant mass m subjected to a force F. Then, by Newton’s second law we have

F = m\frac{dv}{dt} = ma.  (1)

Here the force F and the velocity v are vectors. Therefore, they can be represented as

F = [ F_{x_{1}} F_{x_{2}} F_{x_{3}} ]^{⊤} and v = [ v_{x_{1}} v_{x_{2}} v_{x_{3}} ]^{⊤} = [ ẋ_{1} ẋ_{2} ẋ_{3} ]^{⊤} .

Using the above notation, we represent equation (1) as

[ \begin{matrix}F_{x_{1}} \\ F_{x_{2}} \\ F_{x_{3}} \end{matrix} ] = m [ \frac{dv_{x_{1}}}{dt} ; \frac{dv_{x_{2}}}{dt} ; \frac{dv_{x_{3}}}{dt} ] .  (2)

Suppose now that a force F is acting on a particle located at a point A and the particle moves to a point B. The work δW done by F along infinitesimally small distance δs is

δW = F^{⊤}δs
= [ F_{x_{1}} F_{x_{2}} F_{x_{3}} ] [ \begin{matrix}δx_{1} \\ δx_{2} \\ δx_{3} \end{matrix} ]
= F_{x_{1}}δx_{1} + F_{x_{2}}δx_{2} + F_{x_{3}}δx_{3}.

where δs = [ δx_{1} δx_{2} δx_{3} ]^{⊤}. The work W_{AB} done on the path from A to B is obtained by integrating the above equation,

W_{AB} = ∫_{A}^{B} (F_{x_{1}}dx_{1} + F_{x_{2}}dx_{2} + F_{x_{3}}dx_{3}) .

We now would like to establish a relation between work and kinetic energy. For this observe that

ẍ dx = \frac{dẋ}{dt}dx = ẋdẋ.

Using the above relation and Newton’s second law, we express W_{AB} in a different way as

W_{AB} = ∫_{A}^{B} (mẍ_{1}dx_{1} + mẍ_{2}dx_{2} + mẍ_{3}dx_{3})
= ∫_{A}^{B} m (ẋ_{1}dẋ_{1} + ẋ_{2}dẋ_{2} + ẋ_{3}dẋ_{3})
= m (\frac{ẋ_{1}^{2}}{2} + \frac{ẋ_{2}^{2}}{2} + \frac{ẋ_{3}^{2}}{2} )|_{A}^{B}
= m (\frac{∥v_{B}∥^{2}}{2} − \frac{∥v_{A}∥^{2}}{2} )
= \frac{m∥v_{B}∥^{2}}{2} − \frac{m∥v_{A}∥^{2}}{2} .  (3)

where v_{A} = [ ẋ_{1A} ẋ_{2A} ẋ_{3A} ]^{⊤}, and v_{B} is defined in an analogous manner, and ∥v∥^{2} = v^{⊤}v. Observing that \frac{m∥v_{B}∥^{2}}{2} is the kinetic energy of the particle at the point B and \frac{m∥v_{A}∥^{2}}{2} is its kinetic energy at the point A, we obtain

W_{AB} = K_{B} − K_{A}.

The above relation can be used to define the kinetic energy of a particle as the work required to change its velocity from some value v_{A} to a final value v_{B}. The relation, W_{AB} = K_{B} − K_{A} = ΔK is also known as the work-energy theorem for a particle. A force is conservative if the work done by the force on a particle that moves through any round trip is zero. In other words, a force is conservative if the work done by it on a particle that moves between two points depends only on these points and not on the path followed. We now review the notion of potential energy or what is the same, the energy of configuration. Recall that if the kinetic energy K of a particle changes by ΔK, then the potential energy U must change by an equal but opposite amount so that the sum of the two changes is zero, that is,

ΔK + ΔU = 0.

This is equivalent to saying that any change in the kinetic energy of the particle is compensated for by an equal but opposite change in the potential energy U of the particle so that their sum remains constant, that is,

K + U = constant.

The potential energy of a particle represents a form of stored energy that can be recovered and converted into the kinetic energy. If we now use the work-energy theorem and the equation ΔK + ΔU = 0, then we obtain

W = ΔK = −ΔU.

The work done by a conservative force depends only on the starting and the end points of motion and not on the path followed between them. Therefore, for motion in one dimension,

ΔU = −W = − ∫_{x_{0}}^{x} F (s)ds.  (4)

We thus can write

F (x) = −\frac{dU(x)}{dx}.  (5)

A generalization of equation (4) to motion in three dimensions yields

ΔU = − ∫_{x_{10}}^{x_{1}} F_{x_{1}}ds − ∫_{x_{20}}^{x_{2}} F_{x_{2}}ds − ∫_{x_{30}}^{x_{3}} F_{x_{3}}ds,

and a generalization of equation (5) to motion in three dimensions is

F (x) = − [ \frac{∂U}{∂x_{1}}  \frac{∂U}{∂x_{2}}  \frac{∂U}{∂x_{3}} ]^{⊤} = −∇U (x) ,  (6)

where x = [ x_{1} x_{2} x_{3} ]^{⊤}.

In a conservative vector field the work done by it on a particle that moves between two points depends only on these points and not on the path followed.

We now represent Newton’s equation in an equivalent format that establishes a connection with the Lagrange equations of motions to be discussed in the following section. To proceed, note that

K = \frac{m \left∥{ \atop ẋ}\right∥^{2}}{2} = \frac{mẋ^{⊤}ẋ}{2}.  (7)

Hence,

\frac{∂K}{∂ẋ_{i}} = mẋ_{i}, i = 1, 2, 3.  (8)

Combining (6) and (8), we represent Newton’s equations as

\frac{d}{dt} (\frac{∂K}{∂ẋ_{i}} ) + (∇U)_{i} = 0, i = 1, 2, 3.  (9)

Let

L = K − U.  (10)

The function L defined above is called the Lagrangian function or just the Lagrangian. Note that

\frac{∂L}{∂ẋ_{i}} = \frac{∂K}{∂ẋ_{i}} ,  (11)

and

\frac{∂L}{∂x_{i}} = −\frac{∂U}{∂x_{i}} ,  (12)

Substituting (11) and (12) into (9), we obtain

\frac{d}{dt} (\frac{∂L}{∂ẋ_{i}} ) − \frac{∂L}{∂x_{i}} = 0, i = 1, 2, 3.  (13)

Equations (13) are called the Lagrange equations of motions, in Cartesian coordinates, for a single particle. They are just an equivalent representation of Newton’s equations given by (2).

3 The Lagrange Equations of Motion

A single particle with no constraints is free to take up any position in the three dimensional space. This is the reason we need three coordinates to describe its position. In the case of a system of p particles there are 3p degrees of freedom. In general, if there are constraints imposed on the particles, then the number of the degrees of freedom is

n = 3p − number of constraints.

For example, if the system consists of only one particle in space and if the particle is restricted to move along a line, then n = 3 − 2 = 1. We can use many different coordinate systems to describe a location of a particle in space. We will use the letter q as a symbol for a coordinate regardless of its nature. A generalized coordinate may represent, for example, angular displacement or electric charge. This indicates that the results we are going to obtain will be of general nature and will not be restricted to mechanical systems. We refer to q as a generalized coordinate. Thus, if

\begin{matrix}x & = & x(q_{1}, q_{2}, q_{3}) \\ y & = & y(q_{1}, q_{2}, q_{3}) \\ z & = & z(q_{1}, q_{2}, q_{3}) \end{matrix} \}

then

ẋ = \frac{∂x}{∂q_{1}} q̇_{1} + \frac{∂x}{∂q_{2}} q̇_{2} + \frac{∂x}{∂q_{3}} q̇_{3},

etc. We assume that the functions x, y, z ∈ C^{2}, which implies that their corresponding mixed second partial derivatives are equal. Suppose now that a particle is confined to move on a surface. Therefore, to describe the location of the particle, we need only two generalized coordinates. The position of the particle on a surface in a three dimensional space is then described as

x = x(q_{1}, q_{2}),
y = y(q_{1}, q_{2}),
z = z(q_{1}, q_{2}),

and hence

\begin{matrix}δx & = & \frac{∂x}{∂q_{1}} δq_{1} + \frac{∂x}{∂q_{2}} δq_{2} \\ δy & = & \frac{∂y}{∂q_{1}} δq_{1} + \frac{∂y}{∂q_{2}} δq_{2} \\ δz & = & \frac{∂z}{∂q_{1}} δq_{1} + \frac{∂z}{∂q_{2}} δq_{2}. \end{matrix} \}  (14)

We know that the work δW done by the force F is δW = F_{x}δx + F_{y}δy + F_{z}δz. On the other hand, this work is equal to the change in the kinetic energy. Therefore,

m (ẍδx + ÿδy + z̈δz) = F_{x}δx + F_{y}δy + F_{z}δz  (15)

The above equation is known as D’Alembert’s equation.

If we have a system of particles acted upon by forces F_{1}, …, F_{p} and if each of the particles undergoes displacements δs_{1}, …, δs_{p}, then the total work done by the forces is

δW = ∑_{i=1}^{p} (F_{x_{i}}δx_{i} + F_{y_{i}}δy_{i} + F_{z_{i}}δz_{i}) .  (16)

Substituting (14) into D’Alembert’s equation and rearranging, we obtain

δW = m (ẍ\frac{∂x}{∂q_{1}} + ÿ\frac{∂y}{∂q_{1}} + z̈\frac{∂z}{∂q_{1}} ) δq_{1} + m (ẍ\frac{∂x}{∂q_{2}} + ÿ\frac{∂y}{∂q_{2}} + z̈\frac{∂z}{∂q_{2}} ) δq_{2}
= (F_{x} \frac{∂x}{∂q_{1}} + F_{y} \frac{∂y}{∂q_{1}} + F_{z} \frac{∂z}{∂q_{1}} ) δq_{1} + (F_{x} \frac{∂x}{∂q_{2}} + F_{y} \frac{∂y}{∂q_{2}} + F_{z} \frac{∂z}{∂q_{2}} ) δq_{2}. (17)

Let

F_{q_{1}} = F_{x} \frac{∂x}{∂q_{1}} + F_{y} \frac{∂y}{∂q_{1}} + F_{z} \frac{∂z}{∂q_{1}} ,

and

F_{q_{2}} = F_{x} \frac{∂x}{∂q_{2}} + F_{y} \frac{∂y}{∂q_{2}} + F_{z} \frac{∂z}{∂q_{2}} .

Then, (17) can be represented as,

δW = F_{q_{1}}δq_{1} + F_{q_{2}}δq_{2}.

To proceed further, note that

\frac{d}{dt} (ẋ\frac{∂x}{∂q_{1}} ) = ẍ\frac{∂x}{∂q_{1}} + ẋ\frac{d}{dt} (\frac{∂x}{∂q_{1}} ) .

Hence

ẍ\frac{∂x}{∂q_{1}} = \frac{d}{dt} (ẋ\frac{∂x}{∂q_{1}} ) − ẋ\frac{d}{dt} (\frac{∂x}{∂q_{1}} )  (18)

The time derivative of x = x(q_{1}, q_{2}) is

ẋ = \frac{∂x}{∂q_{1}} q̇_{1} + \frac{∂x}{∂q_{2}} q̇_{2}.  (19)

Differentiating the above partially with respect to q̇_{1} yields

\frac{∂ẋ}{∂q̇_{1}} = \frac{∂x}{∂q_{1}}  (20)

Finally, we will show that

\frac{d}{dt} (\frac{∂x}{∂q_{1}} ) = \frac{∂ẋ}{∂q_{1}}  (21)

Indeed, because x = x(q_{1}, q_{2}) and the partial derivative ∂x/∂q_{1} is in general a function of q_{1} and q_{2}, we have

\frac{d}{dt} (\frac{∂x}{∂q_{1}} ) = \frac{∂}{∂q_{1}} (\frac{∂x}{∂q_{1}} ) q̇_{1} + \frac{∂}{∂q_{2}} (\frac{∂x}{∂q_{1}} ) q̇_{2}.  (22)

Taking the partial derivative of (19) with respect to q_{1}, we get

\frac{∂ẋ}{∂q_{1}} = \frac{∂}{∂q_{1}} (\frac{∂x}{∂q_{1}} ) q̇_{1} + \frac{∂}{∂q_{1}} (\frac{∂x}{∂q_{2}} ) q̇_{2}
= \frac{∂}{∂q_{1}} (\frac{∂x}{∂q_{1}} ) q̇_{1} + \frac{∂}{∂q_{2}} (\frac{∂x}{∂q_{1}} ) q̇_{2},  (23)

because by assumption x, y, z ∈ C^{2} and hence

\frac{∂}{∂q_{1}} (\frac{∂x}{∂q_{2}} ) = \frac{∂}{∂q_{2}} (\frac{∂x}{∂q_{1}} ) .

With the above in mind, comparing the right hand sides of (22) and (23) yields (21). We now use (18), (20), and (21) to obtain the Lagrange equation of motion for a single particle. We first substitute (20) and (21) into (18) to get

ẍ\frac{∂x}{∂q_{1}} = \frac{d}{dt} (ẋ\frac{∂ẋ}{∂q̇_{1}} ) − ẋ\frac{∂ẋ}{∂q_{1}} .  (24)

Note that

\frac{∂}{∂q̇_{1}} (\frac{ẋ^{2}}{2} ) = ẋ\frac{∂ẋ}{∂q̇_{1}} ,  (25)

and

\frac{∂}{∂q_{1}} (\frac{ẋ^{2}}{2} ) = ẋ\frac{∂ẋ}{∂q_{1}} .  (26)

Substituting (25) and (26) into (24) gives

ẍ\frac{∂x}{∂q_{1}} = \frac{d}{dt} (\frac{∂ (\frac{ẋ^{2}}{2} )}{∂q̇_{1}} ) − \frac{∂ (\frac{ẋ^{2}}{2} )}{∂q_{1}}.

We can obtain similar expressions for y and z. Taking this into account, the expression (17) when δq_{2} = 0 becomes,

δW = (\frac{d}{dt} (\frac{∂}{∂q̇_{1}} \frac{m (ẋ^{2} + ẏ^{2} + ż^{2})}{2}) − \frac{∂}{∂q_{1}} \frac{m (ẋ^{2} + ẏ^{2} + ż^{2})}{2}) δq_{1}
= (F_{x} \frac{∂x}{∂q_{1}} + F_{y} \frac{∂y}{∂q_{1}} + F_{z} \frac{∂z}{∂q_{1}} ) δq_{1}
= F_{q_{1}}δq_{1}.  (27)

Let K = \frac{1}{2}m (ẋ^{2} + ẏ^{2} + ż^{2}) denote the kinetic energy of the particle. Then, we represent (27) as

\frac{d}{dt} (\frac{∂K}{∂q̇_{1}} ) − \frac{∂K}{∂q_{1}} = F_{q_{1}}  (28)

This equation is called the Lagrange equation of motion for the q_{1} coordinate. Using the same arguments as above, we can derive the Lagrange equation of motion for the q_{2} coordinate. In general there are as many the Lagrange equations of motion as there are degrees of freedom of the particle.

Recall that the expression δW = F_{x}δx + F_{y}δy + F_{z}δz is a work done by a force F = [ F_{x} F_{y} F_{z} ]^{⊤} for a general displacement δs = [ δx δy δz ]^{⊤} . In terms of generalized coordinates and when, say, δq_{2} = 0, we have

δW = (F_{x} \frac{∂x}{∂q_{1}} + F_{y} \frac{∂y}{∂q_{1}} + F_{z} \frac{∂z}{∂q_{1}} ) δq_{1}
= F_{q_{1}}δq_{1}
= δW_{q_{1}}.

A generalized force F_{q_{r}} is of such nature that the product F_{q_{r}}δq_{r} is the work done by driving forces when q_{r} alone is changed by δq_{r}. We mention that a generalized force does not have to be a force in the usual sense. For example, if q_{r} is an angle, then F_{q_{r}} must be a torque in order that F_{q_{r}}δq_{r} be work.

We now derive a more general version of the Lagrange equations of motion for a single particle. In our derivation we will need the notion of nonconservative forces that we discuss next. First recall that we interpret the kinetic energy of a particle as its ability to do work by virtue of its motion. This motion is the result of conservative and nonconservative forces acting upon the particle. A force is nonconservative if the work done by the force on a particle that moves through any round trip is not zero. Thus a force is nonconservative if the work done by that force on a particle that moves between two points depends on the path taken between these points.

Let us now suppose that in addition to the conservative forces, a single nonconservative force due to friction acts on the particle. Denote by ∑ W_{c} the sum of the work done by the conservative forces, and by W_{f} the work done by friction. We then have

W_{f} + ∑ W_{c} = ΔK.  (29)

Thus, if a nonconservative force, like friction, acts on a particle, then the total mechanical energy is not constant, but changes by the amount of work done by a nonconservative force according to the relation (29). Let us now denote by E the final mechanical energy of the particle and by E_{0} its initial mechanical energy. Then, we can write

E − E_{0} = ΔE = W_{f}.

The work done by friction, the nonconservative force, is always negative. Hence, the final mechanical energy E is less than the initial mechanical energy. The “lost” mechanical energy is transformed into internal energy U_{int}. In other words, the loss of mechanical energy equals the gain in internal energy and we write

ΔE + U_{int} = 0.

In general, we have

ΔK + ∑ ΔU + U_{int} + (change in other forms of energy) = 0.

This means that the total energy, kinetic plus potential plus internal plus all other forms of energy, does not change. This is the principle of the conservation of energy. To proceed further, recall that any generalized force F_{q_{i}} acting on our particle can be expressed as

F_{q_{i}} = F_{x} \frac{∂x}{∂q_{i}} + F_{y} \frac{∂y}{∂q_{i}} + F_{z} \frac{∂z}{∂q_{i}} .

Assuming that the forces are conservative, we get

F_{q_{i}} = − (\frac{∂U}{∂x} \frac{∂x}{∂q_{i}} + \frac{∂U}{∂y} \frac{∂y}{∂q_{i}} + \frac{∂U}{∂z} \frac{∂z}{∂q_{i}} ) .

We thus have

F_{q_{i}} = −\frac{∂U}{∂q_{i}} .

Therefore, the Lagrange equations of motion for a single particle for conservative forces can be represented as

\frac{d}{dt} (\frac{∂K}{∂q̇_{i}} ) − \frac{∂K}{∂q_{i}} = −\frac{∂U}{∂q_{i}} .

We can rearrange the above equations as

\frac{d}{dt} (\frac{∂K}{∂q̇_{i}} ) − \frac{∂}{∂q_{i}} (K − U) = 0.

We now introduce the Lagrangian function L, defined as

L = K − U.

Note that because U is not a function of q̇_{i}, we have

\frac{∂K}{∂q̇_{i}} = \frac{∂L}{∂q̇_{i}}.

Taking the above into account, we can represent the Lagrange equations of motion for conservative forces as

\frac{d}{dt} (\frac{∂L}{∂q̇_{i}}) − \frac{∂L}{∂q_{i}} = 0.

If some of the forces acting on the particle are nonconservative, the Lagrange equations of motion for our particle will take the form

\frac{d}{dt} (\frac{∂L}{∂q̇_{i}}) − \frac{∂L}{∂q_{i}} = F̃_{q_{i}}  (30)

where F̃_{q_{i}} are non conservative forces acting on the particle.

Example 1 Consider the simple pendulum shown in Figure 2. The simple pendulum is an idealized body consisting of a point mass M, suspended by a weightless inextensible cord of length l. The simple pendulum is an example of a one degree of freedom system with a generalized coordinate being the angular displacement θ. The pendulum kinetic energy is

K = \frac{1}{2}Ml^{2}θ̇^{2}.

We assume that the potential energy of the pendulum is zero when the pendulum is at rest. Hence, its potential energy is

U = Mgl(1 − cos θ),

Figure 2: A simple pendulum.

where g = 10 m/sec^{2} is the acceleration due to gravity. The Lagrangian function is

L = K − U
= \frac{1}{2}Ml^{2}θ̇^{2} − Mgl(1 − cos θ)
= 6θ̇^{2} − 60(1 − cos θ).

The Lagrange equation describing the pendulum motion is

\frac{d}{dt} (\frac{∂L}{∂θ̇}) − \frac{∂L}{∂θ} = 0.

In our example the above equation evaluates to

12θ̈ + 60 sin θ = 0,

or equivalently

θ̈ = −5 sin θ.

4 Double Pendulum Modeling

We construct a model of the double pendulum shown in Figure 3. This double pendulum

Figure 3: Double pendulum.

Table 1: Parameter numerical values of double pendulum system shown in Figure 3.

m_{1} = m_{2} 1 kg
l_{1} = l_{2}  0.5 m
g  9.8 m/sec^{2}

can be viewed as a model of a double-link robot manipulator. We assume that the input torques, T_{1} and T_{2} are applied at each joint, ω_{i} = θ̇_{i}, i = 1, 2, are the angular velocities of the pendulum arms, Δθ = θ_{2} − θ_{1}, and m_{i} and l_{i} are the masses and arm lengths, respectively. The parameter values are given in Table 1. We write the Lagrange equations of motion for the system and represent them in state-space format, where the state variables are defined as, x_{1} = θ_{1}, x_{2} = θ̇_{1}, x_{3} = θ_{2}, and x_{4} = θ̇_{2}. Lagrange Equations of Motion of the Double Pendulum System

Recall that the Lagrangian L is defined as L = K − U, where K is the kinetic energy of the double pendulum system shown in Figure 3 and U is the potential energy of the system. The kinetic energy is

K = K_{m_{1}} + K_{m_{2}}
= \frac{1}{2}m_{1}∥v_{m_{1}}∥^{2} + \frac{1}{2}m_{1}∥v_{m_{2}}∥^{2},

where v_{m_{i}}, i = 1, 2, are the velocities of the point masses. The velocity vectors can be computed from the position of the point masses. The position of the first mass is described by

x_{1} = l_{1} sin(θ_{1}), y_{1} = l_{1} cos(θ_{1}).  (31)

The position of m_{2} is determined by

x_{2} = l_{1} sin(θ_{1}) + l_{2} sin(θ_{2}), y_{2} = l_{1} cos(θ_{1}) + l_{2} cos(θ_{2}).  (32)

We next find time derivatives of the expressions in (31) and (32). The time derivatives of x_{1} and y_{1} are,

ẋ_{1} = l_{1}θ̇_{1} cos(θ_{1}), ẏ_{1} = −l_{1}θ̇_{1} sin(θ_{1}).

The time derivatives of x_{2} and y_{2} are,

ẋ_{2} = l_{1}θ̇_{1} cos(θ_{1}) + l_{2}θ̇_{2} cos(θ_{2}), ẏ_{2} = −l_{1}θ̇_{1} sin(θ_{1}) − l_{2}θ̇_{2} sin(θ_{2}).

Using the above expressions, we compute

∥v_{m_{1}}∥^{2} = ẋ_{1}^{2} + ẏ_{1}^{2} = (l_{1}θ̇_{1})^{2}
∥v_{m_{2}}∥^{2} = ẋ_{1}^{2} + ẏ_{1}^{2}
= (l_{1}θ̇_{1})^{2} cos^{2}(θ_{1}) + 2l_{1}l_{2}θ̇_{1}θ̇_{2} cos(θ_{1}) cos(θ_{1}) + (l_{2}θ̇_{2})^{2} cos^{2}(θ_{2})
+ (l_{1}θ̇_{1})^{2} sin^{2}(θ_{1}) + 2l_{1}l_{2}θ̇_{1}θ̇_{2} sin(θ_{1}) sin(θ_{1}) + (l_{2}θ̇_{2})^{2} sin^{2}(θ_{2})
= (l_{1}θ̇_{1})^{2} + (l_{2}θ̇_{2})^{2} + 2l_{1}l_{2}θ̇_{1}θ̇_{2} cos(θ_{2} − θ_{1}),

where we used the well known trigonometric identity,

cos(θ_{1} − θ_{1}) = cos(θ_{2}) cos(θ_{1}) + sin(θ_{2}) sin(θ_{1}).

In view of the above, the kinetic energy of the system is

K = K_{m_{1}} + K_{m_{2}}
= \frac{1}{2}m_{1}∥v_{m_{1}}∥^{2} + \frac{1}{2}m_{1}∥v_{m_{2}}∥^{2}
= \frac{1}{2}m_{1} (l_{1}θ̇_{1})^{2} + \frac{1}{2}m_{2} ((l_{1}θ̇_{1})^{2} + (l_{2}θ̇_{2})^{2} + 2l_{1}l_{2}θ̇_{1}θ̇_{2} cos(θ_{2} − θ_{1})) .

To compute the potential energy, we observe the positive direction of the y-coordinate and the direction of the gravity force are the same. Therefore, the potential energy is

U = −m_{1}gy_{1} − m_{2}gy_{2}
= −m_{1}gl_{1} cos(θ_{1}) − m_{2}g (l_{1} cos(θ_{1}) + l_{2} cos(θ))
= −(m_{1} + m_{2})gl_{1} cos(θ_{1}) − m_{2}gl_{2} cos(θ_{2}).

Figure 4: Angles in the double pendulum system of Figure 3.

Hence, the Lagrangian function of the pendulum system is

L = K − U
= \frac{1}{2}m_{1} (l_{1}θ̇_{1})^{2} + \frac{1}{2}m_{2} ((l_{1}θ̇_{1})^{2} + (l_{2}θ̇_{2})^{2} + 2l_{1}l_{2}θ̇_{1}θ̇_{2} cos(θ_{2} − θ_{1}))
+(m_{1} + m_{2})gl_{1} cos(θ_{1}) + m_{2}gl_{2} cos(θ_{2}).

We next derive the Lagrange equations of motion of the system. For this, we need to determine the generalized forces that constitute the right hand sides of the Lagrange equations of motion. Referring to Figure 4, we write,

δW = T_{1}δθ_{1} + T_{2}δ(θ_{2} − θ_{1})
= (T_{1} − T_{2})δθ_{1} + T_{2}δθ_{2},  (33)

where T_{1} and T_{2} are the torques applied at each joint. Hence the generalized forces are

F_{θ_{1}} = T_{1} − T_{2} and F_{θ_{2}} = T_{2}.  (34)

To proceed, we need to find partial derivatives of the Lagrangian function. We first evaluate ∂L/∂θ̇_{1} to get

\frac{∂L}{∂θ̇_{1}} = m_{1}l_{1}^{2}θ̇_{1} + m_{2}l_{1}^{2}θ̇_{1} + m_{2}l_{1}l_{2}θ̇_{2} cos(θ_{2} − θ_{1}).  (35)

Next, we find

\frac{d}{dt} \frac{∂L}{∂θ̇_{1}} = m_{1}l_{1}^{2}θ̈_{1} +m_{2}l_{1}^{2}θ̈_{1} +m_{2}l_{1}l_{2}θ̈_{2} cos(θ_{2} −θ_{1})−m_{2}l_{1}l_{2}θ̇_{2}(dottheta_{2} −θ̇_{1}) sin(θ_{2} −θ_{1}). (36)

The partial derivative ∂L/∂θ_{1} is

\frac{∂L}{∂θ_{1}} = m_{2}l_{1}l_{2}θ̇_{1}θ̇_{2} sin(θ_{2} − θ_{1}) − (m_{1} + m_{2})gl_{1} sin(θ_{1}).  (37)

Combining (34), (36), and (37), we obtain the first Lagrange equation of motion,

(m_{1}+m_{2})l_{1}^{2}θ̈_{1}+m_{2}l_{1}l_{2}θ̈cos(θ_{2}−θ_{1})−m_{2}l_{1}l_{2}θ̇_{2}^{2} sin(θ_{1}−θ_{1})+(m_{1}+m_{2})gl_{1} sin(θ_{1}) = T_{1}−T_{2}. (38)

We next compute partial derivatives of L that we use to write the second Lagrange equation of motion. We first find ∂L/∂θ̇_{2},

\frac{∂L}{∂θ̇_{2}} = m_{2}l_{2}^{2}θ̇_{2} + m_{2}l_{2}l_{2}θ̇_{1} cos(θ_{2} − θ_{1}).  (39)

The time derivative of the above expression is

\frac{d}{dt} \frac{∂L}{∂θ̇_{2}} = m_{2}l_{2}^{2}θ̈_{2} + m_{2}l_{1}l_{2}θ̈_{1} cos(θ_{2} − θ_{1}) − m_{2}l_{1}l_{2}(θ̇_{2} − θ̇_{1})θ̇_{1} sin(θ_{1} − θ_{1}).  (40)

The partial derivative ∂L/∂θ_{2} is

\frac{∂L}{∂θ_{2}} = −m_{2}l_{1}l_{2}θ̇_{1}θ̇_{2} sin(θ_{2} − θ_{1}) − m_{2}gl_{2} sin(θ_{2}).  (41)

Combining (34), (40), and (41), we obtain the second Lagrange equation of motion,

m_{2}l_{2}^{2}θ̈_{2} + m_{2}l_{1}l_{2}θ̈_{1} cos(θ_{2} − θ_{1}) + m_{2}l_{1}l_{2}θ̇_{1}^{2} sin(θ_{2} − θ_{1}) + m_{2}gl_{2} sin(θ_{2}) = T_{2}.  (42)

State-Space Model of the Double Pendulum System Taking into account the definition of the state variables, we obtain

[ ẋ_{1} ; ẋ_{2} ; ẋ_{3} ; ẋ_{4} ] = [ x_{2} ; θ̈_{1} ; x_{4} ; θ̈_{2} ] .  (43)

Thus, to obtain a state-space model of the double pendulum system, we need to find the acceleration of each link. For this, we represent Lagrange equations (38) and (42) in matrix form,

[ \begin{matrix}(m_{1} + m_{2})l_{1}^{2} & m_{2}l_{1}l_{2} cos(θ_{2} − θ_{1}) \\ m_{2}l_{1}l_{2} cos(θ_{2} − θ_{1}) & m_{2}l_{2}^{2} \end{matrix} ] [ θ̈_{1} ; θ̈_{2} ]
= [ \begin{matrix}T_{1} − T_{2} + m_{2}l_{1}l_{2}θ̇_{2}^{2} sin(θ_{1} − θ_{1}) − (m_{1} + m_{2})gl_{1} sin(θ_{1}) \\ T_{2} − m_{2}l_{1}l_{2}θ̇_{1}^{2} sin(θ_{2} − θ_{1}) − m_{2}gl_{2} sin(θ_{2}) \end{matrix} ] . (44)

Note that the square matrix, call it M, on the left-hand side of the above equation is invertible because its determinant is

det M = (m_{1} + m_{2})m_{2}l_{1}^{2}l_{2}^{2} − m_{2}^{2}l_{1}^{2}l_{2}^{2} cos^{2}(θ_{2} − θ_{1}) > 0.

Hence we can solve (44) for θ̈_{1} and θ̈_{2}. We have

θ̈_{1} = \frac{m_{2}l_{2}^{2} (T_{1} − T_{2} + m_{2}l_{1}l_{2}θ̇_{2}^{2} sin(θ_{2} − θ_{1}) − (m_{1} + m_{2})gl_{1} sin(θ_{1})) −m_{2}l_{1}l_{2} cos(θ_{2} − θ_{1}) (T_{2} − m_{2}l_{1}l_{2}θ̇_{1}^{2} sin(θ_{2} − θ_{1}) − m_{2}gl_{2} sin(θ_{2}))}{(m_{1} + m_{2})m_{2}l_{1}^{2}l_{2}^{2} − m_{2}^{2}l_{1}^{2}l_{2}^{2} cos^{2}(θ_{2} − θ_{1})}  (45)

and

θ̈_{2} = \frac{−m_{2}l_{1}l_{2} cos(θ_{2} − θ_{1}) (T_{1} − T_{2} + m_{2}l_{1}l_{2}θ̇_{2}^{2} sin(θ_{2} − θ_{1}) − (m_{1} + m_{2})gl_{1} sin(θ_{1})) +(m_{1} + m_{2})l_{1}^{2} (T_{2} − m_{2}l_{1}l_{2}θ̇_{1}^{2} sin(θ_{2} − θ_{1}) − m_{2}gl_{2} sin(θ_{2}))}{(m_{1} + m_{2})m_{2}l_{1}^{2}l_{2}^{2} − m_{2}^{2}l_{1}^{2}l_{2}^{2} cos^{2}(θ_{2} − θ_{1})} . (46)

The input vector is defined as

u = [ u_{1} u_{2} ]^{⊤} = [ T_{1} T_{2} ]^{⊤} .

Taking into account the above definition of the input, the definition of the state variables given by (43) and relations (45) and (46), we obtain the following state-space model of the double pendulum system

[ ẋ_{1} ; ẋ_{2} ; ẋ_{3} ; ẋ_{4} ] = f(x, u)

= [ x_{2} \frac{m_{2}l_{2}^{2} (u_{1} − u_{2} + m_{2}l_{1}l_{2}x_{4}^{2} sin(x_{3} − x_{1}) − (m_{1} + m_{2})gl_{1} sin(x_{1})) −m_{2}l_{1}l_{2} cos(x_{3} − x_{1}) (u_{2} − m_{2}l_{1}l_{2}x_{2}^{2} sin(x_{3} − x_{1}) − m_{2}gl_{2} sin(x_{3}))}{(m_{1}+m_{2})m_{2}l_{1}^{2}l_{2}^{2}−m_{2}^{2}l_{1}^{2}l_{2}^{2} cos^{2}(x_{3}−x_{1})} x_{4} \frac{−m_{2}l_{1}l_{2} cos(x_{3} − x_{1}) (u_{1} − u_{2} + m_{2}l_{1}l_{2}x_{4}^{2} sin(x_{3} − x_{1}) − (m_{1} + m_{2})gl_{1} sin(x_{1})) +(m_{1} + m_{2})l_{1}^{2} (u_{2} − m_{2}l_{1}l_{2}x_{2}^{2} sin(x_{3} − x_{1}) − m_{2}gl_{2} sin(x_{3}))}{(m_{1}+m_{2})m_{2}l_{1}^{2}l_{2}^{2}−m_{2}^{2}l_{1}^{2}l_{2}^{2} cos^{2}(x_{3}−x_{1})} ] .

References

[1] K. D. Elworthy, W. N. Everitt, and E. B. Lee, editors. Differential Equations, Dynamical Systems, and Control Science: A Festschrift in Honor of Lawrence Markus. Marcel Dekker, Inc., New York, 1994. Lecture Notes in Pure and Applied Mathematics, Vol. 152.

[2] J. R. Leigh. Modelling and Simulation. IEE Topics in Control Series 1. Peter Peregrinus Ltd., London, 1983.

[3] Reader’s Digest Practical Problem Solver. Substitutes, Shortcuts, and Ingenious Solutions for Making Life Easier. The Reader’s Digest Association, Inc., Pleasantville, New York, 1991.