The Equilibrium and Form-Finding of General

Tensegrity Systems with Rigid Bodies

Shuo Ma, Muhao Texas A&m University Chen, Zhangli Peng, Xingfei Yuan,
Robert E Skelton

To cite this version:

Shuo Ma, Muhao Texas A&m University Chen, Zhangli Peng, Xingfei Yuan, Robert E Skelton. The
Equilibrium and Form-Finding of General Tensegrity Systems with Rigid Bodies. Engineering Struc-
tures, In press. �hal-03569335�

HAL Id: hal-03569335

[Link]
Submitted on 12 Feb 2022

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.

Distributed under a Creative Commons Attribution - NonCommercial - NoDerivatives| 4.0

International License
 The Equilibrium and Form-Finding of General Tensegrity
Systems with Rigid Bodies
Shuo Maa,b, Muhao Chenc,*, Zhangli Pengd, Xingfei Yuane,*, Robert E. Skeltonc
a College of Civil Engineering, Zhejiang University of Technology, Hangzhou, Zhejiang 310023, China

b Key Laboratory of Space Structures of Zhejiang Province, Hangzhou, Zhejiang 310058, China

c Department of Aerospace Engineering, Texas A & M University, College Station, Texas 77840, United States

d Department of Biomedical Engineering, University of Illinois at Chicago, Chicago, Illinois 60607, United
States

e College of Civil Engineering and Architecture, Zhejiang University, Hangzhou, Zhejiang 310058, China

Abstract
We develop a general approach to study the equilibrium and form-finding of any general tensegrity systems with
rigid bodies. The equilibrium equations are derived in an explicit form in terms of a nodal coordinate and
orientation parameter as the minimal coordinate. The nodal vector consists of nodes (either free or pinned) in
the pure bar-string tensegrity network and nodes on the rigid bodies (those connected to the pure bar-string
tensegrity network). Based on the Lagrangian method, the nonlinear statics of the general tensegrity system in
terms of the minimal coordinate is first given. Then, we linearize the statics equation and obtain its equivalent
form, in terms of the force vector of the compressive and tensile members, for the analysis of structure
equilibrium configurations and prestress modes. To study the system's stability and have a comprehensive insight
into the materials and structure members, we present the tangent stiffness matrix as a combination of prestress,
material, and geometric information of the structure. It is also shown that without rigid bodies, the governing
equations of the general tensegrity system yields to the classical tensegrity structure (pure string-bar network).
Form-finding of general tensegrity is implemented based on solving the nonlinear equilibrium equation, where
the modification of tangent stiffness matrix and line search algorithm is used. Numerical examples are given to
demonstrate the capability of our developed method in finding the feasible prestress modes, conducting form-
finding and prestress designs, and checking the structural robustness of any tensegrity systems with rigid bodies.

KEYWORDS: Generalized tensegrity, Rigid body, Tensegrity equilibrium, Minimal coordinate, Form-finding

*Corresponding authors.

E-mail addresses: yuanxf@[Link] (X. Yuan), muhaochen@[Link] (M. Chen)

1 Introduction
Tensegrity is a conjunction of two words (tension and integrity) which was first proposed by Buckminster Fuller
[1] for the art form by Ioganson (1921) and Snelon (1948) [2]. In their work, they never assumed there were no
rigid bodies in the tensegrity structures. And in fact, the tensegrity sculpture built by Snelson in 1948 is two X-
shape rigid bodies stabilized by several cables. However, it is probably because bars and strings are more
efficient in taking compression, provide more accurate models (uncertainty is only along with the axially loaded
members), and it is complicated to model the irregular shape of the rigid bodies, most of the literature focus on
pure stable bar-strings networks.

Indeed, after decades of study, the pure bar-string tensegrity structures have shown their many advantages in
lightweight structure topology design [3-6], engineering structures [7,8,5], soft robotics [9,10], deployable
structures [11-13], energy absorption [14-16], etc. But for many engineering structures, we must include the
rigid bodies, i.e., the deck of the bridges, the roof of the shelters, the shell of cable domes, the D-section of the
airfoils, and the shield of space structures. To deal with these rigid bodies in their tensegrity structure design,
many researchers have proposed their compromised solutions to the rigid body tensegrities. For example,
Carpentieri et al. [17] separated the minimal mass design of the tensegrity bridge structure and its deck. Laccone
et al. [18] analyzed the cable-tensioned dome and its glass shell by the nonlinear finite element analysis software
Straus7. Levin et al. [19] studied the rigid body spine mechanics based on the tensegrity-truss model. Chen and
Jiang [20] used parallel mechanism theory to compute the stiffness of a fish, made of a set of rigid ribs stabilized
by strings. Chen et al. [21] decoupled the force analysis of a tensegrity space habitat and its shield. However,
none of these approaches started from the fundamental governing equations of the whole system and developed
a general approach to the analysis of tensegrity systems with rigid bodies. It is also worth mentioning that few
software packages have the compatibility of simulating tensegrity systems with rigid bodies. For example, Wang
et al. [22] modeled tensegrity swimmer and rigid bodies in the MuJoCo simulator and studied the data-based
control methods. Sun et al. [23] studied a tensegrity foot with a rigid board and universal joint in ADAMS.
Pajunen et al. [16] implemented ABAQUS to analyze the 3D-printable tensegrity lander with rigid joints.
However, these commercial packages are costly, require much experience, and the insight of the algorithm is
not clear.

In the past years, a few attempts have been made to study tensegrity with rigid body models analytically. For
example, for the static analysis, Hangai and Wu [24] proposed kinematics and equilibrium equations to study
the behaviors of a hybrid structure that consists of cables and rigid structures. Wang et al. [25] derived the statics
equilibrium equation of general tensegrity and used the mixed-integer linear programming method for the
topology design. Chen and Jiang [26] derived the total stiffness of a general tensegrity structure in an explicit
form and developed a set of sufficient and necessary conditions to guarantee the stability of the tensegrity
structures. For the dynamics analysis, Nagase and Skelton [27] used non-minimal coordinates to write the
2
dynamics equations of tensegrity by assuming the compression members are rigid bodies. Kan et al. [28,29]
studied the nonlinear dynamics of clustered tensegrity with rigid bodies by using the configuration of the attached
rigid bodies as the generalized coordinate. Li et al. [30] studied the kinodynamic planning of cable-driven
tensegrity manipulators composed of clustered cables and rigid bodies. However, the equilibrium theory in most
of the work is in a complicated form and limited to structures with small deformations. Moreover, there is an
increasing interest in using tensegrity structures to build robotics due to the many advantages of tensegrity
structure, i.e., mass saving, control energy efficiency, abundant equilibrium states, etc. In many of the tensegrity
robot applications, rigid bodies cannot be avoided. The current equilibrium theories and form-finding methods
of the tensegrity system with rigid bodies are still limited. It is critical to have an efficient form-finding approach
to find the configurations of the whole system to enlarge the applications of tensegrity systems. To this end, we
derived a general approach to the nonlinear equilibrium equations and proposed a corresponding form-finding
method to the tensegrity system with rigid bodies. In this study, the tensegrity with pure axial form elements is
referred to as the traditional tensegrity, while the tensegrity with rigid bodies is called the general tensegrity.

The paper is structured as follows. Section 2 presents the tensegrity and rigid body notations. Section 3 derives
the kinematics of the system. Section 4 gives the nonlinear and linearized statics equations. Section 5 shows the
form-finding approach to the tensegrity systems with rigid bodies. Section 6 summarizes the conclusions.

2 Notations of tensegrity systems with rigid bodies

Figure 1 Diagram of tensegrity with rigid bodies, bar (𝑏) and string (𝑠) vectors are marked in black and red.

3
2.1 Nodal coordinates of the system and its components
The tensegrity system with rigid bodies is composed of bars, strings, and rigid bodies, as shown in Figure 1. The
rigid bodies in the tensegrity structures are connected by the strings and bars nodes on the rigid bodies. We name
the nodes on the rigid body as rigid body nodes. The nodes only on the bars and strings are free tensegrity nodes,
and the other nodes in the fixed point are the pinned tensegrity nodes. The position of all the nodes can be
expressed in any frame, and we choose to label them in the Cartesian coordinates in an inertially fixed frame by
a nodal vector. Assume there are 𝑛𝑛 number of nodes, the 𝑋, 𝑌, and 𝑍-coordinates of the 𝑖th node 𝒏𝑖 ∈ ℝ3 in
the vector form is 𝒏𝑖 = [𝑥𝑖 𝑦𝑖 𝑧𝑖 ]𝑇 . By stacking 𝒏𝑖 for 𝑖 = 1, 2, ⋯ , 𝑛𝑛 together, we can get the nodal vector
𝒏 ∈ ℝ3𝑛𝑛 for the whole structure:
𝑇
𝒏 = [𝒏1𝑇 𝒏𝑇2 ⋯ 𝒏𝑇𝑛𝑛 ] , (1)

and its equivalent matrix form[31] 𝑵 ∈ ℝ3×𝑛𝑛 is:

𝑵 = [𝒏1 𝒏2 ⋯ 𝒏𝑛𝑛 ]. (2)

Note that one can simply obtain the nodal coordinate vector 𝒏 by vectorizing the nodal coordinate matrix 𝑵:

𝒏 = 𝑣𝑒𝑐(𝑵) = 𝑵(: ), (3)

where 𝑣𝑒𝑐(𝑵) is an operator that stacks all the columns of matrix 𝑵 into one vector. Normally, the positions of
some of the nodes in the structure are fixed/pinned in certain directions. Let there be 𝑛𝑎 degree of freedom of
free tensegrity nodes, 𝑛𝑏 degree of freedom of pinned tensegrity nodes, and 𝑚 rigid bodies with a total number
of 𝑛𝑞 degree of freedom of the rigid nodes. Suppose there are 𝑧𝑖 number of nodes in the 𝑖th rigid body. To deal
with the constraints, we distinguish the free tensegrity nodes, pinned tensegrity nodes, and the jth node in the ith
rigid body by introducing three kinds of vectors 𝒂 ∈ ℝ𝑛𝑎 , 𝒃 ∈ ℝ𝑛𝑏 , and 𝒒𝑖𝑗 ∈ ℝ𝑛𝑞 :

