A comprehensive dynamic model for class-1 tensegrity systems based on quaternions

Author's personal copy A comprehensive dynamic model for class-1 tensegrity systems based on quaternions Massimo Cefalo a, ⇑ , Josep M. Mirats-Tur b,1 a Institut de Robòtica i Informàtica Industrial, LLorens i Artigas 4-6 2nd Floor, 08028 Barcelona, Spain b Cetaqua, Passeig dels Til.lers, 3, 08034 Barcelona, Spain article info Article history: Received 3 March 2010 Received in revised form 26 September 2010 Available online 16 November 2010 Keywords: Dynamic model Tensegrity Quaternions abstract In this paper we propose a new dynamic model, based on quaternions, for tensegrity systems of class-1. Quaternions are used to represent orientations of a rigid body in the 3-dimensional space eliminating the problem of singularities. Moreover, the equations based on quaternions allow to perform more precise calculations and simulations because they do not use trigonometric functions for the representation of angles. We present a thorough introduction of tensegrities and the current state of research. We also introduce the quaternions and provide in the appendix some important details and useful properties. Applying the Euler–Lagrange approach we derive a comprehensive dynamic model, first for a simple rigid bar in the space and, at last, for a class-1 tensegrity system. We present two model forms: a matrix and a vectorial form. The first more compact and easier to write, the latter more suitable to apply the tools and the theory based on vector fields. Ó2010 Elsevier Ltd. All rights reserved. 1. Introduction Tensegrities are mechanical systems born in the art community in the early 60s (Snelson, 1965), at the beginning studied by architects and successively widely used in engineering applications. Among several definitions we can say that a tensegrity system is an aggregation of mechanical elements, carrying either compression or tension but not both, for which at least one equilibrium state exists. Mixing together the terms tensile and integrity, which emphasize the main characteristic of such structures, the name tensegrity was coined by Buckminister Fuller in the early 60s (Fuller, 1962). An interesting and quite general definition was given by Pugh (1976):‘‘A tensegrity system is established when a set of discontinuous compressive components interacts with a set of continuous tensile components to define a stable volume in space’’. Elements working under compression are called struts while the elements working in tension are called tensors. Usually the struts are rigid bars and the tensors are cables or springs, but several variations may be allowed. Actuating bars and cables allow to dynamically modify their lengths and, hence, change the configuration of the system. To build a tensegrity structure, several approaches have been proposed. Between others, one of the most powerful techniques is based on the construction of a base module and of a design pattern. The overall structure is the result of assembling the base module along the design pattern following proper rules to ensure the existence of, at least, one equilibrium point. So, usually, the most famous tensegrity structures have some degree of geometric symmetry in their equilibrium positions. These systems have a lot of interesting properties. It has been shown, for example, that carefully designing the appropriate net of connections between rigid elements, it is possible to provide to the structure the desired rigidity (within the limits of the material employed). Compared with traditional structures, tensegrities can be much more stiff, much more light and occupy much less space and volume. The main applications in the robotic field, tend to exploit the capability of such systems to be extensible and redundant (in case of one, or more, struts or tensors fail, it is possible to guarantee the functionality of the system by means of other elements). See, for instance, Aldrich (2004), Paul et al. (2006), Masic and Skelton (2004) and Mirats-Tur (2010) for some application concerning manipulators and mobile robotics. For several years they have been studied only from a static point of view, mainly because there were no possibilities to face with the enormous quantity of computations required to model their dynamics. The static analysis of tensegrity has reached a certain level of maturity, with lots of contributions by different authors and fields of study. See for instance Roth and Whiteley (1981), Connelly (1999) for a mathematic perspective, Calladine and Pellegrino (1992), Motro et al. (1986), Hanaor (1988) for studies about the self-stress states of the structure, or Tarnai (1989), Vassart et al. 0020-7683/$ - see front matter Ó2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijsolstr.2010.11.015 ⇑ Corresponding author. Tel.: +39 0649766812; fax: +39 064957647. E-mail addresses: [email protected] (M. Cefalo), [email protected] (J.M. Mirats-Tur). 1 Tel.: +34 933124879; fax: +34 933124801. International Journal of Solids and Structures 48 (2011) 785–802 Contents lists available at ScienceDirect International Journal of Solids and Structures journal homepage: www.elsevier.com/locate/ijsolstr Author's personal copy (2000) talking about the prestressability problem. A recent review about tensegrity statics was given in Hernandez and Mirats-Tur (2008) presenting different existing definitions for tensegrity structures, as well as their main properties. The dynamics were first studied by Motro et al. (1986). Kanchanasaratool and Williamson (2002) studied dynamic particle models while considering the bars to be massless; other studies, as those by Skelton et al. (2001) or Sultan (1999) consider mass on bars. Also non-linear models and their linearizations have been considered by Sultan et al. (2002). A recent survey on the dynamics of tensegrity structures including current open research problems was given in Mirats-Tur and Hernández (2009). In the following we shall focus on the so called class-1 tensegrity structures. A definition of class-ktensegrity structures can be found in Skelton and Oliveira (2009):a tensegrity system with k rigid components in contact, with a frictionless ball joint, at a given node. So a class-1 tensegrity is a structure in which there is only one rigid element per node, or in other words (as captured in Pugh definition): systems where all tensors constitute a continuous set of elements and no rigid bodies are in contact. This is the simplest category, yet useful enough to allow the study of general dynamical models and their properties. To obtain a set of ordinary differential equations of a complex system three different methods are commonly used: Newton– Euler, Lagrange and Hamilton approaches. The Newton–Euler approach is based on the two cardinal equations of the classical mechanics and is particularly well suited for those situations in which it is possible to exploit the recursive structure of the system so that the equations can be written in algorithmic way. It is the case, for example, of robot manipulators. For more details on this approach read Arnold (1978) and Siciliano et al. (2009) for applications to robotics. In general, for complex systems, Lagrange and Hamilton approaches allow writing the equations of motion in a more straightforward way. Both of them are based on systematic equations which directly lead to the closed form equations of motion once the energy of the system has been determined. These methods are a direct consequence of a variational principle (Hamilton’s principle), which states, in a nutshell, that the trajectories of a system can only be solutions when they are stationary points of a special integral function defined by the energy of the system (named the action integral). For the systems studied here, these stationary points should be minimums of this function, involving boundary conditions not directly tested with the stationary point of the action integral. We propose to explore the Lagrange approach for a general system composed of nrigid bodies that are linked between them by means of massless connections. Such connections allow forces to be exchanged between bodies. To represent the dynamics of the system without singularities we will introduce and make use of the quaternions. This paper is organized as follows: We begin by a section in which an algorithm to associate a quaternion to a bar starting from a vector which locates one endpoint of the bar with respect to the other (the bar vector) is presented. In order to maintain the paper self-contained, further details and properties about quaternions exploited to obtain the dynamic model can be found in the Appendix 1. In Section 3we apply the Euler–Lagrange equations to a class-1 tensegrity structure using quaternions. Section 3.1 introduces the notation utilized and the framework of reference; in the Sections 3.2 and 3.3 we study, respectively, the kinetic and the potential energy. Finally, in Section 3.4 we deploy the equations of motion, first for a single bar and then for a complete tensegrity structure, differentiating two forms: a vectorial and a matrix form (some calculations and useful properties are collected in Appendix 2). In Section 4we present some simulation results aimed to show possible applications of the model proposed and to provide an empirical proof of its correctness. Last, Section 5 contains some conclusions and hints about possible future developments. 2. Three-dimensional rotations In this paragraph we will focus on unit quaternions, with the aim to give an example of how they can be exploited to define the orientation of a rigid body. Please refer to Appendix 1 for further details on quaternions and used notation. A unit quaternion ecan be written as follows: e¼ e 0 e 1 e 2 e 3 0 B B B @1 C C C A¼cosðh=2Þ sinðh=2Þu  ð1Þ where urepresents a versor (unit vector). From Eq. (1), it is easy to see that a unit quaternion contains information about a versor and an angle, which may be thought as a rotation angle around that versor (h2[ p , p ]). Roughly speaking, a unit quaternion specifies an axis and the rotation around it, hence we can also say that it represents a new reference framework once a base reference frame (with respect to which the vector uis expressed) is fixed. The framework it represents is the result of a rotation of the base frame around the uvector of the hangle. Let now take into account a rigid body in the space. Say, for simplicity, and to anticipate the discussions of the next sections concerning tensegrity structures, this body is a simple bar. We want to represent the orientation of a bar in the space by means of a quaternion. We will neglect rotations of the bar around its axis (so we finally have five degrees of freedom in a 5-dimensional coordinate space). To associate a quaternion to a bar means associating it to a relative reference frame: a problem with infinite solutions. Hence, there are infinite ways to associate a quaternion to a bar, while there is only one solution to the reverse problem (a quaternion uniquely identify one orientation). Let r 0 be a vector quaternion represented in the inertial reference frame Ox 0 y 0 z 0 (the base frame from now on), and let r 1 be the representation of the same vector quaternion in another frame Ox 1 y 1 z 1 . Let ebe the quaternion representing the orientation of the frame 1 with respect to the base frame, then, it is known that r 0 ¼er 1 e  ð2Þ The quaternion e*represents the orientation of the base frame in the frame Ox 1 y 1 z 1 . Let P 1 and P 2 be the coordinates of the extremes of the bar expressed in the base frame. The bar geometric vector, in the base frame, is defined by: ~ b¼~ p 2 ~ p 1 ð3Þ where ~ p 1 and ~ p 2 are the geometric vectors associated to the points P 1 and P 2 . The quaternion associated to the bar in the base frame is: b 0 ¼0 b  ¼ 0 b 1 b 2 b 3 0 B B B @1 C C C Að4Þ where bis the vector associated to the bar meant as element of a vector space (see Fig. 1). Suppose now to associate a frame to the bar with the origin in the center and one of the main axis directed like the bar; let it be, 786 M. Cefalo, J.M. Mirats-Tur / International Journal of Solids and Structures 48 (2011) 785–802 Author's personal copy for example, the xaxis. In this frame the vector quaternion associated to the bar is: b 1 ¼ 0 l 0 0 0 B B B @1 C C C Að5Þ where lis length of the bar: kp 2 p 1 k. From Eq. (2), we can write: b 0 ¼eb 1 e  ¼E þ b þ1 e  ¼E þ E T b 1 ¼E T E þ b 1 ð6Þ The matrices E þ and Eare defined like Q þ and Qin (A.1). Eq. (6), in matrix form becomes: b 0 ¼ 0 b 1 b 2 b 3 0 B B B @1 C C C A¼ 1000 0 0 0 E þ E T 0 B B B @1 C C C A 0 l 0 0 0 B B B @1 C C C Að7Þ where E þ ¼e 1:3 ðe 0 I 3 þe EÞ hi ¼e 1 e 0 e 3 e 2 e 2 e 3 e 0 e 1 e 3 e 2 e 1 e 0 0 B @1 C Að8Þ E T ¼e 1:3 ðe 0 I 3 e EÞ hi T ¼ e 1 e 2 e 3 e 0 e 3 e 2 e 3 e 0 e 1 e 2 e 1 e 0 0 B B B @1 C C C Að9Þ are two special matrices associated to the quaternion that represents the bar orientation. Observe that E þT and E T can be obtained, respectively, from Eand E þ suppressing the first column: EE þ ## E þT E T ð10Þ From Eq. (7) we can write b 1 b 2 b 3 0 B @1 C A¼E þ E T l 0 0 0 B @1 C Að11Þ which is equivalent to the following non-linear system: b 1 ¼lðe 2 0 þe 2 1 e 2 2 e 2 3 Þ b 2 ¼2lðe 0 e 3 þe 1 e 2 Þ b 3 ¼2lðe 1 e 3 e 0 e 2 Þ 8 > < > :ð12Þ System (12) allows to determine the components of bgiven a quaternion e. To find a solution of the reverse problem, e=f(b), we must also add the constraint equation e 2 0 þe 2 1 þe 2 2 þe 2 3 ¼1ð13Þ The resulting system is made of four equations in four unknowns, is non-linear and generally non-invertible, unless some hypothesis are made (for the reasons above discussed the problem has infinite solutions). A way to find a solution is to put one component of the vector part of eto zero. This fact corresponds to fix the vector uof Eq. (1) in a plane defined by two main axis of the base frame. Adding the constraint u x =0oru y =0oru z = 0, we obtain a system of five equations in four unknowns and it may be proven that it always admits one solution, with the exception of one case. 2 Such case is a singularity of the general solution and corresponds to a special orientation of the bar. There is a different singularity for every choice of the additional constraint on u. For all singularities, the unit quaternion is however easily determinable: all components of ewill be zero, unless one, which will be 1 or 1. Let us now show the algorithm introduced with an example. Consider the system (12) with the constraint (13). Suppose to add the constraint equation e 1 = 0. It corresponds to put u x =0, i.e., to choose the vector u, which defines the direction around which rotate the base frame to obtain the frame associated to the bar, in the yz plane. Furthermore, suppose to fix the origin of the bar frame on the center of the bar and to orient the xaxis along the bar itself. The resulting system becomes: b 1 ¼le 2 0 e 2 2 e 2 3  b 2 ¼2le 0 e 3 b 3 ¼2le 0 e 2 1¼e 2 0 þe 2 2 þe 2 3 8 > > > > > < > > > > > : ð14Þ Fig. 2 shows the solution for a given orientation of a bar. There are infinite frames with the xaxis oriented along the bar and the origin coinciding with one of the extreme points. One of such frames corresponds to a rotation of the base frame around a vector in the yz plane (of the base plane). In the figure, the frame O 1 x 1 y 1 z 1 (in red) 3 is the frame associated to the red versor u. The analytical solution can be found as follows. From the first and the fourth of (14) we have 2e 2 0 1¼b 1 l which implies ð15Þ From the second we have e 3 ¼b 2 2le 0 Fig. 1. A single bar in the space. 2 Observe that a non-linear algebraic system, in general, may admit several solutions even if the number of unknowns is less than the number of equations. 3 For interpretation of color in Figs. 1-9,11,12, the reader is referred to the web version of this article. M. Cefalo, J.M. Mirats-Tur / International Journal of Solids and Structures 48 (2011) 785–802 787 Author's personal copy which, with Eq. (15), gives ð16Þ Finally, from the third we obtain e 2 ¼ b 3 2le 0 which, together with Eq. (15), gives ð17Þ Eqs. (15), (17), (16) with ð18Þ are a solution of the system made of (12) and of Eq. (13): it can be directly proven substituting the solution found into the equations. The solution found has a singularity for b 1 =l, i.e., for the case in which the bar is oriented in the opposite direction of the xaxis of the base frame. In such case, the orientation of the bar frame is represented by the following unit quaternion: e¼ 0 0 0 1 0 B B B @1 C C C A There is a singularity in any case, i.e., whichever is the component of the vector part of the unit quaternion that we choose to put equal to zero. In general, if we choose to align the xaxis of the bar fixed frame with the bar itself, depending on which component of e we put to zero, we have the following cases: e 1 ¼0)b¼l 0 0 0 B @1 C A e 2 ¼0;e 3 ¼0)b¼ l 0 0 0 B @1 C A The singularity depends on the way we choose to orient the bar fixed frame with respect to the bar and on which component, of the vector part of e, we put equal to zero. If we choose a bar fixed frame oriented with the yaxis directed along the bar, i.e. b 1 =(00l0) T the system (12) becomes: b 1 ¼2lðe 1 e 2 e 0 e 3 Þ b 2 ¼lðe 2 0 þe 2 2 e 2 1 e 2 3 Þ b 3 ¼2lðe 0 e 1 þe 2 e 3 Þ 8 > < > :ð19Þ and the singularities occur for: e 1 ¼0)b¼ 0 l 0 0 B @1 C A e 1 ¼0;e 3 ¼0)b¼ 0 l 0 0 B @1 C A 3. Euler–Lagrange equations for a tensegrity structure In this section we will focus on the so called class-1 tensegrity structures. We will deploy a dynamic model making use of quaternions to describe the orientation of the bars. Fig. 2. Correspondence between rotations around the x 1 axis of the frame O 1 x 1 y 1 z 1 and rotations around the versor uof the frame O 0 x 0 y 0 z 0 . 788 M. Cefalo, J.M. Mirats-Tur / International Journal of Solids and Structures 48 (2011) 785–802 Author's personal copy 3.1. The framework and the notations Consider a class-1 tensegrity system where the struts are rigid bars and the tensors are cables. We admit, for simplicity, that the bars have a radial dimension negligible with respect to their length. On the basis of this hypothesis, we will model the bars as straight lines in the space. For this purpose we will need three independent parameters to describe the position and two independent parameters to describe the orientation. Before introducing the generalized coordinates and the notation, let us insert here a brief digression concerning terminology and conventions. To avoid confusion and help the reader, we recall some definitions from geometry and algebra. With the term vector we mean an element of a vector space. To avoid confusion, the oriented segment joining two points in R 3 will be called geometric vector. The set of all geometric vectors, considered applied at the origin of a reference frame, together with the product of a geometric vector by a scalar and the sum of two geometric vectors is also a vector space. Often, it is easy to get confusion between the Euclidean space and the vector space of geometric vectors. Points of the Euclidean space and geometric vectors representing such points are different entities. Displacements in R 3 are geometric vectors. Usually, position vectors (and therefore displacements) are considered as applied vectors, while velocities and accelerations (both, angular or linear) are not. A position vector is said to be referred to a frame when it is meant as applied to the origin of such frame. The transformation applied by a rotation matrix to a vector changes the frame in which it is expressed, but not the frame which the vector refers to. Consider Fig. 3, the vector p 1 is the position vector associated to the point P, referred to the frame 1and expressed in the same frame. If R 0 1 is the rotation matrix representing the orientation of the frame 1with respect to the frame 0;R 0 1 p 1 represents the position vector of the point Pexpressed in the frame 0 but still referring the frame 1. Besides, p 0 ¼o 0 1 þR 0 1 p 1 is the vector position associated to the point Pexpressed in the frame 0 and referred in the same frame. The quantities p 0 and p 1 are meant as elements of the vector space of all geometric vectors of R 3 . This is the reason why we write them without the arrow. The two reference frames, meant coincident in the origin, are the geometric equivalent of a base change in the vector space. In Fig. 3 we associate to the geometric vectors (the arrows) their equivalent vectors meant as elements of a vector space. A superscript is typically used to put in evidence the frame in which a variable is expressed, but sometimes, to simplify the notation, we will omit the superscript to indicate that a variable is expressed in the base frame. Let us now turn talking about tensegrities. To describe the system we will start introducing the description of a single bar. Let Ox 0 y 0 z 0 be an inertial reference frame (the base frame) and let P 1 and P 2 be the coordinates of the extremes of the bar and ~ p 1 and ~ p 2 be the geometric vectors associated to the points P 1 and P 2 (see Fig. 1). Let us now associate a frame to the bar like described in Section 2(see also Fig. 2): the origin is in the center and the xaxis is directed along the bar itself. Let ebe the unitary quaternion representing the orientation of the bar-fixed frame (and hence also of the bar). The bar may be completely described by means of eand by the position of anyone of its points. We choose the center of mass P c (whose associated vector will be denoted with c) which can be easily computed starting from P 1 and P 2 : P c ¼1 2ðP 1 þP 2 Þ)c¼1 2ðp 1 þp 2 Þ On the base of what introduced for a single bar, we define in the following the notation needed to describe a class-1 tensegrity system. We will introduce a double notation, to be able to present two forms of the dynamic model: a matrix form (compact form) and a classical vector form, both based on the Lagrange equations. The former is useful because compact and elegant, the latter, allows to manage the model with vector fields. Therefore, for each quantity we introduce a matrix notation and a vector notation. For the positions of the nodes and the bars we define N¼ðN 1 N 2 Þ¼ðn 1 n 2n b Þnodes row matrix,32n b where n b is the number of the bars. Bars are defined as b i ¼n iþn b n i where n iþn b 2N 2 and n i 2N 1 m ¼vecðNÞ¼ m 1 m 2  ¼ n 1 . . . n 2n b 0 B @1 C A nodes column matrix,6n b 1 B¼ðb 1 b n b Þ¼NI n b I n b  bar row matrix,3n b b¼vecðBÞ¼ðI 3n b I 3n b Þ m bar column matrix,3n b 1 In the above definitions, the symbol vec represents the operator which put in one column all the elements of its column arguments (see (A.27)). To describe the connections made by the cables (also called strings) we define C con =[c i,j ] c i;j ¼ 1 if the string terminates on the node n j 1 if the string originates from the node n j 0 if the string does not concern node n j 8 < : string connectivity matrix, n s 2n b where n s is the number of the strings and, to handle the connections in the vector form, we also define C con ¼½c 0 i;j  c 0 i;j ¼ I 3 if the string terminates on the node n j I 3 if the string originates from the node n j 0 33 if the string does not concern node n j 8 < : general connectivity matrix, 3n s 6n b The cables are represented by a string vectors: S¼ðs 1 s n s Þ¼NC T con string matrix,3n s r ¼vecðSÞ¼ s 1 . . . s n s 0 B @1 C A¼C con m 3n s 1 Fig. 3. Reference frames. M. Cefalo, J.M. Mirats-Tur / International Journal of Solids and Structures 48 (2011) 785–802 789 Author's personal copy To model the effects of undesired interactions with the environment, we consider the presence of some disturbances applied at the extreme of the bars: D=[d i,j ] d i;j ¼1 if the disturbance is applied to the node n j 0 otherwise disturbance connectivity matrix, n w 2n b where n w is the number of external forces applied to the nodes of the structure. Each row only contains one nonzero element. D¼d 0 i;j hi d 0 i;j ¼I 3 if disturbance is applied to the node n j 0 33 otherwise general disturbance connectivity matrix, 3n w 6n b W¼ðw 1 w n w Þdisturbances matrix, 3n w ; the set of all external forces acting on the system f¼vecðWÞ¼ w 1 . . . w n w 0 B @1 C A 3n w 1 Fig. 4 shows some of the quantities above introduced. We also define the tension vectors, which represent the forces applied by the cables to the bars. We have a tension vector for each couple of cable-bar connected. Ultimately, we have 2n s tensions applied to the bars. From the concept of ‘force coefficient’ we define as follows the tension applied by the cable ito the node j: t i;j ¼Cði;jÞ c i ks i ks i ks i kð20Þ The force coefficients c i , for an element in tension, must be always positive, hence, in Eq. (20) it is necessary to take into account the coefficient C(i,j) to coherently modify the sign of the vector. This comes from the fact that, for each cable, the strings have a fixed direction (even if arbitrary): the same for both nodes connected by the tendon. The coefficient C(i,j) may be +1 or 1, depending on the fact that the string associated to the cable, respectively, terminates or originates from the node jand this is why we need to correct the sign of a tension vector with respect to the string vector which it refers to. Let us define C¼diagf c 1  c n s gforce density matrix, n s n s ; where c i are the force density coefficients: forces for length unit C 0 ¼diagf c 1 ; c 1 ; c 1 ... c n s ; c n s ; c n s ggeneral force density matrix, 3n s 3n s F¼ðF 1 F 2 Þ¼ðf 1 f 2n b Þforces matrix,32n b ; the set of all resultant forces applied to the nodes: include tensions and external forces /¼vecðFÞ¼ / 1 / 2  ¼ f 1 . . . f 2n b 0 B @1 C A 6n b 1 It is easy to prove that: F¼NC T con CC con þWDð21Þ /¼C T con C 0 C con m þD T fð22Þ Observe the definition of the external forces and of the corresponding matrix D: the external forces may be directed anyway, not necessary along the cables. Adopting the elastic model for the cables, from the equations of a spring, we can also write: t i;j ¼k i Cði;jÞðks i kks i0 kÞ s i ks i kð23Þ where k i is the stiffness of the cable, ks i kthe actual length and ks i0 k the rest length. Fig. 4. Definitions of bar vectors, node vectors and forces. 790 M. Cefalo, J.M. Mirats-Tur / International Journal of Solids and Structures 48 (2011) 785–802 Author's personal copy Comparing Eq. (20) with Eq. (23) we obtain: k i ðks i kks i0 kÞ ¼ c i ks i kð24Þ In (21) the matrix C con is constant and well known (it depends on the mechanical connections). W and Dare well known as well, even if they may be time varying. Hence, in (21), the unknown variables are the positions of the bars nodes and C , which is function of the cables length and, therefore, function of N. From (24) it results c i ¼k i 1ks i0 k ks i k  ð25Þ Finally we define K= diag{k i } diagonal matrix of all cables stiffness coefficients, n s n s S 0 = diag{ks i0 k} matrix of the cables rest lengths, n s n s In the following we will treat the forces exerted by the cables like non conservative forces, even if we defined them on the base of the elastic model (they come from a potential function). The reason for this choice is twofold: first, in this way the final equations remain valid even if the model adopted to describe the cables changes, and, second, we will be able to consider the cables (or some of them) as actuated without the need to change anything. In this case the tensions on the cables do not necessary obey more to a model, they are directly imposed, for example, by means of engines and therefore they really are non-conservative. 3.2. Kinetic energy The kinetic energy of a tensegrity is the sum of the kinetic energy of the single bars. Accordingly, we will write first the kinetic energy for a single bar, which is, itself, sum of two components: the energy due to the translational motion and the component related to a rotational motion. To write the kinetic energy as a function of the center of mass and the quaternion representing the orientation of the bar, we need to introduce first the complementary kinetic energy. Let c i be the center of mass (expressed in the base frame) and x 1 i be the angular velocity vector (expressed in the body-fixed frame) of the ith bar. The complementary kinetic energy is defined as (A.26) T  i ¼1 2m i _ c T i _ c i þ1 2 x 1 i T J 1 3;i x 1 i ð26Þ ¼1 2m i _ c T i _ c i þ1 2 x 1 i T J 1 4;i x 1 i ð27Þ where m i is the mass of the bar and J 1 3;i is the central tensor of inertia (i.e. the tensor of inertia expressed in a frame positioned in the center of mass). The subscript ‘3’ is used to distinguish this matrix from the version valid with quaternions (J 1 4;i ), which is just a 4 by 4 matrix. Eq. (27) allows to express the complementary kinetic energy making use of quaternions in place of 3-dimensional vectors. The matrix J 1 4;i is defined as (Chou, 1992): J 1 4;i ¼ 1 2 traceðJ 1 3;i Þ0 0J 1 3;i ! where the operator trace(A) is the sum of the elements on the main diagonal of the matrix A. At last, in the hypothesis that the bar is a solid, regular cylinder (not hollow), we easily compute the central tensor of inertia, considering that all the elements outside the main diagonal (products of inertia) are null and that the elements on the main diagonal simply are the moments of inertia of the bar with respect to the main axis (of the bar-fixed frame): J 1 3;i ¼ 1 2 m i r 2 i 00 0 1 12 m i l 2 i 0 00 1 12 m i l 2 i 0 B B @1 C C Að28Þ where r i is the radius of the bar section and l i is the length of the bar. Eq. (27) is commonly called kinetic energy but strictly speaking it corresponds to the definition of the complementary kinetic energy (see Tabarrok and Rimrott, 1994; Rimrott et al., 1993). For mechanical scleronomic systems 4 it happens that the two functions have the same values when computed on the same state, even if this is not true in general. The complementary function may be useful to compute some partial differentiations with respect to variables that do not explicitly appear in the kinetic energy function. For our purposes, it will be useful to write the Lagrange’s equations in function of the quaternions e. Substituting (A.26) in (26), (A.22) in (27) we obtain T  i ¼T  i ð_ c i ;e i ;_ e i Þ: T  i ¼1 2m i _ c T i _ c i þ2_ e T i E T i J 1 3;i E i _ e i ¼1 2m i _ c T i _ c i þ2e T i _ E T i J 1 3;i _ E i e i ¼1 2m i _ c T i _ c i þ2_ e T i E þi J 1 4;i E þT i _ e i ¼1 2m 1 _ c T i _ c i þ2e T i _ E þi J 1 4;i _ E þT i e i The kinetic energy can be obtained from the complementary function applying the Legendre transformation. Let define the generalized momenta p i ¼@T  i @_ c i  T ¼m i _ c i ð29Þ g i ¼@T  i @_ e i  T ¼4E þi J 1 4;i E þT i _ e¼4E T i J 1 3;i E i _ e i ð30Þ Let q i be the vector of the generalized coordinates: q i ¼c i e i  ð31Þ It results T  i ¼T  i ðq i ;_ q i Þ. Consider now the complementary kinetic energy only as function of _ q:T  ¼T  ð_ qÞ. The definition of the kinetic energy (coming from the application of the Legendre transformation) is (see Arnold, 1978; Shivarama and Fahrenthold, 2004): T¼x T _ qT  ð_ qÞ where x¼ @T  @_ q  T . Finally, we have T i ¼T i ðp i ;e i ;g i Þ¼p T i _ c i þg T i _ e i T  i ð32Þ ¼1 2m i p T i p i þ1 8g T i E T i J 1 3;11 E i g i ¼1 2m i p T i p i þ1 8g T i E þi J 1 4;i1 E þT i g i ð33Þ This last expression, is particularly useful, since it allows to explicit the dependency of the kinetic energy from the quaternion e. At last, if n b is the number of the bars of a tensegrity structure, then the total kinetic energy is T¼X n b i¼1 T i ð34Þ 4 Are said scleronomics the systems for which the transformation equations between cartesian coordinates and generalized coordinates do not explicitly depends on time. A scleronomic constraint is a time-independent constraint, opposite to a rheonomic constraint, which is a time varying constraint. M. Cefalo, J.M. Mirats-Tur / International Journal of Solids and Structures 48 (2011) 785–802 791 Author's personal copy Let us compute now the total differential of both T i and T  i : dT  i ¼@T  i @_ c i  d_ c i þ@T  i @e i  de i þ@T  i @_ e i  d_ e i ¼p T i d_ c i þ@T  i @e i  de i þg T i d_ e i ð35Þ differentiating (32) and exploiting Eq. (35) we have dT i ¼p T i d_ c i þ_ c T i dp i þg T i d_ e i þ_ e T i dg i dT  i ¼_ c T i dp i þ_ e T i dg i @T  i @e i  de i ð36Þ while differentiating (33) we find dT i ¼@T i @p i  dp i þ@T i @e i  de i þ@T i @g i ! dg i ð37Þ Finally, comparing (36) with (37) we obtain a very important property, that we will exploit later onto determine the Lagrange’s equations: @T i @e i ¼@T  i @e i ð38Þ 3.3. Potential energy and non-conservative forces From now on we will consider the actions exerted by the cables like external non-conservative forces applied to the extremes of the bars (see Subsection 3.1). The potential energy is hence due only to the gravitational field. If the zaxis of the base frame is oriented in the vertical direction (upward direction), the potential energy function of a single bar is: V i ¼V i ðc i Þ¼m i g T 0 c i ð39Þ where g T 0 ¼00gðÞand the potential function of the overall system is: V¼X n b i¼1 V i ð40Þ To be able to apply the Euler–Lagrange’s equations we also need to define the non-conservative forces and moments acting on the system. Referring to Fig. 5, in hypothesis that the only external force is the friction, here modeled as viscous friction, the total force applied to the center of mass of a single bar is the vectorial sum of all forces applied at the extremes plus the gravity and the friction: f i ¼f 1;i þf 2;i @V i @c i  T b t;i _ c i ð41Þ where f 1 and f 2 are the sum of the forces applied, respectively, to the first and to the second endpoint of the bar and b t,i is the coefficient of the viscous friction associated to the translational motion of the ith bar. In term of f 1 and f 2 , the resulting total moment responsible for rotations is: l i ¼1 2b i f 1;i þ1 2b i f 2;i b r;i x i ¼1 2b i ðf 2;i f 1;i Þb r;i x i ð42Þ where b r,i is the coefficient of the viscous friction associated to the rotational motion and x i is the angular velocity of the bar computed in the base frame. Eqs. (41) and (42) are expressed in the base frame and refers to the cartesian coordinates. Starting from (42) we can determine now the generalized torques acting on quaternions. Let q i be the quasi-coordinates 5,6 associated with the components of the angular velocity vector: _ q 0 i ¼ x 0 i The virtual work made by the non-conservative forces is dW i ¼f 1;i þf 2;i b t;i _ c i  T dc i þ l T i dq 0 i ð43Þ where dq 0 i ¼dq 0 i ¼ x 0 i dt ¼2E þi de i ð44Þ From (43) and (44) we have dW i ¼f 1;i þf 2;i b t;i _ c i  T dc i þ2E þT i l i  T de i ð45Þ hence the generalized torque acting on e i is 2E þT i l i . Similar calculus can be done using quaternions in place of moment vectors for the rotational part of the kinetic energy. The variables can be referred to the base frame or to the body fixed frame: dW rot;i ¼ l 0 i T dq 0 i ¼ l 1 i T dq 1 i ð46Þ where l i and q i are the quaternions representing, respectively, the vectors l i and q i . Eq. (46) can be proven as follows. Considering that dq 0 i ¼ x 0 i dt ¼2E i ðe i Þ T de i ð47Þ dq 1 i ¼ x 1 i dt ¼2E þi ðe i Þ T de i ð48Þ we have dW rot;i ¼2E i ðe i Þ l 0 i  T de i ¼2E þi ðe i Þ l 1 i  T de i ð49Þ Fig. 5. Forces equivalent frames. 5 Quasi-coordinates were introduced to derive equations of motion using the Boltzmann–Hamel equations. They are velocity coordinates which are not simply the time derivative of position coordinates (see, for example, Nejmark and Fufaev, 1972 or Papastavridis, 2002 for more details). 6 Note the differences between qand q, which is the symbol used in (31) to refer to the set of generalized coordinates of a bar. 792 M. Cefalo, J.M. Mirats-Tur / International Journal of Solids and Structures 48 (2011) 785–802 Author's personal copy Recalling now that l 1 i ¼A 1 4;i l 0 i ¼E þi ðe i Þ T E i ðe i Þl 0 i , it can be easily shown that the last two expressions of Eq. (46) are the same. The expression referring to the base frame is, typically, more useful, because it requires to express the torque acting on the bar with respect to the base frame, in which it is naturally deduced from the forces expression. At last, observe that in Eq. (49),E i ðe i Þis the square version of E þT i in Eq. (45): the latter can be extracted from the former suppressing the first column (recall property (10)). 3.4. The Lagrange equations In this section we will obtain the equations of motion for a class-1 tensegrity system starting from the equations of motion of a single bar. Let L i ¼T i V i be the Lagrangian function of a single bar, the Euler–Lagrange equations are d dt @L @_ c i  T @L @c i  T ¼n c;i d dt @L @_ e i  T @L @e i  T ¼n e;i ð50Þ where n c,i and n e,i are the vectors of the generalized forces acting, respectively, on the state variables c i and e i . To apply Eq. (50),we start computing the following terms: @L i @c i  T ¼ @V i @c i  T ¼m i g 0 ð51Þ @L i @_ c i  T ¼@T i @p i @p i @_ c i  T ¼m i _ c i ð52Þ @L i @e i  T ¼@T i @e i  T ¼ @T  i @e  T ¼4 _ E þi J 1 4;i _ E þT i e i ð53Þ @L i @_ e i  T ¼@T @g i @g i @_ e i ! T ¼4E þi J 1 4;i E þT i _ e i ð54Þ where to obtain Eq. (53) the property (38) has been exploited. Finally, the equations of motion of a single bar are: € c i ¼1 m i ðf 1;i þf 2;i þm i g 0 b t;i _ c i þf d;i Þð55Þ d dt 4E þi J 1 4;i E þT i _ e i  þ4 _ E þi J 1 4;i _ E þT i e i ¼n e;i ð56Þ where f d,i is a generic disturbance force acting on the motion of the center of mass. To determine n e,i , we can consider the virtual work made by the torques applied to the bar. Recalling Eqs. (42) and (49) we write: n e;i ¼2E i 0 1 2 ðb i ðf 2;i f 1:i ÞÞb r;i x 0 i þ s 0 d;i "# ð57Þ where s 0 d;i represents a generic disturbance torque vector expressed with respect to the base frame. Computing the derivative in (56) we obtain: 4 _ E þi J 1 4;i E þT i _ e i þ4E þi J 1 4;i _ E þT i _ e i þ4E þi J 1 4;i E þT i € e i þ4 _ E þi J 1 4;i _ E þT i e i ¼n e;i ð58Þ which, exploiting (A.13) becomes 4E þi J 1 4;i E þT i € e i þ4E þi J 1 4;i _ E þT i _ e i ¼n e;i )€ e i ¼1 4E þi J 1 4;i1 E þT i n e;i E þi _ E þT i _ e i ð59Þ and, finally, applying (A.5) we have € e i ¼1 4E þi J 1 4;i1 E þT i n e;i E þi 1 0 0 0 0 B B B @1 C C C A _ e T i _ e i )€ e i ¼1 4E þi J 1 4;i1 E þT i n e;i e i _ e T i _ e i ð60Þ The system made of (55) and (56) becomes: € c i ¼1 m i f 1;i þf 2;i þm i g 0 b t;i _ c i þf d;i  ð61Þ € e i ¼1 4E þi J 1 4;i1 E þT i E i 0 b i ðf 2;i f 1;i Þb r;i x i þ s d;i  e i _ e T i _ e i ¼1 4E þi J 1 4;i1 E þT i E i 0 b i ðf 2;i f 1;i Þþ s d;i  þb r;i 2e_ e   e i _ e T i _ e i ð62Þ or, € c i ¼1 m i ðf 1;i þf 2;i þm i g 0 b t;i _ c i þf d;i Þð63Þ € e i ¼1 4E þi J 1 4;i1 0 l 0 0 0 B @1 C AA 1 3;i ðf 2;i f 1;i Þþ s 1 d;i 2b r;i E_ e i 2 6 6 6 43 7 7 7 5e i _ e T i _ e i ð64Þ where s 1 d;i is the disturbance torque expressed in the body fixed frame. Here the contribution of the friction, referring to the case of viscous friction, has been described in detail, to provide an example of how such kind of forces act on the generalized coordinate chosen. In the following, to simplify the notation, we will consider the friction included in the generic disturbance forces w i introduced in Section 3.1. We are now ready to write the equations of motion of a class-1 tensegrity system. We will shape two models: a matrix form model and a vector form model. The former is easier to obtain. It is also more compact and elegant, hence more convenient to write. The latter allows to quickly write the equations of motion in term of the state variables in the classical form _ x¼fðx;uÞ. This last form is more suitable for studying properties of a tensegrity by means of the classical analysis tools for non-linear systems based on the vectors fields. We start with the model in matrix form. Define the following set of generalized coordinates for a class-1 tensegrity structure: C¼ðc 1 c n b Þð65Þ E¼ð e 1 e n b Þð66Þ Observe that the total number of variables used to describe the system is 7n b , while the number of variables strictly necessary is 5n b , considering that the orientation of each bar may be represented only with two coordinates (the rotations around their symmetry axis will be neglected). Therefore (65) and (66) constitute a non-minimum set of coordinates. In what follows we will make extensive use of the notation introduced in Section 3.1. To completely explicit the equations of motion in function of the centers of mass and of the unit quaternions, we need to explicit Nand C in Eq. (21). The nodes coordinates can be written as follows: M. Cefalo, J.M. Mirats-Tur / International Journal of Solids and Structures 48 (2011) 785–802 793 Author's personal copy _ Q þT _ q¼ 1 0 0 0 0 B B B @1 C C C A _ q T _ qðA:5Þ where Qis defined like Ein Eq. (8) for unit quaternions. Differentiating once and twice Eq. (A.3) we also obtain, respectively Q_ qþ_ Qq ¼0 31 ðA:6Þ Q€ qþ€ Qq ¼0 31 ðA:7Þ Limiting now the attention to unit quaternions (shown below with the symbol eto differentiate them from the generic quaternions),fromthedefinitionand thepropertyofunitarynorm, wehave E þT E þ _ e¼ðI 4 ee T Þ_ e¼_ eeð_ e T eÞ T ¼_ e E T E  ¼E  E T ¼I 4 In the last equation, the signs + and over the matrices, must be used coherently: the property is valid only using the same sign over the two matrices involved in the product. Differentiating the unitary norm constraint with respect to time we can derive some other useful properties e T e¼1ðA:8Þ _ e T e¼0ðA:9Þ € e T eþ_ e T _ e¼0ðA:10Þ The unitary norm constraint can also be exploited to write E þT e¼ 1 0 0 0 0 B B B @1 C C C AðA:11Þ EE T ¼I 3 ðA:12Þ At last, differentiating (A.11) and (A.12), we obtain respectively: E þT _ eþ _ E þT e¼0 41 ðA:13Þ E_ E T þ_ EE T ¼0 31 ðA:14Þ Observe that E T E–I 4 and that E T Eis not invertible since qðE T EÞ6qðEÞ63. Besides, called X¼E T E, it is easy to prove that EX ¼E. Let us now summarize some properties concerning finite rotations of quaternions. Let r 0 and r 1 be, respectively, a vector quaternion before and after the rotation induced by the unit quaternion e. They can be interpreted as the same quaternion represented in the base frame and in the body fixed frame (whose orientation is precisely identified by e). The following equations are valid: r 0 ¼A 4 r 1 ðA:15Þ A 4 ¼ D E þ ðeÞEðe  Þ¼E þ ðeÞEðeÞ T ¼EðeÞ T E þ ðeÞ¼E þ ðe  Þ T Eðe  Þ ¼10 T 3 0 3 E þ E T ! ðA:16Þ r 1 ¼A 1 4 r 0 ðA:17Þ A 1 4 ¼E þ ðe  ÞEðeÞ¼E þ ðe  ÞEðe  Þ T ¼E þ ðeÞ T EðeÞ ¼EðeÞE þ ðeÞ T 10 T 3 0 3 EE þT ! ðA:18Þ From Eqs. (A.16) and (A.18) we also define the 3D rotation matrix and its inverse: A 3 ,E þ E T ðA:19Þ A 1 3 ¼EE þT ðA:20Þ Observe that A 1 4 ¼A T 4 as well as A 1 3 ¼A T 3 . Finally let x 0 be the quaternion associated to the angular velocity x expressed in base frame and x 1 be the same quaternion expressed in the body fixed frame. We have x 0 ¼2_ ee  ¼2e_ e  ðA:21Þ x 1 ¼2_ e  e¼2e  _ eðA:22Þ x 0 ¼A 4 x 1 ðA:23Þ x 1 ¼A T 4 x 0 ðA:24Þ besides, let x 0 and x 1 be the angular velocity vectors expressed, respectively, in the base frame and in the body fixed frame. It can be proven that x 0 ¼2E þ _ eðA:25Þ x 1 ¼2E_ eðA:26Þ Observe also that x =(0 x T ) T . For the proofs of the above mentioned properties and more details on quaternions (see Chou, 1992; Shivarama and Fahrenthold, 2004). Appendix B. Some properties on the matrix calculus It follows a brief list of definitions and properties of matrix operators recalled in this paper. 1. The operator vecð v 1  v n Þ where v 1 v n are vectors, defines a new vector made of all the input vectors: vecð v 1  v n Þ¼ v 1 . . . v n 0 B B @1 C C AðA:27Þ The same operator applied to a matrix is meant as applied to the column vectors of such matrix. 2. The Kronecker product between two matrices is defined as follows: KroneckerðA;BÞ¼AB¼ a 1;1 B... a 1;n B . . ... . a n;1 Ba n;m B 0 B B @1 C C AðA:28Þ If [A]=mnand [B]=rs, then the Kronecker product between A and Bis a matrix mrns. Some important properties are: (a) A(B+C)=AB+AC (b) (A+B)C=AC+BC (c) If A,B,C, and Dare four matrices which allow to compute the products ACand BD, we have ðABÞðCDÞ¼ðACÞðBDÞ (d) If Aand Care two matrices which allow to compute the products ACwe have: 800 M. Cefalo, J.M. Mirats-Tur / International Journal of Solids and Structures 48 (2011) 785–802 Author's personal copy ðABÞC¼ðACÞB and if it possible to compute BCwe can also write: ðABÞC¼AðBCÞ (e) If Aand Bare two invertible matrices, we have ðABÞ 1 ¼A 1 B 1 3. Given two 3-dimensional vectors u 0 and v 0 represented in the same frame, let u 1 be the representation of u 0 in a different frame and Rbe the rotation matrix which define the orientation of the frame 1 with respect to the frame 0: u 0 =Ru 1 . Then the following holds: u 0  v 0 ¼e U 0 v 0 ¼Ru 1  v 0 ¼Rðu 1 R T v 0 Þ¼Re U 1 R T v 0 from which we have e U 0 ¼Re U 1 R T ðA:29Þ 4. Let Q¼ 10 0 0 0 zfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflffl{ ncomponents 0 0 01 0  0 |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} ncomponents  0 . . ... . 00 01 0 B B B B B B B B B @ 1 C C C C C C C C C A ðA:30Þ dimensionally it is [Q]=nn 2 . If e i are m1 vectors, it results: diagðe T 1 e 1 e T n e n Þ¼ e T 1 e 1 ... 0 . . ... . 0e T n e n 0 B B @1 C C A ¼QðI n EÞ T ðI n EÞQ T ¼QðI n E T ÞðI n EÞQ T ¼Q½I n ðE T EÞQ T where E=(e 1 e n ) 5. Let E¼ðe 1 e n Þ and e v ¼vecðEÞ where e i are m1 vectors. We have: E¼I m 0 m ...0 m |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} nblocks ... 0 m ...I m ...0 m ... 0 m ...0 m I m |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ntimes 0 B @1 C A ðI n e v Þ ¼10...0... 0...1...0... 0...01ðÞI m ½ðI n e v Þ where 0 m =0 mm . 6. The matrix diagðe T 1 e 1 e T n e n Þcan also be computed starting from e v instead of E. Consider that e i ¼ e i;1 . . . e i;m 0 B B @1 C C A¼0 m ... I m ... 0 m ðÞe v where I m is the ith block mm. Transposing we have e T i ¼e T v 0 m . . . I m . . . 0 m 0 B B B B B B @ 1 C C C C C C A Finally, we can write diagðe T 1 e 1 ...e T n e n Þ¼ðI n e T v ÞQ 0 Q 0T ðI n e v Þ where Q 0 ¼ I m 0 mm . . . 0 mm 9 > > > > = > > > > ; nblocks mm 0 mm I m . . . 0 mm . . ... . 0 mm . . . 0 mm I m 0 B B B B B B B B B B B B B B B B B B B B B B B B B B B B B @ 1 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C A ðA:31Þ dimensionally it is [Q 0 ]=n 2 mnm. Observe that Q 0 =Q T I m . References Aldrich, J., 2004. Control synthesis for a class of light and agile robotic tensegrity structures. PhD thesis, University of California. Arnold, V.I., 1978. Mathematical Methods of Classical Mechanics. Springer-Verlag, USA. Calladine, C., Pellegrino, S., 1992. Further remarks on first order infinitesimal mechanisms. International Journal of Solids and Structures 29, 2119–2122. Chou, J.C.K., 1992. Quaternion kinematic and dynamic differential equations. IEEE Transaction on Robotics and Automation 8 (1), 53–64. Connelly, R., 1999. Rigidity Theory and Applications. In: Chapter. Tensegrity Structures: Why are They Stable? Kluwer Academic/Plenum, pp. 47–54. Connelly, R., Back, A., 1998. Mathematics and tensegrity. American Scientist 86, 142151. Fuller, R., 1962. Tensile-integrity structures, United States Patent 3063521. Hamilton, W., 1847. On quaternions. Proceedings of the Royal Irish Academy 3, 1–16. Hamilton, W., 1866. Elements of Quaternions. Longmans, Green, & co.. Hanaor, A., 1988. Prestressed pin-jointed structures-flexibility analysis and prestress design. Computers and Structures 28 (6), 757–769. Hernandez, S., Mirats-Tur, J., 2008. Tensegrity frameworks: static analysis review. Mechanism and Machine Theory vol. 43 (7), 859–881. Kanchanasaratool, N., Williamson, D., 2002. Modelling and control of class NSP tensegrity structures. International Journal of Control 75, 123–139. Masic, M., Skelton, R., 2004. Open-loop control of class-2 tensegrity towers. Proceedings of SPIE 5383, 298–308. Mirats-Tur, J., Hernández, S., 2009. Tensegrity frameworks: dynamic analysis review and open problems. International Journal of Mechanism and Machine Theory 44 (1), 1–18. Mirats-Tur, J., 2010. On the movement of tensegrity structures. Internationa; Journal of Space Structures 25 (1), 1–14. Motro, R., Najari, S., Jouanna, P., 1986. Static and dynamic analysis of tensegrity systems. In: Proceedings of the ASCE International Symposium on Shell and Spatial Structures: Computational Aspects. Springer, pp. 270–279. Nejmark, J., Fufaev, N., 1972. Dynamics of Nonholonomic Systems. Translations of Mathematical Monographs, AMS, 33. Papastavridis, J., 2002. Analytical Mechanics, A Ccomprehensive Treatise on the Dynamics of Constrained Systems. Oxford University Press. Paul, C., Valero-Cuevas, F., Lipson, H., 2006. Design and control of tensegrity robots for locomotion. IEEE Transactions on Robotics 22 (5), 944–957. Pugh, A., 1976. An Introduction to Tensegrity. University of California Press, USA. Rimrott, F.P.J., Szczygielski, W.M., Tabarrok, B., 1993. Kinetic energy and complementary kinetic energy in gyrodynamics. Journal of Applied Mechanics, Transaction of the ASME 60, 398–405. Roth, B., Whiteley, W., 1981. Tensegrity frameworks. Transaction of the American Mathematical Society 265, 419–446. Shivarama, R., Fahrenthold, E.P., 2004. Hamilton’s equations with Euler parameters for rigid body dynamics modeling. Journal of Dynamic Systems, Measurement, and Control, Transaction of the ASME 126, 124–130. Siciliano, B., Sciavicco, L., Villani, L., Oriolo, G., 2009. Robotics: Modelling, Planning and Control. Springer-Verlag, USA. Skelton, R., Oliveira, M.D., 2009. Tensegrity Systems. Springer. M. Cefalo, J.M. Mirats-Tur / International Journal of Solids and Structures 48 (2011) 785–802 801 Author's personal copy Skelton, R., Pinaud, J., Mingori, D., 2001. Dynamics of the shell class of tensegrity structures. Journal of the Franklin Institute 338, 255–320. Snelson, K., 1965. Discontinuous compression structures. United States Patent 3169611. Sultan, C., 1999. Modeling, design and control of tensegrity structures with applications. PhD thesis, Purdue Univeristy. Sultan, C., Corless, M., Skelton, R., 2002. Linear dynamics of tensegrity structures. Engineering Structures vol. 24. Tabarrok, B., Rimrott, F.P.J., 1994. Variational Methods and Complementary Formulations in Dynamics. Springer. Tarnai, T., 1989. Higher-order infinitesimal mechanisms. Acta Technical Academy Science 102, 363–378. Vassart, N., Laporte, R., Motro, R., 2000. Determination of mechanism’s order for kinematically and statically indeterminate systems. International Journal of Solids and Structures 37, 3807–3839. 802 M. Cefalo, J.M. Mirats-Tur / International Journal of Solids and Structures 48 (2011) 785–802