𝒂 = [𝑎1 𝑎2 ⋯ 𝑎𝑛𝑎 ]T, (4)

𝒃 = [𝑏1 𝑏2 ⋯ 𝑏𝑛𝑏 ]T, (5)

𝒒𝑖𝑗 = [𝑞𝑖𝑗𝑥 𝑞𝑖𝑗𝑦 𝑞𝑖𝑗𝑧 ]T , (𝑖 = 1, 2, ⋯ , 𝑚; 𝑗 = 1,2, ⋯ , 𝑧𝑖 ), (6)

where the values of 𝑎𝛼 (𝛼 = 1, 2, ⋯ , 𝑛𝑎 ) , 𝑏𝛽 (𝛽 = 1, 2, ⋯ , 𝑛𝑏 ) and 𝑞𝑖𝑗𝑥 , 𝑞𝑖𝑗𝑦 , 𝑞𝑖𝑗𝑧 (𝑖 = 1, 2, ⋯ , 𝑚; 𝑗 =

1,2, ⋯ , 𝑧𝑖 ) are the indices of the entries corresponding to the free tensegrity nodes, pinned tensegrity nodes, and
the 𝑋, 𝑌, 𝑍 freedom of the jth node in the ith rigid body in the nodal vector 𝒏. We use vectors 𝒏𝑎 , 𝒏𝑏 , and 𝒏𝑞𝑖𝑗

to label the nodal coordinate of the free node, pinned node, and the jth node in the ith rigid body. And 𝑬𝑛𝑎 ∈
ℝ3𝑛𝑛 ×𝑛𝑎 , 𝑬𝑛𝑏 ∈ ℝ3𝑛𝑛 ×𝑛𝑏 , and 𝑬𝑛𝑞𝑖𝑗 ∈ ℝ3𝑛𝑛 ×3 are the location matrices to extract vectors 𝒏𝑎 , 𝒏𝑏 , and 𝒏𝑞𝑖𝑗
from the vector 𝒏:

4
 𝑬𝑛𝑎 (: , 𝑘) = 𝑰3𝑛𝑛 (: , 𝑎𝑘 ), 𝑬𝑛𝑏 (: , 𝑘) = 𝑰3𝑛𝑛 (: , 𝑏𝑘 ), 𝑬𝑛𝑞𝑖𝑗 = 𝑰3𝑛𝑛 (: , [𝑞𝑖𝑗𝑥 𝑞𝑖𝑗𝑦 𝑞𝑖𝑗𝑧 ]), (7)

where 𝑰3𝑛𝑛 is the identity matrix in 3𝑛𝑛 order. Thus, we have the following:

𝒏𝒂 = 𝑬𝑻𝑛𝑎 𝒏, 𝒏𝒃 = 𝑬𝑇𝑛𝑏 𝒏, 𝒏𝑞𝑖𝑗 = 𝑬𝑇𝑛𝑞𝑖𝑗 𝒏. (8)

The nodal coordinate of the whole structure is obtained by summing all the free tensegrity nodes, pinned
tensegrity nodes, and rigid body nodes:
𝑚 𝑧𝑖

𝒏 = 𝑬𝑛𝑎 𝒏𝑎 + 𝑬𝑛𝑏 𝒏𝑏 + ∑ ∑ 𝑬𝑛𝑞𝑖𝑗 𝒏𝑞𝑖𝑗 . (9)

𝑗=1 𝑘=1

The 𝑖th (𝑖 = 1, 2, ⋯ , 𝑚) rigid body nodal coordinate vector is obtained by stacking the nodal coordinate of the
𝑧𝑖 rigid-body nodes:
𝒏𝑞𝑖1
𝒏𝑞𝑖2
𝒏𝑞 𝑖 = ⋮ . (10)
[𝒏𝑞𝑖𝑧 ] 𝑖

The location matrix corresponding to the 𝑖th (𝑖 = 1, 2, ⋯ , 𝑚) rigid body nodes is:

𝑬𝑛𝑞𝑖 = [𝑬𝑛𝑞𝑖1 𝑬𝑛𝑞𝑖2 ⋯ 𝑬𝑛𝑞𝑖𝑧 ].

𝑖 (11)

Then, the nodal coordinate vector of the 𝑖th (𝑖 = 1, 2, ⋯ , 𝑚) rigid body can be calculated by:

𝒏𝑞𝑖 = 𝑬𝑇𝑛𝑞𝑖 𝒏, (12)

2.1.1 Connectivity matrix

A connectivity matrix provides the connection pattern of all the nodes in the structure. Let 𝑪 ∈ ℝ𝑛𝑒×𝑛𝑛 be the
connectivity matrix of the tensegrity systems with rigid bodies, where 𝑛𝑒 is the total number of axially loaded
members (bars and strings). The 𝑖th (𝑖 = 1,2, ⋯ , 𝑛𝑒 ) row of 𝑪, denoted as 𝑪𝑖 = [𝑪](𝑖,:) ∈ ℝ1×𝑛𝑛 , represents

connectivity information of the 𝑖th element in the structure. Suppose the 𝑖th member is from the 𝑗th node to the
𝑘th node. The 𝑟th (𝑖 = 1,2, ⋯ , 𝑛𝑛 ) entry the 𝑖th row of 𝑪 satisfies:

−1 , 𝑟=𝑗
[𝑪]𝑖𝑟 = { 1 , 𝑟=𝑘 . (13)
0 , 𝑟 = 𝑒𝑙𝑠𝑒

2.1.2 The geometry of axial elements

An axial element vector denotes the start and end nodes of an axial element (bar or string). For example, the 𝑖th
axial element vector 𝒉𝑖 ∈ ℝ3×1 is:

5
 𝒉𝑖 = 𝒏𝑘 − 𝒏𝑗 = 𝑪𝑖 ⨂𝑰3 𝒏. (14)

where the symbol ⨂ represents the Kronecker product. Stacking all the axial elements into a structure element
matrix 𝑯 ∈ ℝ3×𝑛𝑒 , we get:

𝑯 = 𝑵𝑪𝑇 . (15)

The present length of the 𝑖th axial element is:

1
𝑙𝑖 = ‖𝒉𝑖 ‖ = (𝒏𝑇 (𝑪𝑇𝑖 𝑪𝑖 ) ⊗ 𝐈3 𝒏)2 . (16)

Rest length is the length of an axial element with no tension or compression. We use the subscript 0 to denote
the rest length of an axial element, i.e., the rest length of the 𝑖th axial element is 𝒍0𝑖 . The length and rest length
vector of all the axial elements are:

𝒍0 = [𝑙01 𝑙02 ⋯ 𝑙0𝑛𝑒 ]𝑇 , (17)

𝒍 = [𝑙1 𝑙2 ⋯ 𝑙𝑛𝑒 ]𝑇 . (18)

2.1.3 Stiffness of axial elements

Let the cross-sectional area, secant modulus, and tangent modulus of the 𝑖 th element be 𝐴𝑖 , 𝐸𝑖 , and 𝐸𝑡𝑖 ,
respectively. Then, the cross-sectional area, secant modulus, and tangent modulus vector of the structure
𝑨, 𝑬, 𝑬𝑡 ∈ ℝ𝑛𝑒 can be written as:

𝑨 = [𝐴1 𝐴2 ⋯ 𝐴 𝑛𝑒 ] 𝑇 , (19)

𝑬 = [𝐸1 𝐸2 ⋯ 𝐸𝑛𝑒 ]𝑇 , (20)

𝑬𝑡 = [𝐸𝑡1 𝐸𝑡2 ⋯ 𝐸𝑡𝑛𝑒 ]𝑇 . (21)

The internal force of the 𝑖th element is 𝑡𝑖 = 𝐴𝑖 𝜎𝑖 = 𝐸𝑖 𝐴𝑖 (𝑙𝑖 − 𝑙0𝑖 )/𝑙0𝑖 , the internal force vector of the structure
𝒕 ∈ ℝ𝑛𝑒 can be written as:

𝒕 = [𝑡1 𝑡2 ⋯ ̂𝑨
𝑡𝑛𝑒 ]𝑇 = 𝑬 ̂ 𝒍̂−1
0 (𝒍 − 𝒍0 ), (22)

̂ is an operator that converts vector 𝑬 into a diagonal matrix.

where 𝑬

2.2 Notations of the rigid bodies

2.2.1 Orientation matrix of rigid bodies

Unlike the bars and strings in the rigid body tensegrity, one can use the nodal vector to describe the exact attitude
of these axial elements. To describe the attitude of a rigid body, an orientation matrix must be included to show
the transition process. There are many approaches to achieve this goal, i.e., Euler angle, Euler principal axis, and

6
quaternion. We chose the Euler angle approach because it is a minimal coordinate method to describe the attitude
of rigid bodies. In this problem, we implemented a simple (1-2-3) orientation set, which means to rotate α, β, γ
about the principal axis of 𝒃1 , 𝒃2 , 𝒃3 in sequence in the body-fixed frame. The attitude parameter 𝝋 is the vector
composed of Euler angle:
𝛼
𝝋 = [𝛽 ]. (23)
𝛾

The attitude matrix is [32]:

𝑹(𝛼, 𝛽, 𝛾) = 𝑹3 (𝛾)𝑹2 (𝛽)𝑹1 (𝛼) =

cos 𝛾 cos 𝛽 cos 𝛾 sin 𝛽 sin 𝛼 + sin 𝛾 cos 𝛼 − cos 𝛾 sin 𝛽 cos 𝛼 + sin 𝛾 sin 𝛼 (24)
[− sin 𝛾 cos 𝛽 − sin 𝛾 sin 𝛽 sin 𝛼 + cos 𝛾 cos 𝛼 sin 𝛾 sin 𝛽 cos 𝛼 + cos 𝛾 sin 𝛼 ].
sin 𝛽 − cos 𝛽 sin 𝛼 cos 𝛽 cos 𝛼

Even though the Euler angle has kinematic singularities for the value of 𝛽 = 0, this is only a problem in
calculating the velocity of orientation parameters from angular velocities. For solving the static equilibrium and
form-finding of general tensegrities, there is no such problem using the Euler angle as the orientation parameter.

2.2.2 Mass center of rigid body

Let the mass center of the 𝑖th rigid body be 𝒏𝑐𝑖 ∈ ℝ3×1 . Normally, the position of the mass center can be given
by measuring the mass distribution of the rigid body in an experiment. However, in the static analysis, the
equilibrium of total force and moment is independent of the choice of the mass center. For simplicity, we can
directly use the geometry center of the 𝑖th rigid body nodes as the mass center:
1 (25)
𝒏𝑐 𝑖 = 𝑰 ⨂𝑰3𝒏𝑞𝑖 ,
𝑧𝑖 1,𝑧𝑖

where 𝑰1,𝑧𝑖 ∈ ℝ1×𝑧𝑖 is an all-ones vector with 𝑧𝑖 columns, and 𝑧𝑖 is the number of rigid body nodes in the ith
rigid body. Substitute Eq.(12) into Eq. (25), one can compute the mass center from the nodal coordinate vector
of the structure:

𝒏𝑐𝑖 = 𝑬𝑇𝑛𝑐 𝒏, (26)

𝑖

where 𝑬𝑛𝑐 is:

𝑖

1 (27)
𝑬𝑛𝑐 = 𝑧 𝑬𝑛𝑞 𝑰𝑧𝑖 ,1 ⨂𝑰3.
𝑖 𝑖 𝑖

2.2.3 Nodal coordinate of rigid bodies

If there is translation or rotation of the rigid bodies, the nodal coordinate of the 𝑗th node on the 𝑖th rigid body
𝒏𝑞𝑖𝑗 is:

7
 𝒏𝑞𝑖𝑗 = 𝒏𝑐𝑖 + 𝒓𝑖𝑗 , (28)

𝒓𝑖𝑗 = 𝑹𝑇𝑖 (𝒏𝑞𝑖𝑗0 − 𝒏𝑐𝑖0 ) = (𝑬𝑇𝑛𝑐𝑖 − 𝑬𝑇𝑛𝑞𝑖𝑗 ) 𝒏, (29)

where 𝒓𝑖𝑗 is the vector from the center of mass 𝒏𝑐𝑖 to the 𝑗th node in the 𝑖th rigid body, 𝒏𝑞𝑖𝑗0 and 𝒏𝑐𝑖0 is the
nodal coordinate vector of the 𝑗th node and the mass center of the 𝑖th rigid body in the body-fixed frame. 𝑹𝑖 is
the attitude matrix of the 𝑖th rigid body.

2.3 Minimal coordinate of the system

The minimal coordinate 𝑼 ∈ ℝ𝑛𝑈 is used to represent the position of the free tensegrity nodes and the rigid
bodies:
𝒏𝑎
𝑼1
U= 𝑼2 , (30)
⋮
[ 𝑼𝑚 ]

where 𝑼𝑖 is the minimal coordinate for the 𝑖th rigid body, including the position of the mass center 𝒏𝑐𝑖 ∈ ℝ3 and
attitude parameter 𝝋𝑖 ∈ ℝ3:

𝒏𝑐
𝑼𝑖 = [ 𝝋 𝑖 ], (31)
𝑖

The location matrix is used to locate minimal coordinate of free tensegrity nodes and rigid bodies:
𝒏𝑎
𝒏𝑐1
[𝝋 ]
1
𝒏𝑐2
𝑼 = [𝑬𝑈𝑎 [𝑬𝑈𝑐1 𝑬𝑈𝝋1 ] [𝑬𝑈𝑐2 𝑬𝑈𝝋2 ] ⋯ [𝑬𝑈𝑐𝑚 𝑬𝑈𝝋𝑚 ]] [ ] .
𝝋2 (32)
⋮
𝒏𝑐𝑚
[[ 𝝋𝑚 ]]

The nodal coordinate vector of free nodes, mass center, Euler angle, and minimal coordinate of the 𝑖th rigid body
is:
𝒏𝑎 = 𝑬𝑇𝑈𝑎 𝑼, 𝒏𝑐𝑖 = 𝑬𝑇𝑈𝑐𝑖 𝑼, 𝝋𝑖 = 𝑬𝑇𝑈𝝋 𝑼, 𝑼𝑖 = 𝑬𝑇𝑈𝑖 𝑼. (33)
𝑖

𝑬𝑈𝑐 and 𝑬𝑈𝜑 is used to extract the mass center and Euler angle information of all rigid bodies:

𝑬𝑈𝑐 = [𝑬𝑈𝑐1 𝑬𝑈𝑐2 ⋯ 𝑬𝑈𝑐𝑚 ], 𝑬𝑈𝝋 = [𝑬𝑈𝝋1 𝑬𝑈𝝋2 ⋯ 𝑬𝑈𝝋𝑚 ]. (34)

𝑬𝑈𝑖 is used to extract the minimal coordinate of the 𝑖th rigid body:

8
 𝑬𝑈𝑖 = [𝑬𝑈𝑐𝑖 𝑬𝑈𝝋 ].
𝑖
(35)

3 Kinematics of the rigid body

3.1 Attitude kinematics

The angular velocity vector of the 𝑖th rigid body in the inertial frame is [33]:
𝜔1
𝝎𝑖 = [𝜔2 ] = [𝑩𝑖 ]𝝋̇𝑖 , (36)
𝜔3

The 𝑩𝑖 matrix for the Euler angle (1-2-3) orientation set is:

1 0 sin 𝛽
𝑩𝑖 = [0 cos 𝛼 − cos 𝛽 sin 𝛼 ]. (37)
0 sin 𝛼 cos 𝛼 cos 𝛽

3.2 Transformation matrix

The velocity vector of the 𝑗th node on the 𝑖th rigid body is:

𝒏̇ 𝑞𝑖𝑗 = 𝒏̇ 𝑐𝑖 + 𝝎𝑖 × 𝒓𝑖𝑗 . (38)

Substitute Eq.(36) into Eq. (38), we will have:

d𝒏𝑞𝑖𝑗 d𝒏𝑐𝑖
= − 𝒓𝑖𝑗 × 𝝎𝑖
d𝑡 d𝑡
(39)
d𝒏𝑐𝑖 d𝝋𝑖
= − 𝒓𝑖𝑗 × 𝑩𝑖 .
d𝑡 d𝑡

where 𝒓𝑖𝑗 × is the anti-symmetric matrix of the vector 𝒓𝑖𝑗 . Eliminate the time derivative part, and the above
equation can be written as:

d𝒏𝑞𝑖𝑗 = d𝒏𝑐𝑖 − 𝒓𝑖𝑗 × 𝑩𝑖 d𝝋𝑖 . (40)

So, the partial derivative of 𝒏𝑞𝑖𝑗 to 𝑼𝑖 is:

̅ 𝑖𝑗 = 𝜕𝒏𝑞𝑖𝑗
𝑮 = [𝑰3 −𝒓𝑖𝑗 × 𝑩𝑖 ], (41)
𝜕𝑼𝑇 𝑖

𝜕𝒂 𝜕𝒃𝑇
where 𝜕𝒃𝑇 and 𝜕𝒂
represent the partial derivative of vector 𝒂 to vector 𝒃 in numerator layout, respectively. The

partial derivative of 𝒏𝑞𝑖𝑗 to the minimal coordinate 𝑼 is:

𝜕𝒏𝑞𝑖𝑗 𝜕𝒏𝑞𝑖𝑗 𝜕𝑼𝑖

𝑮𝑖𝑗 = = ̅ 𝑖𝑗 𝑬𝑇𝑈 .
=𝑮 (42)
𝜕𝑼𝑇 𝜕𝑼𝑇
𝑖 𝜕𝑼
𝑇 𝑖

9
The transformation matrix 𝑮 of the entire structure is:
𝑧
𝜕𝒏 𝜕(𝑬𝑛𝑎 𝒏𝒂 +∑𝑚 𝑖
𝑖=1 ∑𝑗=1 𝑬𝑛𝑞𝑖𝑗 𝒏𝑞𝑖𝑗 ) 𝑧
𝑮= 𝜕𝑼𝑇
= 𝜕𝑼T
= 𝑬𝑛𝑎 𝑬𝑇𝑈𝑎 + ∑𝑚
𝑖=1 ∑𝑗=1 𝑬𝑛𝑞𝑖𝑗 𝑮𝑖𝑗 ,
𝑖 (43)

which maps the difference of nodal coordinate 𝒏 to the difference of minimal coordinate 𝑼.

4 Equilibrium equation

4.1 The Lagrangian method

The general form of the Lagrangian equation is:

d ∂𝐿 ∂𝐿
˙
− = 𝑸𝑛𝑝 ,
d𝑡 ∂𝑼 ∂𝑼 (44)

where 𝐿 = 𝑇 − 𝑉 is the Lagrangian function, 𝑇 and 𝑉 are the kinetic energy and potential energy of the system,
𝑸𝑛𝑝 is the non-potential force vector of the general tensegrity structures, 𝑼 is the minimal coordinate of the
system. For the statics problem, the kinetic energy 𝑇 is zero in this study, and we study the potential energy of
the system. For statics problem, the Lagrangian method degenerates to:

∂𝑉
= 𝑸𝑛𝑝 .
∂𝑼 (45)

Note that Eq.(45) is consistent with the principle of stationary total potential energy and the principle of virtual
work. However, using the Lagrangian method to derive the equilibrium equation will make it easy to extend to
the future study of the dynamic problem. It is required in the Lagrangian method to use minimal coordinate as
the variable, which is critical for the derivation. Note that if we use variables with overparameterization like the
Euler parameter, modified Rodrigues parameters, etc., there will be an issue in violation of the constraints of the
variables.

4.2 Energy function

The total potential energy 𝑉 of the tensegrity system with the rigid body is composed of strain potential energy
𝑉𝑒 and gravitational potential energy 𝑉𝑔 :

𝑉 = 𝑉𝑒 + 𝑉𝑔 . (46)

4.2.1 Strain potential energy

There is no deformation in a rigid body, so the strain potential energy for a rigid body is zero. The strain potential
energy is only stored in the axial members:

10
 𝑛𝑖 𝑙
𝑉𝑒 = ∑𝑖=1
𝑒
∫𝑙 𝒕𝑖 d𝑥. (47)
0𝑖

𝜕𝑉𝑒
From the statics equation of traditional tensegrity [34], we can compute the partial derivative of 𝑉𝑒 to 𝑼, 𝜕𝒏
:

𝜕𝑉𝑒 𝜕𝒏𝑇 𝜕𝑉𝑒

= = 𝑮𝑻 (𝑪𝑇 𝒍̂−1 𝒕̂𝑪) ⊗ 𝑰3 𝒏. (48)
𝜕𝑼 𝜕𝑼 𝜕𝒏

4.2.2 Gravitational potential energy

The gravitational potential energy is relative to any member that has mass. In tensegrity with a rigid body, all
axial members, point mass, and rigid body will contribute to gravitational potential energy:

𝑉𝑔 = 𝑉𝑔𝑒 + 𝑉𝑔𝑝 + 𝑉𝑔𝑟 . (49)

The gravitational potential energy corresponding to the axial elements 𝑉𝑔𝑒 is:

𝑛𝑒 𝑥𝑗𝑖 + 𝑥𝑘𝑖
𝑚𝑒𝑖
𝑉𝑔𝑒 = ∑ [𝑎𝑥 𝑎𝑦 𝑎𝑧 ] [𝑦𝑗𝑖 + 𝑦𝑘𝑖 ]
2
𝑖=1 𝑧𝑗𝑖 + 𝑧𝑘𝑖
(50)
𝑛 𝑚𝑒𝑖
= ∑𝑖=1
𝑒
[𝑎𝑥 𝑎𝑦 𝑎𝑧 ] |𝑪𝑖 | ⊗ 𝑰3 𝒏
2

1
= (𝒎𝑇𝑒 |𝑪|) ⊗ [𝑎𝑥 𝑎𝑦 𝑎𝑧 ]𝒏.
2

where 𝑚𝑒𝑖 is the mass of the ith axial element, and 𝒎𝑒 is the mass vector of all axial elements. 𝑎𝑥 , 𝑎𝑦 , 𝑎𝑧 are the
gravitational acceleration in the 𝑋, 𝑌, and 𝑍-axis, respectively. The gravitational potential energy corresponding
to point mass 𝑉𝑔𝑝 is:
𝑛𝑛 𝑥𝑖
𝑉𝑔𝑝 = ∑ 𝑚𝑝𝑖 ⊗ [𝑎𝑥 𝑎𝑦 𝑎𝑧 ] [𝑦𝑖 ]
𝑖=1 𝑧𝑖 (51)

= 𝒎𝑇𝑝 ⊗ [𝑎𝑥 𝑎𝑦 𝑎𝑧 ]𝒏.

where 𝑚𝑝𝑖 is the mass of the ith node, and 𝒎𝑝 is the node mass vector. The gravitational potential energy
corresponding to rigid body 𝑉𝑔𝑟 is:
𝑛𝑞

𝑉𝑔𝑟 = ∑ 𝑚𝑞𝑖 ⊗ [𝑎𝑥 𝑎𝑦 𝑎𝑧 ] 𝒏𝑐𝑖

𝑖=1
(52)
𝒏𝑐1
= 𝒎𝑇𝑞 ⊗ [𝑎𝑥 𝑎𝑦 𝑎𝑧 ] [ ⋮ ]
𝒏𝑐𝑚

11
 = 𝒎𝑇𝑞 ⊗ [𝑎𝑥 𝑎𝑦 𝑎𝑧 ]𝑬𝑇𝑈𝑐 𝑼.

where 𝑚𝑞𝑖 is the mass of the ith rigid body, and 𝒎𝑞 is the mass vector rigid bodies. The partial derivative of 𝑉𝑔
to 𝒏 is:

𝜕𝑉𝑔 𝜕𝒏𝑇 𝜕𝑉𝑔𝑒 𝜕𝑉𝑔𝑚 𝜕𝑉𝑔𝑟

= ( + )+
𝜕𝑼 𝜕𝑼 𝜕𝒏 𝜕𝒏 𝜕𝑼
(53)
1
= {𝑮𝑻 (2 |𝑪|𝑇 𝒎𝑒 + 𝒎𝑝 ) + 𝑬𝑈𝑐 𝒎𝑞 } ⊗ [𝑎𝑥 𝑎𝑦 𝑎𝑧 ]𝑇 = 𝒈,

where 𝒈 is the gravitational force vector.

4.3 Nonlinear equilibrium equation

The statics equation of tensegrity with the rigid body is calculated by the partial derivative of V with respect to
𝑼:

𝜕𝑉 𝜕𝑉𝑒 𝜕𝑉𝑔
= + = 𝑸𝑛𝑝 . (54)
𝜕𝑼 𝜕𝑼 𝜕𝑼

Substitute the Eq.(48) and Eq.(53) into Eq.(54), we will have:

𝑮𝑻 (𝑪𝑇 𝒍̂−1 𝒕̂𝑪) ⊗ 𝑰3 𝒏 = 𝑸𝑛𝑝 − 𝒈. (55)

Eq.(55) is the static equilibrium equation of the general tensegrity system with rigid bodies. The second part
(𝑪𝑇 𝒍̂−1 𝒕̂𝑪) ⊗ 𝑰3 𝒏 is the collection of inner force of members in nodes, which is identical to 𝑲𝒏 in traditional
tensegrity structure [34]. Note that (𝑪𝑇 𝒍̂−1 𝒕̂𝑪) ⊗ 𝑰3 is a nonlinear function of nodal coordinate, so Eq. (55) is
nonlinear. The first part 𝑮𝑻 transforms the nodal force from the node space to body space, which is identity to
the generalized force. Eq.(55) can be written into a simple form:

𝑲𝑟 𝒏 = 𝑸𝑛𝑝 − 𝒈, (56)

where 𝑲𝑟 is the stiffness matrix of general tensegrity with nodal coordinate vector 𝒏 as the variable:

𝑲𝑟 = 𝑮𝑻 (𝑪𝑇 𝒍̂−1 ̂𝒕𝑪) ⊗ 𝑰3 . (57)

The right part of Eq. (54) is the generalized force 𝑸𝑛𝑝 , which can be calculated by using the transformation
matrix [32]:
𝑛𝑒 𝑛𝑒
𝜕𝒏𝑇 𝜕𝒏𝑇𝑐𝑖 𝜕𝝎𝑇𝑖
𝑸𝑛𝑝 = 𝒇+ ∑ 𝒇𝑐𝑖 + ∑ 𝒎𝑐𝑖
𝜕𝑼 𝜕𝑼 𝜕𝑼̇ (58)
𝑖=1 𝑖=1

𝑛𝑒 𝑛𝑒
= 𝑮𝑻 𝒇 + ∑𝑖=1 𝑬𝑈𝑐𝑖 𝒇𝑐𝑖 + ∑𝑖=1 𝑬𝑈𝝋𝑖 𝑩𝑇𝑖 𝒎𝑐𝑖

12
 = 𝑮𝑻 𝒇 + 𝑬𝑈𝑐 𝒇𝑐 + 𝑬𝑈𝝋 𝑩𝑇 𝒎𝑐 ,

where 𝒇 is the non-potential external force vector exerted on the tensegrity node. 𝒇𝑐𝑖 and 𝒎𝑐𝑖 is the total force
and moment exerted on the 𝑖th rigid body. 𝒇𝑐 and 𝒎𝑐 are the collection of force and moment of all rigid bodies.

𝒇𝑐1 𝒎𝑐1
𝒇𝑐 = [ ⋮ ] , 𝒎𝑐 = [ ⋮ ]. (59)
𝒇𝑐𝑚 𝒎𝑐𝑚

𝑩 matrix is defined as:

𝑩1
𝑩=[ ⋱ ]. (60)
𝑩𝑚

4.4 Linearized equilibrium equation

4.4.1 Linearized equilibrium equation with minimal coordinate as the variable

Using Taylor expansion of Eq. (56) about a configuration 𝒏𝑖 in the ith iteration step, we have the linearized
equilibrium equation:

𝑲𝑟 |𝒏𝑖 𝒏𝑖 + 𝑲 𝑇𝑟 (𝑼𝑖+1 − 𝑼𝑖 ) = 𝑸𝑛𝑝 − 𝒈, (61)

where 𝑲 𝑇𝑟 is the tangent stiffness matrix of the structure, 𝑼𝑖 is the minimal coordinate corresponding to 𝒏𝑖 .
𝑲𝑟 |𝒏𝑖 is the stiffness matrix in 𝒏𝑖 configuration. By solving Eq.(61), we can obtain a new configuration 𝑼𝑖+1 in
the i+1 iteration step, which is closer to the equilibrium configuration. The out-of-balance forces of the system
is defined as:

𝑷𝑖 = 𝑸𝑛𝑝 − 𝒈 − 𝑲𝑟 |𝒏𝑖 𝒏𝑖 . (62)

The difference of the minimal coordinate can be simply computed by:

d𝑼𝑖 = 𝑲−1 𝑖
𝑇𝑟 𝑷 . (63)

The above three equations can be used in solving nonlinear equilibrium equations based on an iteration method.

4.4.2 Linearized equilibrium equation in terms of the member force

The Eq.(55) can be written linearly in terms of the member force 𝒕:

𝑨𝑟 𝒕 = 𝑸𝑛𝑝 − 𝒈, (64)

where 𝑨𝑟 ∈ ℝ𝑛𝑈×𝑛𝑒 is the equilibrium matrix for tensegrity with rigid bodies:

𝑨𝑟 = 𝑮𝑻 𝑨2 . (65)

13
where 𝑨2 is the equilibrium equation of traditional tensegrity [34]:

𝑨2 = 𝑪𝑇 ⊗ 𝑰3 b. d. (𝑯). (66)

where b. d. (𝑯) is the block diagonal matrix of 𝑯. Note that the equilibrium matrix for tensegrity with rigid
bodies 𝑨𝑟 is identical to the 𝑪 matrix in Wang et al. [25]. The singular value decomposition of the equilibrium
matrix 𝑨𝑟 reveals the self-stress mode and mechanism mode of the structure [35]:

𝑾𝟐 ] [𝚺0 𝟎 𝑽1𝑇
𝑨𝑟 = 𝑾𝚺𝑽𝑇 = [𝑾𝟏 ] [ ], (67)
𝟎 𝟎 𝑽𝑇2

where 𝑾 ∈ 𝑹𝑛𝑈 ×𝑛𝑈 , and 𝑽 ∈ 𝑹𝑛𝑒×𝑛𝑒 are orthogonal matrices. Let 𝑟 = rank(𝑨𝑟 ) be the rank of 𝑨𝑟 . 𝑽1 ∈
𝑹𝑛𝑒 ×r , 𝑽2 ∈ 𝑹𝑛𝑒 ×(𝑛𝑒 −r) is respectively the row space and null space of 𝑨𝑟 , and 𝑾1 ∈ 𝑹𝑛𝑈×r , 𝑾2 ∈ 𝑹𝑛𝑈 ×(𝑛𝑈 −r)
is respectively the column space and left null space of 𝑨𝑟 . 𝑨𝑟 𝑽2 = 𝟎 and 𝑨𝑇𝑟 𝑾2 = 𝟎, 𝑽2 and 𝑾2 are the self-
stress mode and mechanism mode of the tensegrity structure, respectively.

4.4.3 Compatibility equation

The compatibility equation is the relation between d𝑼 and d𝒍 that guarantees the structural deformations are
physically valid. The compatibility equation of the structure is:

𝑩𝑟 d𝑼 = d𝒍, (68)

where 𝑩𝒓 ∈ ℝ𝑛𝑒×𝑛𝑈 is the compatibility matrix:

𝜕𝒍 𝜕𝒍 𝜕𝒏
𝑩𝑟 = = = 𝑨𝑇2 𝑮. (69)
𝜕𝑼𝑇 𝜕𝒏𝑇 𝜕𝑼𝑇

It can be found that the compatibility matrix is the transpose of the equilibrium matrix:

𝑩𝑟 = 𝑨𝑇𝑟 . (70)

This can also be proved by the principle of virtual work.

4.5 Tangent stiffness matrix

Refer to the derivation of tangent stiffness in Chen and Jiang [26], the tangent stiffness matrix of the general
tensegrity with the rigid body is:

𝜕(𝑩𝑇
𝑟 𝑡) 𝜕𝒕 𝜕𝑩𝑇
𝑲 𝑇𝑟 = 𝜕𝑼𝑇
= 𝑩𝑇𝑟 𝜕𝑼𝑇 + 𝜕𝑼𝑟𝑇 𝒕 = 𝑲𝐸 + 𝑲𝐺 . (71)

The first part of Eq.(71) is the material stiffness 𝑲𝐸 caused by the difference of member force:

𝜕𝒕 𝜕𝒍
𝑲𝐸 = 𝑩𝑇𝑟 ̂ 𝑩𝑟 = 𝑨𝑟 𝒌
= 𝑩𝑇𝑟 𝒌 ̂𝑨𝑇𝑟 , (72)
𝜕𝒍𝑇 𝜕𝑼𝑇

14
 ̂𝑨
where 𝒌 = 𝑬 ̂ 𝒍−1
0 is the stiffness of the axial members. The second part of Eq.(71) is the geometry stiffness 𝑲𝐺

caused by the difference of structural shape:

𝜕𝑩𝑇 𝑛𝑠
𝑲𝐺 = 𝑟
𝜕𝑼𝑇
𝒕 = 𝛀𝑇 𝒕 = ∑𝑖=1 𝛀𝑇𝑖 𝒕𝑖 . (73)

where the Hessian matrix 𝛀 ∈ ℝ𝑛𝑒×𝑛𝑈 ×𝑛𝑈 is expressed as:

𝜕𝑩𝑟 𝑇
𝛀= = [𝛀1𝑇 ⋯ 𝛀𝑇𝑖 ⋯ 𝛀𝑇𝑛𝑒 ] . (74)
𝜕𝑼

𝜕𝑩𝒓𝑖
where 𝛀𝑖 = 𝜕𝑼
∈ ℝ𝑛𝑈 ×𝑛𝑈 is the 𝑖th member’s Hessian matrix. Note that the explicit formulation of 𝛀 is vital

to calculate the geometry stiffness matrix. Fortunately, 𝛀𝑖 can be obtained by calculating and comparing two
equivalent expressions of the ith cable’s acceleration 𝑙𝑖̈ . The Eq. (68) is equivalent to:

𝑙𝑖̇ = 𝑩𝑟𝑖 𝑼̇, (75)

𝜕𝑩𝑟𝑖 𝜕𝑼𝑇 𝜕𝑩𝑟𝑖

Using = = 𝑼̇𝑇 𝛀𝑖 , the 𝑖th cable’s acceleration 𝑙𝑖̈ is:
𝜕𝑡 𝜕𝑡 𝜕𝑼

𝜕𝑩
𝑙𝑖̈ = 𝑩𝑟𝑖 𝑼̈ + 𝜕𝑡𝑟𝑖 𝑼̈ = 𝑩𝑟𝑖 𝑼̈ + 𝑼̇𝑇 𝛀𝑖 𝑼̇. (76)

From the derivation in Appendix, the ith cable’s acceleration 𝑙𝑖̈ is expressed as:

𝑷ℎ𝑛𝑖
𝑙𝑖̈ = 𝑩𝑟𝑖 𝑼̈ + 𝑼̇𝑇 (𝑮𝑇 (𝑪𝑇𝑖 ⊗ 𝑰3 ) 𝑙𝑖
(𝑪𝑖 ⊗ 𝑰3 )𝑮 + 𝑭𝑖 )𝑼̇. (77)

Comparing Eq. (76) with Eq. (77), the matrix 𝛀𝑖 is written as:

𝑷ℎ𝑛𝑖
𝛀𝑖 = 𝑮𝑇 (𝑪𝑇𝑖 ⊗ 𝑰3 ) (𝑪𝑖 ⊗ 𝑰3 )𝑮 + 𝑭𝑖 . (78)
𝑙𝑖

where 𝑷ℎ𝑛𝑖 = 𝑰3 − 𝒉𝑛𝑖 𝒉𝑇𝑛𝑖 ∈ ℝ3×3 denotes the projector to the plane with the normal vector 𝒉𝑛𝑖 , in which
𝒉𝑖
𝒉 𝑛𝑖 = 𝑙𝑖
is the ith cable’s unit vector. From the derivation in Appendix, the matrix 𝑭𝑖 ∈ ℝ𝑛𝑈×𝑛𝑈 is written as:

𝑧𝑖 𝟎 𝟎
𝑭𝑖 = ∑𝑚 𝑇
𝑗=1 ∑𝑘=1 𝑬𝑈𝑗 [𝟎 𝑩𝑇 𝒛× 𝒓× 𝑩 ] 𝑬𝑈𝑗 . (79)
𝑗 𝑖𝑗𝑘 𝑗𝑘 𝑗

in which 𝒛𝑖𝑗𝑘 ∈ ℝ3 is:

𝑇
𝒛𝑖𝑗𝑘 = (𝒉𝑇𝑛𝑖 (𝑪𝑇𝑖 ⊗ 𝑰3 )𝑬𝑛𝑞𝑗𝑘 ) . (80)

Note that the tangent stiffness is a general form of classical tensegrity. That is, if there is no rigid body, the
tangent stiffness will degenerate to a classical tensegrity [34]. Also, note that the above derivation is generally
consistent with the formulation in Chen and Jiang [26]. The difference is that the proposed formulation in this

15
paper is capable of considering free and pinned tensegrity nodes in the general tensegrity system, and the use of
location matrix makes the formulation in Eq.(79) be expressed in a more simple and neat form.

5 Form-finding of tensegrity systems with rigid bodies

In this section, we formulate the form-finding method for tensegrity systems with rigid bodies. Three numerical
examples are carried out to illustrate the accuracy and efficiency of the proposed form-finding method.

5.1 Form-finding method

5.1.1 Form-finding procedure

The form-finding method is basically solving the nonlinear equilibrium equation. However, the self-equilibrated
tensegrities lacking proper constraints have several problems in solving its equilibrium equation [36]. Firstly,
the rigid body mode will lead to a singular tangent stiffness matrix. Newton's method is not able to solve the
equation with a singular Hessian matrix. Secondly, the tangent stiffness matrix may have a negative eigenvalue,
and the result of solving the nonlinear equilibrium equation will converge to an unstable equilibrium
configuration. To ensure the result is stable equilibrium, modification of the tangent stiffness matrix to positive
definite is necessary. Thirdly, an appropriate optimization objective needs to be defined to guarantee that the
result approaches the equilibrium configuration. The form-finding procedure consists of the following main steps,
as shown in Figure 2.

Inputs:

(1) Specify the basic data of a tensegrity system with rigid bodies, including the minimal coordinate 𝑼0 ,
connectivity matrix 𝑪, axial stiffnesses vector 𝑬, cross-section area vector 𝑨, rest length vector 𝒍0, location
matrix 𝑬na , 𝑬𝑛𝑏 , 𝑬𝑛𝑞𝑖𝑗 , 𝑬Ua , 𝑬𝑈𝑖 , 𝑬𝑈𝜑 , coefficient 𝑢 and ε. Compute the nodal coordinate 𝒏0, stiffness matrix

𝑲𝑟 |𝒏0 , out-of-balance forces 𝑷0 in the initial configuration.

Iteration:

(2) Compute the tangent stiffness 𝑲 𝑇𝑟 for the structure in the current configuration. Compute the minimal
eigenvalue of the 𝑲 𝑇𝑟 as λ.
(3) Check whether the tangent stiffness matrix is positive definite or not. Use the method in Section 5.1.2 to
modify the stiffness matrix such that it is positive definite.
(4) Solve the difference of minimal coordinate d𝑼𝑖 , employ the line search algorithm in Section 5.1.3 to
calculate the updated minimal coordinate 𝑼𝑖 .
(5) Calculate the nodal coordinate 𝒏𝑖 , stiffness matrix 𝑲𝑟 |𝒏𝑖 and out-of-balance forces 𝑷𝑖 . Check whether the
current configuration is in equilibrium or not. If not, set 𝑖 ← 𝑖 + 1 and go to step (2).

16
(6) Terminate the iteration when an equilibrium configuration has been obtained.

Figure 2 Flow chart of the form-finding algorithm.

5.1.2 Modification of tangent stiffness matrix

To guarantee the form-finding result converges to a stable equilibrium. The positive definiteness of the tangent
stiffness matrix 𝑲 𝑇𝑟 should be examined and modified. For the configuration 𝑼𝑖 at an iteration step, if the
minimal eigenvalue of the tangent stiffness matrix λ is negative, a sufficiently large identity matrix (|λ| + 𝑢)𝑰
̃ 𝑇𝑟 , where 𝑢 is a positive coefficient to
will be added to 𝑲 𝑇𝑟 to obtain the modified tangent stiffness matrix 𝑲
guarantee the modified tangent stiffness matrix is not seriously ill. Otherwise, 𝑢𝑰 will be added to the tangent
stiffness matrix:

𝑲 + (𝑢 + ‖𝜆‖)𝑰, 𝜆 < 0
̃ 𝑇 = { 𝑇𝑟
𝑲 . (81)
𝑟 𝑲 𝑇𝑟 + 𝑢𝑰, 𝜆 > 0

From experience, in this paper, we set 𝑢 = 0.1. Using the modified tangent stiffness matrix, the increment of the
generalized coordinate vector d𝑼 can be obtained from Eq. (63):

17
 ̃ −1
d𝑼𝑖 = 𝑲 𝑖
𝑇𝑟 𝑷 . (82)

5.1.3 Line search algorithm

To increase the convergence speed of solving the nonlinear equilibrium equation. We use a line search algorithm
[37,36] in each iteration step to minimize the total potential energy of the system. In the ith step, we update the
minimal coordinate vector 𝑼𝑖 from that in step 𝑖 − 1 by:

𝑼𝑖 = 𝑼𝑖−1 + 𝑥d𝑼𝑖 . (83)

where the coefficient 𝑥 is determined by the following optimization problem of single-variable function on the
fixed interval:

min 𝑉(𝑥)
(84)
s. t. 0 < 𝑥 ≤ 1.

Given 𝑼𝑖 , the nodal coordinate vector 𝒏𝑖 can be calculated by Eqs.(9), (24), and (28). And the total potential
energy can be calculated by Eqs.(46) to (52). The line search algorithm can be simply implemented by the
‘fminbnd’ function in MATLAB.

5.2 Numerical examples

In this section, four examples are studied to demonstrate the accuracy and efficiency of the proposed form-
finding method for tensegrity with rigid bodies. Different examples are carefully chosen to represent generalized
tensegrity with one or multiple rigid bodies, with or without free nodes and pinned nodes. In these examples, the
equilibrium configurations and prestress are tuned by varying the rest length of the strings in the structure. The
tangent modulus and cross-sectional area of the strings in all the examples are set to be 7.6 × 1010 Pa and
1 × 10−4 m2.

5.2.1 Patio shade cover

This example presents a structure composed of a rigid triangle piece, five strings, a free node, and four pinned
nodes. The index of nodes and elements are marked in black numbers and blue numbers in circles, respectively,
as shown in Figure 3. And Figure 3 is the initial configuration of the generalized tensegrity. To generate the
prestress of the structure, the rest length of strings is set to be 0.3 times the present length in the initial
configuration, which is 𝒍0 = 0.3𝒍. Figure 4 gives the equilibrium configuration of the form-finding result.

18
 Figure 3. The initial configuration Figure 4. The equilibrium configuration

The nodal coordinate matrix 𝑵 ∈ ℝ3×8 in the equilibrium configuration in the form of Eq. (3) is given as:

0.8955 −0.8365 0.0310 2.0000 −2.0000 0.0094 1.0000 −1.0000 0.0300

𝑵 = [0.3221 0.3204 −1.1787 1.0000 1.0000 −2.3218 −3.0000 −3.0000 −0.1787]. (85)
0 0 0 0 0 0 0 0 0

From Eq. (13), the connectivity matrix 𝑪 ∈ ℝ5×8 in initial configuration is:

−1 0 0 1 0 0 0 0 0
0 −1 0 0 1 0 0 0 0
𝑪𝑠 = 0 0 −1 0 0 1 0 0 0. (86)
0 0 0 0 0 −1 1 0 0
[0 0 0 0 0 −1 0 1 0]

From Eq.(15), one can get the structure element matrix 𝑯 ∈ ℝ3×5:

1.1045 −1.1635 −0.0216 0.9906 −1.0094

𝑯 = [0.6779 0.6796 −1.1430 −0.6782 −0.6782]. (87)
0 0 0 0 0

From Eq. (65), the equilibrium matrix for tensegrity with rigid bodies 𝑨𝑟 ∈ ℝ9×5 can be calculated:

0 0 −0.0189 −0.8251 0.8300

0 0 −0.9998 0.5649 0.5577
0 0 0 0 1
−0.8523 0.8635 0.0189 0 0
𝑨𝑟 = −0.5231 −0.5044 0.9998 0 0 . (88)
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
[−0.0259 0.0060 0.0199 0 0 ]

Singular value decomposition of equilibrium matrix reveals the rank of 𝑨𝑟 is 𝑟 = 4. That is to say, the structure
has s = 5 − r = 1 self-stress mode and m = 9 − 𝑟 = 5 mechanism modes. The null space of the equilibrium
matrix 𝑨𝑟 gives the self-stress mode 𝑽2 of the system:

19
 𝑽2 = [−0.4678 −0.4514 −0.4725 −0.4166 −0.4250]𝑇 . (89)

The left null space of the equilibrium matrix 𝑨𝑟 gives the mechanism modes 𝑾2 of the system:

0 0 0 0 0
0 0 0 0 0
−0.6024 0 0 0 −0.7982
0.0147 0 0 0 −0.0111
𝑾2 = 0.0156 0 0 0 −0.0118 . (90)
0 1.0000 0 0 0
0 0 1.0000 0 0
0 0 0 1.0000 0
[−0.7970 0 0 0 0.6021 ]

Each column of 𝑾2 represent a mechanism mode. The five mechanism mode shapes are plotted in Figure 5,
where the dashed line and solid lines are the equilibrium configuration and the deformed shape of the mechanism,
respectively. The 1st, 3rd, and 4th mechanism modes correspond to the rotation motion of the rigid body about
the X, 𝑌, 𝑍-axis. The 2nd mechanism mode contains the translation motion of the rigid body in the 𝑍-axis, and
the 5th mode contains the translation of free node in the 𝑍-axis, translation of mass center in 𝑋, 𝑌-axis, and
rotation of rigid body by 𝑍-axis.

Figure 5 The five mechanism modes of the structure.

The mechanism mode is the null space of the material stiffness matrix which means there is no elongation of the
axial member in the mechanism mode. For tensegrity systems, the mechanism mode can be stiffened by prestress,
and the stability of the system can be checked by the product force [38,39] or by the positive-definite of tangent
stiffness matrix [25,40]. The eigenvalue of the tangent stiffness matrix 𝑲 𝑇𝑅 is plotted in Figure 6, we can see
that all the eigenvalues of the tangent stiffness matrix are positive, which means the mechanism mode is stiffened
by prestress.

20
 Figure 6 Eigenvalue of the tangent stiffness matrix.

The deformed shape corresponding to the first four eigenvalues of the tangent stiffness matrix is plotted in Figure
7, where the dotted line is the equilibrium configuration, and the solid line is the deformed shape. As we can see,
the 1st and 4th mode shapes contain out-of-plane deformation, while the 2nd and 3rd mode shapes contain pure
planer deformation.

Figure 7 Deformed shapes of the modes corresponding to the first four eigenvalues.

5.2.2 Tensegrity table

This example presents a self-equilibrated tensegrity table composed of two rigid bodies and four strings. Figure
8 is the initial configuration. The rest length of strings is set to be 0.3 times of present length, which is 𝒍0 = 0.3𝒍,
to generate prestress of the structure. Figure 8 shows the equilibrium configuration of the form-finding result.

21
 Figure 8 The initial configuration of the tensegrity table. Figure 9 The equilibrium configuration of the tensegrity
table.

The prestress mode of the equilibrium matrix is:

𝑽2 = [0.2887 0.2887 0.2887 0.8660]𝑇 . (91)

The first three values reveal that the forces of the three long strings are the same. And the fourth value indicates
that the inner force of the short string is three times of the long string at an equilibrium state.

Figure 10 Eigenvalues of the tangent stiffness matrix.

Figure 10 is the eigenvalue of tangent stiffness. The first six eigenvalues correspond to the rigid body modes of
the structure. Figure 11 shows the mode shape of the tensegrity table, where mode 6 is a pure rotational mode
that has zero stiffness. The 7th mode is the most flexible one, which involves the relative rotation of two rigid
bodies around the Z-axis.

22
 Figure 11 The mode shape of the tangent stiffness matrix.

5.2.3 Spherical tensegrity

This example presents a spherical tensegrity composed of multiple rigid bodies in which all nodes lie on the
vertices of a regular polyhedron. Truncated tensegrity is the simplest way to build spherical tensegrities, and
there are a few studies about this topic [41,36,42,43]. In this example, we propose a novel method to build
spherical tensegrity with rigid bodies and study the equilibrium condition of the structure. Here we use the
tetrahedron as an example to illustrate the step-by-step procedure to generate a spherical tensegrity with a rigid
body, and the equilibrium configuration as well as the member force of all the other regular polyhedrons
tensegrity with rigid bodies.

In Figure 12, four rigid bodies are initially placed in the plane of the tetrahedron, and rigid bodies nodes are
placed in the vertices of the tetrahedron. Each rigid body is rotated by an angle φ about the normal line of the
plane to generate a new shape with 12 nodes, as shown in Figure 13.

Figure 12 A tetrahedron build with rigid bodies. Figure 13 Rotation of the rigid bodies.

If we connect the nodes of the rigid bodies in the initial configuration, there will be 12 truncating-edge strings
and 6 vertical strings, as in Figure 14. To prestress the spherical tensegrity, the rest length of the vertical strings
is set to [0.1,0.8] times its present length while truncating-edge strings are identical to its present length. The
form-finding result of a truncated tetrahedral generalized tensegrity is shown in Figure 15. The force density of

23
truncating-edge strings and vertical strings are respectively 𝑞𝑡 and 𝑞𝑣 . We can observe that the force density of
both truncating-edge strings and vertical strings increases as the rest length of vertical strings decreases.

Figure 14 Initial configuration of a tetrahedron tensegrity Figure 15. Form-finding solution of tetrahedron generalized
with rigid bodies tensegrity

The form-finding result for other regular polyhedron shapes, including hexahedron, octahedral, dodecahedral,
and icosahedral generalized tensegrities, are shown in Figure 16 to Figure 19.

Figure 16 Form-finding solution of hexahedron generalized Figure 17 Form-finding solution of octahedral generalized
tensegrities. tensegrities.

24
 Figure 18 Form-finding solution of dodecahedral generalized Figure 19 Form-finding solution of icosahedral generalized
tensegrities. tensegrities.

5.2.4 Tensegrity spine

As the last example, we study a tensegrity spine [44,45] composed of multiple rigid bodies. Figure 20 is the
initial configuration of the tensegrity spine. The tensegrity spine is composed of 10 rigid body units, and the 10
rigid bodies are connected by four groups of vertical side strings and nine groups of diagonal strings.

Figure 20 The initial configuration of the tensegrity spine.

The rest length of all the diagonal strings is set to 0.9 times the present length. The rest length of the three groups
of vertical side strings is set to 0.9 times the present length, while the rest length of the other group of vertical
side strings is set to 0.6 times the present length. The equilibrium configuration calculated by the form-finding
method is shown in Figure 21.

25
The rest length of two groups of vertical side strings in the opposite positions is 0.9 times the length in the initial
shape, while the rest length of the other two groups of vertical side strings varies from 0.5 to 1.1 times the length
in the initial shape. The equilibrium configuration calculated by the form-finding method is shown in Figure 22.

Figure 21 C-shape, achieved by changing the rest length Figure 22 S-shape, achieved by changing the
of the strings on one side linearly. rest length of the strings on two sides
sinusoidally but with different phases.

6 Conclusions
During the past few decades, pure bar-string network tensegrity has shown its great strength in designing efficient
structures in many aspects. However, to embrace a much more general problem of system design using the
tensegrity paradigm, rigid bodies must be included. Aiming at extending the ability to analyze rigid body
tensegrities with analytical tools, this paper formulates the nonlinear equilibrium equation of the rigid body
tensegrity in an explicit form in terms of the minimal coordinate. To get the insight of each structure member,
we derived its equivalent form, which is a linear equation in terms of the force vector. Then, we also provide the
compatibility equation and tangent stiffness matrix of the system for stability analysis. Finally, based on the
equilibrium and stiffness equations, an efficient form-finding method of the rigid body tensegrity is given. In the
proposed form-finding method, modification of tangent stiffness matrix and line search algorithm is used to
guarantee the result to fast converge to a stable equilibrium configuration. It is also shown that without rigid
bodies, the nonlinear equilibrium equations of the general tensegrity degenerate to the ones of the traditional
tensegrity. Four numerical examples are given to prove the accuracy and efficiency of the developed principles.
Results show that the developed principles are capable of dealing with form-finding from a non-equilibrium
state, finding the prestress and mechanism modes, and conducting stiffness studies.

26
Acknowledgment
The research was supported by the Foundation of Key Laboratory of Space Structures of Zhejiang Province
(Grant No. 202102) and the National Natural Science Foundation of China (Grant No. 51878600).

Appendix

Derivation of the cable acceleration

From Eq.(16), the 𝑖th cable’s velocity is:

𝒉𝑇
𝑙𝑖̇ = 𝑖
𝑙𝑖
𝒉̇𝑖 = 𝒉𝑇𝑛𝑖 𝒉̇𝑖 , (92)

𝒉𝑖
where 𝒉𝑛𝑖 = 𝑙𝑖
is the 𝑖th cable’s unit vector. The 𝑖th cable’s acceleration is:

𝑇
𝑙𝑖̈ = 𝒉̇𝑛𝑖 𝒉̇𝑖 + 𝒉𝑇𝑛𝑖 𝒉̈𝑖 . (93)

𝒉̇𝑛𝑖 is the time derivative of the 𝑖th cable’s unit vector, which can be derived as:

𝒉̇𝑖 𝑙𝑖 −𝒉𝑇 ̇
𝑛𝑖 𝒉𝑖 𝒉𝑖 (𝑰3 −𝒉𝑇
𝑛𝑖 𝒉𝑛𝑖 ) 𝑷ℎ𝑛𝑖
𝒉̇𝑛𝑖 = = 𝒉̇𝑖 = 𝒉̇𝑖 , (94)
𝑙𝑖2 𝑙𝑖 𝑙𝑖

where 𝑷ℎ𝑛𝑖 = 𝑰3 − 𝒉𝑛𝑖 𝒉𝑇𝑛𝑖 ∈ ℝ3×3 is a symmetric matrix. Therefore, using Eqs.(43), (14), and (94), the first
term of Eq.(93) can be rewritten as:

𝑷𝑇
ℎ𝑛 𝑷ℎ𝑛
𝒉̇𝑇𝑛𝑖 𝒉̇𝑖 = 𝒉̇𝑇𝑖 𝑖
𝒉̇𝑖 = 𝑼̇𝑇 𝑮𝑇 (𝑪𝑇𝑖 ⊗ 𝑰3 ) 𝑖 (𝑪𝑇𝑖 ⊗ 𝑰3 )𝑮𝑼̇. (95)
𝑙𝑖 𝑙𝑖

The acceleration of the 𝑘th node on the 𝑗th rigid is:

𝒏̈ 𝑞𝑗𝑘 = 𝒏̈ 𝑐𝑗 + 𝝎̇𝑗 × 𝒓𝑗𝑘 + 𝝎𝑗 × (𝝎𝑗 × 𝒓𝑗𝑘 ). (96)

According to Eqs.(9) and (14), the second term of Eq.(93) is:

𝑧
𝒉𝑇𝑛𝑖 𝒉̈𝑖 = 𝒉𝑇𝑛𝑖 (𝑪𝑇𝑖 ⊗ 𝑰3 )𝒏̈ = 𝒉𝑇𝑛𝑖 (𝑪𝑇𝑖 ⊗ 𝑰3 )(𝑬𝑛𝑎 𝒏̈𝑎 + ∑𝑚
𝑗=1 ∑𝑘=1 𝑬𝑛𝑞𝑗𝑘 𝒏̈ 𝑞𝑗𝑘 ).
𝑖
(97)

Substitute Eq.(96) into the second term of Eq.(97), we have:

× 𝒏̈ 𝑐𝑗 𝟎3×3 𝟎3×3 𝒏̇𝑐𝑗

𝒉𝑇𝑛𝑖 (𝑪𝑇𝑖 ⊗ 𝑰3 )𝒏̈𝑞𝑗𝑘 = 𝒛𝑇𝑖𝑗𝑘 [𝑰3 𝑇
(−𝒓𝑗𝑘 ) ] [ 𝝎̇ ] + [𝒏̇𝑐𝑗 𝝎̇𝑗 ] [ × ×][
𝒓𝑗𝑘 𝝎𝑗 ]
=
𝑗 𝟎3×3 𝒛𝑖𝑗𝑘

× 𝑇 𝟎 𝟎 𝑇
𝒛𝑇𝑖𝑗𝑘 [𝑰3 (−𝒓𝑗𝑘 ) 𝑩𝑗 ]𝑼̈𝑗 + 𝑼̇𝑗 [𝟎 𝑩𝑇 𝒛× 𝒓× 𝑩 ] 𝑼̇𝑗 = 𝒉𝑛𝑖 (𝑪𝑇𝑖 ⊗ 𝑰3 )𝑬𝑛𝑞𝑗𝑘 𝑮
̅ 𝑗𝑘 𝑬𝑇𝑈𝑗 𝑼̈ + (98)
𝑗 𝑖𝑗𝑘 𝑗𝑘 𝑗

𝑇 𝟎 𝟎
𝑼̇ 𝑬𝑈𝑗 [𝟎 𝑩𝑇𝑗 𝒛× × ] 𝑬𝑇𝑈𝑗 𝑼̇,
𝑖𝑗𝑘 𝒓𝑗𝑘 𝑩𝑗

27
where 𝒛𝑖𝑗𝑘 is:

𝑇
𝒛𝑖𝑗𝑘 = (𝒉𝑇𝑛𝑖 (𝑪𝑇𝑖 ⊗ 𝑰3 )𝑬𝑛𝑞𝑗𝑘 ) . (99)

Substitute Eq.(98) into Eq.(97), 𝒉𝑇𝑛𝑖 𝒉̈𝑖 can be rewritten explicitly with 𝑼̇ and 𝑼̈:

𝑧𝑖
𝒉𝑇𝑛𝑖 𝒉̈𝑖 = 𝒉𝑇𝑛𝑖 (𝑪𝑇𝑖 ⊗ 𝑰3 ) (𝑬𝑛𝑎 𝑬𝑇𝑈𝑎 + ∑𝑚 ̅ 𝑇 ̈
𝑗=1 ∑𝑘=1 𝑬𝑛𝑞𝑗𝑘 𝑮𝑗𝑘 𝑬𝑈𝑗 ) 𝑼 +

𝟎 𝟎 (100)
𝑧𝑖 ̇𝑇 𝑇 ̇ ̈ ̇𝑇 ̇
∑𝑚
𝑗=1 ∑𝑘=1 𝑼 𝑬𝑈𝑗 [𝟎 𝑩𝑇 𝒛× 𝒓× 𝑩 ] 𝑬𝑈𝑗 𝑼 = 𝑩𝑟𝑖 𝑼 + 𝑼 𝑭𝑖 𝑼,
𝑗 𝑖𝑗𝑘 𝑗𝑘 𝑗

where 𝑭𝑖 ∈ ℝ𝑛𝑈×𝑛𝑈 is:

𝑧𝑖 𝟎 𝟎
𝑭𝑖 = ∑𝑚 𝑇
𝑗=1 ∑𝑘=1 𝑬𝑈𝑗 [𝟎 𝑩𝑇 𝒛× 𝒓× 𝑩 ] 𝑬𝑈𝑗 , (101)
𝑗 𝑖𝑗𝑘 𝑗𝑘 𝑗

With the Eqs.(94)-(101), Eq.(93) can be rewritten as:

𝑷ℎ𝑛𝑖
𝑙𝑖̈ = 𝑩𝑟𝑖 𝑼̈ + 𝑼̇𝑇 (𝑮𝑇 (𝑪𝑇𝑖 ⊗ 𝑰3 ) 𝑙𝑖
(𝑪𝑖 ⊗ 𝑰3 )𝑮 + 𝑭𝑖 )𝑼̇. (102)

Compare Eq.(76) with Eq.(102), the matrix 𝛀𝑖 is:

𝑷ℎ𝑛𝑖
𝛀𝑖 = 𝑮𝑇 (𝑪𝑇𝑖 ⊗ 𝑰3 ) (𝑪𝑖 ⊗ 𝑰3 )𝑮 + 𝑭𝑖 . (103)
𝑙𝑖

References:
[1] Fuller, RB. Synergetics: explorations in the geometry of thinking. Estate of R. Buckminster Fuller 1982.
[2] Lalvani, H. Origins of tensegrity: views of emmerich, fuller and snelson. International Journal of Space Structures.
1996;11(1): 28p.
[3] Wang, Y, Xu, X, Luo, Y. Topology-finding of tensegrity structures considering global stability condition. J Struct
Eng. 2020;146(12).
[4] Wang, Y, Xu, X, Luo, Y. A unifying framework for form-finding and topology-finding of tensegrity structures.
Comput Struct. 2021;247: 106486.
[5] Roth, JK, McCarthy, TJ. Optimizing compressive load capacity for differing tensegrity geometries. Comput Struct.
2021;249: 106523.
[6] Wang, Y, Xu, X, Luo, Y. Minimal mass design of active tensegrity structures. Eng Struct. 2021;234: 111965.
[7] Feng, Y, Yuan, XF, Samy, A. Analysis of new wave-curved tensegrity dome. Eng Struct. 2022;250.
[8] Ma, S, Chen, M, Skelton, RE. Design of a new tensegrity cantilever structure. Compos Struct. 2020;243: 112188.
[9] Liu, Y, Bi, Q, Yue, X, Wu, J, Yang, B, Li, Y. A review on tensegrity structures-based robots. Mech Mach Theory.
2022;168: 104571.
28
[10] Paul, C, Valero-Cuevas, FJ, Lipson, H. Design and control of tensegrity robots for locomotion. Ieee T Robot.
2006;22(5): 944-957.
[11] Veuve, N, Sychterz, AC, Smith, IFC. Adaptive control of a deployable tensegrity structure. Eng Struct. 2017;152:
14-23.
[12] Kan, Z, Peng, H, Chen, B, Zhong, W. Nonlinear dynamic and deployment analysis of clustered tensegrity
structures using a positional formulation fem. Compos Struct. 2018;187.
[13] Hrazmi, I, Averseng, J, Quirant, J, Jamin, F. Deployable double layer tensegrity grid platforms for sea
accessibility. Eng Struct. 2021;231.
[14] Zhang, J, Ohsaki, M, Rimoli, JJ, Kogiso, K. Optimization for energy absorption of 3-dimensional tensegrity
lattice with truncated octahedral units. Compos Struct. 2021;267: 113903.
[15] Rimoli, JJ. A reduced-order model for the dynamic and post-buckling behavior of tensegrity structures. Mech
Mater. 2018;116(SI): 146-157.
[16] Pajunen, K, Johanns, P, Pal, RK, Rimoli, JJ, Daraio, C. Design and impact response of 3d-printable tensegrity-
inspired structures. Mater Design. 2019;182: 107966.
[17] Carpentieri, G, Skelton, RE, Fraternali, F. Minimum mass and optimal complexity of planar tensegrity bridges.
International Journal of Space Structures. 2015;30(3-4): 221-244.
[18] Laccone, F, Malomo, L, Froli, M, Cignoni, P, Pietroni, N, Automatic design of cable‐tensioned glass shells,
Computer Graphics Forum, Wiley Online Library, 2020, pp. 260-273.
[19] Levin, Stephen, M. The tensegrity-truss as a model for spine mechanics: biotensegrity. Journal of Mechanics in
Medicine & Biology. 2002.
[20] Chen, B, Jiang, H. Swimming performance of a tensegrity robotic fish. Soft Robot. 2019;6(4): 520-531.
[21] Chen, M, Goyal, R, Majji, M, Skelton, RE. Design and analysis of a growable artificial gravity space habitat.
Aerosp Sci Technol. 2020;106: 106147.
[22] Wang, R, Parunandi, KS, Yu, D, Kalathil, D, Chakravorty, S. Decoupled data-based approach for learning to
control nonlinear dynamical systems. Ieee T Automat Contr. 2021.
[23] Sun, J, Song, G, Chu, J, Ren, L. An adaptive bioinspired foot mechanism based on tensegrity structures. Soft
Robot. 2019;6(6): 778-789.
[24] Hangai, Y, Wu, M. Analytical method of structural behaviours of a hybrid structure consisting of cables and rigid
structures. Eng Struct. 1999;21(8): 726-736.
[25] Wang, Y, Xu, X, Luo, Y. Topology design of general tensegrity with rigid bodies. Int J Solids Struct. 2020;202:
278-298.
[26] Chen, B, Jiang, H. Instability results from purely rotational stiffness for general tensegrity structure with rigid
bodies. Mech Mach Theory. 2022;167: 104485.
[27] Nagase, K, Skelton, RE. Network and vector forms of tensegrity system dynamics. Mech Res Commun. 2014;59:
14-25.
[28] Kan, Z, Song, N, Peng, H, Chen, B, Song, X. A comprehensive framework for multibody system analysis with
clustered cables: examples of tensegrity structures. Int J Solids Struct. 2021;210-211: 289-309.

29
[29] Kan, Z, Peng, H, Chen, B, Zhong, W. A sliding cable element of multibody dynamics with application to
nonlinear dynamic deployment analysis of clustered tensegrity. Int J Solids Struct. 2018;130-131: 61-79.
[30] Li, F, Peng, H, Yang, H, Kan, Z. A symplectic kinodynamic planning method for cable-driven tensegrity
manipulators in a dynamic environment. Nonlinear Dynam. 2021;106(4): 2919-2941.
[31] Skelton, RE, Oliveira, MCD. Tensegrity systems. Springer, Berlin 2011.
[32] Hurtado, JE. Kinematic and kinetic principles 2016.
[33] Shabana, A. Dynamics of multibody systems. Cambridge university press 2020.
[34] Ma, S, Chen, M, Skelton, RE. Tensegrity system dynamics based on finite element method. Compos Struct.
2022;280: 114838.
[35] Pellegrino, S. Structural computations with the singular value decomposition of the equilibrium matrix.
International Journal of Solids & Structures. 1993;30(21): 3025-3035.
[36] Zhang, L, Li, Y, Cao, Y, Feng, X. Stiffness matrix based form-finding method of tensegrity structures. Eng Struct.
2014;58: 36-48.
[37] Armitage, P, Colton, T. Encyclopedia of biostatistics || optimization and nonlinear equations.
2005;10.1002/0470011815.
[38] Pellegrino, S. Analysis of prestressed mechanisms. Int J Solids Struct. 1990;26(12): 1329-1350.
[39] Calladine, CR, Pellegrino, S. First-order infinitesimal mechanisms. Int J Solids Struct. 1991;27(4): 505-515.
[40] Zhang, JY, Ohsaki, M. Stability conditions for tensegrity structures. International Journal of Solids & Structures.
2007;44(11-12): 3875-3886.
[41] Li, Y, Feng, X, Cao, Y, Gao, H. A monte carlo form-finding method for large scale regular and irregular tensegrity
structures. Int J Solids Struct. 2010;47(14-15): 1888-1898.
[42] Yuan, X, Ma, S, Jiang, S. Form-finding of tensegrity structures based on the levenberg–marquardt method.
Computers and Structures. 2017;192: 171-180.
[43] Zhang, LY, Li, Y, Cao, YP, Feng, XQ, Gao, H. Self-equilibrium and super-stability of truncated regular
polyhedral tensegrity structures: a unified analytical solution. Proceedings of the Royal Society A.
2012;468(2147): 3323-3347.
[44] Zappetti, D, Arandes, R, Ajanic, E, Floreano, D. Variable-stiffness tensegrity spine. Smart Mater Struct.
2020;29(7): 75013.
[45] Sabelhaus, AP, Zhao, H, Zhu, EL, Agogino, AK, Agogino, AM. Model-predictive control with inverse statics
optimization for tensegrity spine robots. Ieee T Contr Syst T. 2020;29(1): 263-277.

30