Dynamics of Tensegrity Robots With Negative Stiffness Elements

Tensegrity structures have unique features such as low mass to payload ratio, strength, and robustness. Therefore, they present great potential in robotics, aerospace, and civil engineering. The dynamics of tensegrity robots is highly nonlinear and constrained. As a result, their modeling, simulation, state estimation, and control are non-trivial. Strings in tensegrity structures are usually modeled as linear springs. Utilization of nonlinear elastic/damping elements in tensegrities would further enrich their dynamics and endow them with additional properties, such as multiple equilibrium configurations. In this paper, our preliminary work on the dynamics of actuated tensegrities with strings containing nonlinear elastic and/or damping elements is presented. At first, the formulation of tensegrity dynamics with general nonlinear elastic/damping elements is explored. Later dynamics of tensegrities with negative stiffness honeycombs incorporated into strings are considered. Simulations are performed on three tensegrity systems: two-bar, three-bar, and six-bar structures. Results demonstrate that negative stiffness honeycombs result in nonlinear steady-state response to constant external force, reduced force magnitudes in strings and bars, and increased range of motion.


I. INTRODUCTION
The term tensegrity was coined by the famous architect Buckminster Fuller as an elision of tension and integrity [1]. Tensegrities are lightweight structures consisting of bars and strings. According to the convention, the term ''bar'' denotes an element which can be either under compression or tension, while ''string'' denotes an element which can be only under tension [1]. The advantages of these structures include scalability, high payload-to-mass ratio, and robustness. Therefore, they have significant potential as robot manipulators. However, this requires further work on a number of critical areas such as modeling, system identification, control, and state estimation of tensegrity systems. The solution of these problems would open the door for widespread utilization of tensegrities in robotics.
The researchers initially investigated the equilibrium and stability conditions of tensegrities [2]- [5]. For example, The associate editor coordinating the review of this manuscript and approving it for publication was Yingxiang Liu . necessary stability condition for tensegrity systems based on positive definiteness of the tangent stiffness matrix was presented in [2]. Quadratic energy functions were utilized to investigate the rigidity of convex polygons composed of bars and strings [3]. Here ''rigidity'' means that a tensegrity structure does not collapse as it undergoes an arbitrary continuous motion or deformation, i.e. tensegrity behaves as a rigid body. Rigidity of prismatic symmetric tensegrities was considered in [4]. Form-finding of tensegrity structures refers to the method of finding their equilibrium configuration [5]. Different methods and algorithms, such as genetic algorithm [6], [7], iterative procedure for statically indeterminate structures [8], Monte Carlo method [9], finite element methods [10] and others, were used to solve a form-finding problem of tensegrities.
Later the emphasis shifted toward understanding and accurate modeling of their dynamics [11]- [16]. Specifically, Newton's second law was used to derive equations of motion of tensegrities in Eulerian (spatial) and Lagrangian (material) formulations [11]. Here Eulerian formulation means that nodal position and velocity vectors with respect to global coordinate system were used in dynamic equations, while Lagrangian formulation means that relative nodal displacement and velocity vectors from initial configuration were utilized. Matrix differential equations for tensegrity dynamics in non-minimal form were developed in [13], [14]. Later this theory was generalized for class-k tensegrities [17], for matrix differential equations in minimal form [17], for stacked tensegrities including damping terms [16], and for tensegrities with non-negligible string masses [18]. These advances paved the way for scientists and engineers to further expand research related to tensegrity structures. For instance, problems such as combining tensegrities with nonlinear structures, utilizing tensegrities in vibration isolation and shock absorption applications, designing tensegrities with new topology and including metamaterials are interesting areas for research.
Conventional systems exhibit a positive correlation between deformation and an internal resistance force, which is equal to the externally applied quasi-static load. Rarely, systems possess a negative stiffness behavior, i.e. increased deformation results in a reduction of the internal reaction force. For example, these kinds of structures were proven effective in acoustic and vibration isolation [19]- [22]. One such system is the buckled beam. If a load is applied perpendicular to the buckled beam, it first resists deformation. As the external load exceeds a certain threshold, sudden buckling occurs, and the beam quickly transitions into an alternative buckled state. If negative stiffness beams are grouped in series and/or in parallel, then negative stiffness honeycombs (NSHs) are formed. NSHs are effective in applications requiring energy absorption, shock isolation, and controllable switching [23]- [26]. Additionally, the deformation of NSHs is fully recoverable, which means that NSHs can be utilized to design effective metamaterials intended for multiple cycles of energy absorption [27], [28].
Tensegrities may also have negative stiffness behavior under certain conditions [29]. This behavior is dependent on the structural parameters of the system and exhibits itself under torsion. In [29], the researchers simulated three groups of prismatic tensegrities with different structural parameters. Each structure was subjected to a rotation from its upper part around its normal axis. The torsional stiffness of the structure was negative when the rotation angle was between the critical and maximum angular position. Tensegrities possessing quasi-zero stiffness property were also reported in the literature [30]. Quasi-zero stiffness implies that a tensegrity has a continuous range of configurations which are at equilibrium. This property is used to design different sensors, such as vibration sensors [31].
Extreme stiffening and softening of tensegrity prisms under axial loads were investigated in [32]. The behavior of two tensegrity prisms, 'thick' and 'slender' variants, was analyzed both numerically and experimentally. The two models differed from each other in terms of dimensional characteristics. The results indicated that the 'thick' prism model possesses softening response in compression at lower pre-strain and a stiffening response in tension at greater prestrain. In contrast, 'slender' tensegrity structures exhibit the opposite response. In essence, softening response is a snap buckling behavior that leads to the complete axial collapse of the system. Bistable behavior of tensegrity structures with different potential energies was studied in [33]. It was shown that the bistable response depends on system geometry and pre-strain values of the strings. The study concluded that the transformation of the system from a higher energy state to a lower one under an elastic deformation causes instability. In this process, the value of the externally applied load would first increase from zero to a positive value and then suddenly decrease, which implies negative stiffness behavior.
Tensegrity structures have strings, which are modeled as linear springs (Fig. 1a). Replacing strings with nonlinear stiffness elements will alter the dynamics of tensegrities. For example, they might obtain quasi-equilibrium states, increased safety, and impact absorption capabilities. Additionally, linear and rotational stiffness of tensegrities might be controlled with the presence of nonlinear stiffness elements.
Recently, researchers focused their attention on tensegrity robots, which are tensegrity structures with incorporated actuation mechanisms [34]- [38]. Compared to traditional robot manipulators, tensegrity robots possess several advantages, such as low weight, relatively large strength, selfdeformation, and folding abilities, easy deployability, and capability to safely survive brief shocks and impacts [35], [36]. Incorporating NSHs into tensegrity robots might further improve their static and dynamic characteristics. For instance, tensegrity robots with NSHs might possess more pronounced nonlinear behavior, enriched variable stiffness behavior, enhanced damping and energy absorption capabilities, and multiple stable, quasi-stable and zero-stiffness configurations. Traditionally, system nonlinearities are viewed as disadvantages. However, nonlinearities present in the system might in some cases be beneficial, such as in vibration isolation with lower energy consumption [39]. Currently, nonlinear behavior of tensegrity structures is due to their topology (different ways of connecting bars and string with each other) and geometry (location and orientation of bars and strings). Our main motivation is that by including nonlinear stiffness/damping elements, dynamics of tensegrity structures can be further enriched and varied. Only tensegrities with linear string behavior were considered before. Motivated by the recent progress in NSH modeling [40], the main purpose of this paper is to investigate the effects of incorporating NSH into tensegrity robots. Specifically, our scientific contributions are twofold: Firstly, the dynamics of tensegrities with nonlinear stiffness elements is formulated. Secondly, a general methodology to model the dynamics of tensegrity robot manipulators is developed.
In this work, we investigate the dynamics of tensegrity robots with nonlinear elastic elements. After covering the general case of nonlinear springs, leveraging the progress in the analytical modeling of NSHs [40], we focus on tensegrities where strings are replaced by NSHs (Fig. 1b-c). To the best of our knowledge, this is the first work incorporating nonlinear stiffness elements into the dynamics formulation of tensegrities. The rest of the paper is organized as follows. Section II reviews the general dynamics of tensegrity structures, discusses the modeling of tensegrities with nonlinear elastic elements, and finally derives equations for tensegrities with NSHs. This is followed by the description of the case study in Section III, which also provides simulation scenarios. Extensive simulation results are presented in Section IV. The paper is concluded in Section V.

A. DYNAMICS OF CLASSICAL TENSEGRITY STRUCTURES
Tensegrity structures consist of long rigid bars interconnected with each other through strings (Fig. 1a). Dynamics of tensegrity structures for different configurations were presented in [16]- [18]. In this subsection, we will briefly summarize the important aspects of tensegrity dynamics formulation to demonstrate the connection between our contribution and the existing theory. Let us assume that the general tensegrity structure consists of µ number of physical nodes, β number of rigid bars, and α number of strings. A node is a point where a string and a bar are connected to each other, and is denoted by a position vector n i ∈ R 3 , i ∈ 1, . . . , µ. Node position vectors are assembled into a node matrix N ∈ R 3×2β = [n 1 , . . . , n 2β ], which has more nodes (columns) than the number of physical nodes. These additional nodes are called ''virtual nodes'' and they are required to formulate the full system dynamics of tensegrity structures. In essence, we assume that the bars are connected to each other only through strings, which results in the so-called class-1 tensegrity structure (More information can be found in [16]- [18]). By defining connectivity matrices of bars C b ∈ R β×2β , strings C s ∈ R α×2β , and center of mass of bars C r ∈ R β×2β , we can find the corresponding bar B = NC T b ∈ R 3×β , string S = NC T s ∈ R 3×α and bar center of mass R = NC T r ∈ R 3×β matrices where T denotes the matrix transpose operation. The bar matrix contains the vector of all bars where each bar vector can be represented as a difference of two node vectors b i = n j − n k ∈ R 3 describing two ends of the bar, where j, k ∈ 1, . . . , 2β and j = k. Similarly, the bar center of mass matrix contains the bar center of mass vectors as R = [r 1 , . . . , r β ], where r i = (n j + n k )/2 ∈ R 3 . The string matrix consists of the string vectors S = [s 1 , . . . , s α ], where a single string vector is s i = n j − n k ∈ R 3 , j, k ∈ 1, . . . , 2β and j = k. Tensegrity dynamics in matrix form can be written asN where M = C T bĴ C b + C T rm C r ∈ R 2β×2β is the matrix which includes mass and inertia of all bars. Operator |ˆ| denotes diagonalization: given an arbitrary vector x ∈ R y , y ∈ Z, this operator outputs a diagonal matrixx ∈ R y×y , and viceversa, given a square matrix, it returns a vector of diagonal elements. J = [J 1 , . . . , J β ] ∈ R β is the vector of moments of inertia of bars. J i is the moment of inertia of each bar at the center of mass and perpendicular to its long axis. For example, if a bar is taken as a long solid bar with mass m i and length l i = ||b i ||, where operator ||.|| denotes Euclidean norm, then J i = m i l 2 i /12. Similarly, m = [m 1 , . . . , m β ] ∈ R β is the vector of bar masses. The matrix K contains string and bar forces as is the vector of internal force density magnitudes in strings. While η = [η 1 , . . . , η α ] ∈ R α is a similar vector of frictional force density magnitudes in strings. The matrix λ ∈ R β×β denotes the force densities in bars as λ = −l −2 (ĴḂ TḂ represents the force vectors acting on each node, including virtual nodes. Please note that the λ term in K contains double hat operator, which means that diagonalization is performed twice. It is equivalent to setting all off-diagonal elements in matrix λ equal to zero. The matrix of external force vectors acting on each node, which, for instance, might include payload weight, is denoted as F e ∈ R 3×2β . While F c ∈ R 3×2β is the matrix of constraint force vectors acting on each node. Constraint forces arise due to constraints on node positions. For example, if bottom nodes are fixed to the ground, then F c contains the corresponding constraint forces which do not allow the bottom nodes to move. More information about constraint forces and their formulation can be found in [17]. T s = Sγ ∈ R 3×α is the matrix of force vectors representing internal forces in the strings. Similarly, D s = Sη ∈ R 3×α is the matrix of damping force vectors within the strings.
For regular strings, it is assumed that where ξ i is the length of the string i (i = 1, . . . , α), ξ i = ξ i − ξ 0 i with ξ 0 i denoting constant initial (unstretched) length of the string, while its linear stiffness coefficient is denoted as k i . If strings have a constant cross-section A and are made out of a homogeneous, isotropic and linear material with Young's modulus E, then k i = EA/ξ 0 i . String length can be found as ξ i = ||s i ||. The damping force densities can be estimated as where d i is the damping coefficient (i = 1, . . . , α) and (4) and the state transition matrix where 0 ∈ R 2β×2β and I ∈ R 2β×2β are, respectively, zero and identity matrices, and a new matrix Then, using (4), (5), and (6), the full tensegrity dynamics can be written as a first-order matrix differential equatioṅ Equation (7) can be integrated to obtain the node coordinates as a function of time.

B. MODELING OF TENSEGRITY ROBOTS
Due to highly nonlinear dynamics, high-dimensional state space representation, and multiple constraints, the control problem of tensegrity robots is a challenging task [41]. Different control methods, such as model-based state feedback, linear quadratic regulator, model-free data-based control, and neural networks were proposed to tackle the control problem [38], [41]. The main physical approach to manipulate tensegrity robots is through controlling the string lengths. In other words, by pulling the strings in or out, we can regulate the configuration, size, and shape of tensegrity robots.
Dynamics of tensegrity robots is similar to the dynamics of tensegrity structures described in subsection II-A. However, the main difference comes from time-varying initial string lengths ξ 0 i (t), which effectively model the tensioning and releasing of the strings using actuators. Additionally, dynamics of force/torque-controlled actuators needs to be taken into account as where J a and d a are, respectively, inertia and damping properties of actuators, k a is the torque constant, i a is the input current, while f a is an external force/torque (in our case f a = γ i ξ i ). If a reduction gear is utilized in the actuator, then the aforementioned actuator parameters should already contain the overall gear ratio. θ denotes the angular position of the rotary or linear position of linear actuators. By incorporating (8) into (7), full state-space dynamics of tensegrity robots with actuator dynamics can be obtained. In this case, state coordinates are formed from node coordinates and velocities, as well as from actuator positions and velocities, while input consists of a set of input currents for each actuator. The obtained state space representation is general, valid for an arbitrary tensegrity robot or structure.

C. TENSEGRITY STRUCTURES WITH NONLINEAR STIFFNESS ELEMENTS
Strings in tensegrities are modeled as linear springs. In order to obtain variable stiffness behavior, nonlinear elastic elements between the nodes are required. For example, if negative stiffness elements are embedded in the strings, then tensegrity would have quasi-steady configurations. Quasi-steady configuration is an equilibrium configuration of tensegrity, which is unstable or neutral. In other words, in steady (stable) configuration, slight disturbances will not change the configuration of tensegrity. However, in quasisteady configuration, small disturbances will often cause transformation to other more stable configuration. By carefully designing these configurations coincident with robot positions in tasks, such as pick and place, energy consumption and material wear can be reduced. Therefore, in this subsection, we describe the general problem of modeling tensegrity structures with nonlinear stiffness elements. The modeling of tensegrity structures with nonlinear elastic/damping elements is not trivial. The numerical integration of (7) follows a three-step sequence. Firstly, for the state matrix X (node coordinates and velocities) at current time, elastic and damping forces of the strings can be computed from (2) and (3). Secondly, these forces can be used to findẊ from (7). Thirdly,Ẋ is integrated to find the new state matrix for the next sampling time. If, instead, nonlinear elastic/damping elements are used, then it might be difficult or impossible to find forces from nodal positions and/or velocities (the first step). Still, there is a large group of nonlinear elastic/damping elements, for which tensegrity dynamics formulation (7) can be used. For example, if nonlinear elastic element is described by explicit nonlinear force-displacement characteristic f = f s (ξ ), where f denotes tensile force, then it results in nonlinear string force density In this case, the tensegrity dynamics formulation in section II-A can still be utilized. The only modification is replacing (2) with nonlinear function (9). Similar reasoning applies for nonlinear damping elements that can be described by explicit nonlinear force-velocity characteristics f = f d (ξ ), which results in the following string force density However, in the general case, there is no explicit relation between string length (node position) and the corresponding tensile force, or between string stretching (node velocity) and the respective damping force. For example, if nonlinear stiffness elements possess hysteresis behavior, then forces can not be expressed as simple functions of node positions. For hysteresis, the time evolution of the elastic element deformation needs to be taken into account. Alternatively, if nonlinear elastic/damping elements are stacked in series, then it becomes impossible to find the tensile force from string displacement. This is due to the fact that for all elastic/damping elements stacked in series, the force will be the same. While the string elongation ξ i will be split unevenly between each nonlinear element in the stack. As a result, we have a dilemma: force cannot be calculated given the current node positions, while force is required to find individual compression of each nonlinear elastic element in the stack. The sum of individual compression of each element would give the total displacement of the end nodes. This problem does not arise if nonlinear elastic/damping elements are stacked in parallel. In this case, the displacement will be the same for all elements, while the total force will be a simple sum of the forces of each element. Also, this problem does not arise if different linear springs are connected in series.

D. NEGATIVE STIFFNESS BEAM
In the special case of negative stiffness elements stacked in series, it is possible to find the general expression for force-displacement characteristics. However, before we model NSH in series, let us briefly describe the model of a single negative stiffness element. There are multiple approaches to how negative stiffness behavior can be achieved, and using buckled beams is one of them. In this subsection, the force-displacement characteristics of a single negative stiffness buckled beam is considered (Fig. 2). Negative stiffness behavior can be bistable or metastable. Bistable behavior means that there are two stable configurations, and changing from one to another requires external force. Metastable behavior means a beam buckles under external force, but as the force is removed the beam returns to its initial configuration. More details about this can be found in [42].
A single beam with initial curve according to w 0 (x ) = h 2 (1 − cos (2π x p )), where h is the initial height of the undeformed beam at the center location (x = 0.5p), p is the horizontal length of the beam, has the following vertical force-displacement characteristics where Q is the constant defined by beam material and geometric properties, ζ is dimensionless design parameter, with the value in the range 2 < ζ < ∞, u = w 0 (0.5p) − w(0.5p) is the vertical displacement of the center point of the beam from its initial position. Parameter Q = Eat 3 hπ 4 12p 3 depends on beam's Young's modulus E, length p, thickness t, width a, and h [40]. w(x) denotes the curve of the beam after deformation. If the design parameter ζ is chosen such that ζ < 8, then the resultant negative stiffness beam will have meta-stable behavior. While ζ > 8 leads to bistable negative stiffness behavior. Since strings in a tensegrity structure should act as springs, negative stiffness beams should be metastable (2 ≤ ζ ≤ 8). Due to the symmetry of the initial shape of the beam, 0 ≤ u h ≤ 2 in (11).

E. NEGATIVE STIFFNESS HONEYCOMBS
Now the question arises, how do we model overall force-displacement characteristics of the metastable NSH composed from multiple different negative stiffness beams? NSH is obtained by stacking negative stiffness beams in series or parallel. We will consider the former system because the latter system has an additional rotational degree of freedom, which is not needed for our task. One can notice that force-displacement relationship for a negative stiffness beam (11) is a cubic function, as shown in Fig. 3 plotted for Q = 6 N, ζ = 5, h = 1.5 cm. As was mentioned above, force is constant for nonlinear springs connected in series. Therefore, the general strategy to find the overall force-displacement characteristics of NSH is the following. For a given force, we find the displacement of each beam, and then add these displacements to find the net displacement at the given force value. This process can be repeated for as many force values as desirable. In the end, we obtain the overall displacement-force characteristic of the NSH. If this characteristic is inverted, then we obtain the overall force-displacement characteristic. However, there are some caveats for NSHs. When one of the beams buckles, then the whole NSH rapidly undergoes displacement collapse equivalent to the value of the corresponding collapse of the single beam. The cubic function (11) has three roots, which can be found as where Given the force value f , the roots u I , u II and u III give the corresponding solutions of displacements in regions I , II and III (Fig. 3). These regions cover the following displacement ranges: 0 < u h < u max h for region I , u max h < u h < u min h for region II , u min h < u h for region III . Real parts of the solutions needs to be considered in (12). From (11), the following extrema and important points can be found (Fig. 3) As a result, the Algorithm 1 can be devised to obtain the net force-displacement behavior for the NSH. Please note that the force-displacement function will have discontinuities after the loop in the Algorithm 1. In order to obtain a smooth function, we have to insert into the net force-displacement function the part of the individual curve below the blue line in Fig. 3 for each beam. Finally, the algorithm outputs the net force-displacement function for NSH with beams connected in series. This function will be nonlinear and will be uniquely defined for each net displacement of NSH. NSH with desired buckling force thresholds and buckling displacements can be designed by varying parameters of beams, such as type of material, thickness, length, and width [40]. The net force-displacement characteristics of NSH is obtained in a table form (as a set of coupled numbers: (force, displacement)), and then the table can be used as a sub-function f ( ξ ) in tensegrity dynamics. In other words, this sub-function would output f i through interpolation given the net NSH displacement ξ i . The displacement ξ i in turn can be found from the current node positions N . By utilizing the obtained net force-displacement characteristics of NSH in (9) and (10), dynamics of tensegrity structures with NSHs instead of springs is obtained (Fig. 1c).  (14); (Fig. 3); end; end for After the loop f will have discontinuities at f = f max values. In this case, we have to define f so that f min < f < f max and for each f apply ξ = ξ + u II for u < u min and ξ = ξ + u III for u > u max in order to obtain the smooth function; return f and ξ

III. SIMULATIONS
Three tensegrity systems will be considered for our simulation case study: 1) 2D tensegrity with two bars, three strings and four nodes, 2) 3D tensegrity with three bars, six strings, and six nodes, and 3) 3D tensegrity with six bars, 21 strings and twelve nodes (Fig. 4). 2D tensegrity was assumed to be in the vertical plane, so that gravity has to be considered, and each of its two vertical strings contains NSHs. Three bar tensegrity also has three NSH along its three vertical strings, VOLUME 8, 2020 while six-bar tensegrity has six NSHs along its upper six vertical strings. If NSHs are connected to strings, then they have to behave like springs, and as a result, NSH should consist of meta-stable beams. Also, NSH operates under compression, so for them to function under tension, string connection to NSH can be simply swapped. For example, in normal case, a string which is on the right side of a simple extension spring is connected to the right end of spring, and vice versa for the left string. However, in our case, we assume that a string that is on the right side of NSH is connected to its left end and similarly for the left string. This way NSH would function under tension. For each tensegrity system, the bottom nodes (n 1 and n 2 for the 2D system and n 1 , n 2 , and n 3 for the 3D systems) are rigidly fixed to the ground. In Fig. 4, a node n i is denoted with a red number i ∈ N. All NSHs are assumed to be identical. NSH consists of four meta-stable beams in series and was designed with the following parameters: Q = [6 8 10 12] N, ζ = [3 4 5 6], and h = [1 1 1.5 1] cm. Utilizing the procedure outlined in the Algorithm 1, the net force-displacement characteristics of the NSH is obtained (Fig. 5). We would like to remind that the Algorithm outputs the net force of NSH as a function of net displacement (displacement is input, while force is output). This is due to the requirement that at a given time step (node position), we need to know the forces acting on the nodes in order to perform one step ahead integration. However, in order to generate the net force-displacement characteristics of NSH, inverse procedure is used: given force, displacements of each negative stiffness beam is found, and then all displacements are added. The maximum compression of NSH is slightly more than 8 cm. We assumed the following parameters for the three systems (for simplicity, all elements are taken as identical): bar masses m i = 1 kg, string stiffnesses k i = 15 · 10 3 N/m for 2D system and k i = 5 · 10 3 N/m for other two systems, string damping coefficients d i = 500 Ns/m. In reality, due to high value of tension forces and unaccounted Coulomb frictional forces, damping coefficient has even larger value [15], [16]. The tensegrity systems' initial node coordinates are given in Table 1. String and bar connectivity matrices are provided in Appendix. Bar lengths for these systems are 1.29 m and 2.07 m for the 2D and 3D tensegrities, respectively. Initial forces in strings and bars, which would satisfy the equilibrium condition of tensegrity systems corresponding to initial node positions, were calculated beforehand.
The simulations were performed in MATLAB with ode45 medium order solver with both relative and absolute tolerance criteria set to 10 −4 , while gross simulation and integrator time steps were equal to 10 −2 and 10 −3 s, respectively. Multiple simulations were performed: I) One set of simulations was performed with gross simulation time equal to 50 s. A constant input force was applied to the top nodes (n 3 , n 4 for the 2D system, n 4 -n 6 for the three-bar system and n 10 -n 12 for the six-bar system) in vertical direction in all three tensegrity configurations. Then, the system was allowed to come to a new steady-state equilibrium configuration under the influence of the external forces.  Simulations were repeated for different values of the input force. The range of input force values (from 10 N to 72 N with step size of 2 N) was chosen to result in snaps of all four arcs of the NSH element. For each value of constant input force, a steady state height of top nodes was recorded. For comparison, identical simulations of simple tensegrity systems without NSHs were also performed. II) Additional simulations on the three-bar tensegrity were performed to investigate in more detail the effect of large external forces on the node displacements and string forces. In these simulations, external periodic input force with magnitudes -400 N and 400 N, period of 1 s with duty cycle of 50% were applied to the node n 4 in the vertical direction. The node positions and elastic string forces were recorded for two tensegrity systems: with and without NSH. In order to test the effect of the damping term on NSH tensegrity dynamics, extra simulation runs were performed. In this case, damping terms were set to d i = 50 Ns/m and d i = 5000 Ns/m values, while all other parameters were not changed.
III) In order to show the actuation principle of the tensegrity structures, additional simulations were performed on each of the three tensegrity systems. Actuation was imposed by varying in a sinusoidal manner the unstretched string length ξ 0 i of the top horizontal strings (a single string connecting nodes n 3 − n 4 for the two-bar system, three strings connecting nodes n 4 − n 5 − n 6 for the three-bar system and similar three strings connecting nodes n 10 − n 11 − n 12 for the six-bar tensegrity). For the three bar and six bar tensegrities, actuation was symmetric, i.e. three strings were actuated by the same function without phase differences. In essence, actuators would pull and release the string. Therefore, for simplicity, we ignored actuator dynamics (8) and modeled input to the system taking place on the string side. Unstretched string length was varied sinusoidally with the amplitude of 0.15 m and period of 2.5 s. The simulation duration was set to 10 s, systems were undisturbed for the initial 2 s, then from 2 s until the end of simulation input to the system was applied.

IV. RESULTS AND DISCUSSION
I) Results of the aforementioned steady-state simulations are shown in Fig. 6. Specifically, steady-state vertical coordinatez of the upper node positions versus external input force f ext is depicted for two-bar, three-bar and six-bar tensegrity structures. It can be observed from Fig. 3 or 5 that if displacement is imposed on the NSH structure, then the smooth trajectory with negative stiffness region is obtained for the force-displacement curve. However, if instead force is imposed on the NSH structure, then horizontal ''displacement jumps'' should be observed at force values equal to f max of each beam. In other words, those snaps occur when the external force is just above the extremum point f max . This phenomenon is also visible from the simulation results, black solid lines in Fig. 6. There are four snaps in the position-force plot of the top nodes and this is due to the four beams present in NSH. Since the same amount of vertical force has been applied to all upper nodes (two nodes for 2D system and three nodes for 3D systems), vertical strings undergo the same amount of tension, therefore, all the top nodes have the same vertical z position and snap at the same time.
Simulation results of simple tensegrity systems without NSH are also shown as black dashed lines in Fig. 6. For these systems, the steady-state positionz as a function of external force f ext shows smooth and almost linear behavior. Strings are modeled as simple linear springs, and as a result, the observed behavior is expected. Additionally, the following qualitative explanation can be given to the obtained results. The following dynamic equation can be written for the lumped parameter model of an arbitrary mechanical system:  where M, B, K are, respectively, mass, damping and stiffness parameters, q is coordinate and f ext is an external force. For the steady-state condition under constant external force, all time differentials disappear, and so we get steady-state positionq as a linear function of external force. As a result, reachable workspace for simple tensegrity systems might be a single continuous domain. However, for nonlinear systems, the dependence ofq on f ext would be nonlinear. Therefore, for tensegrity systems with NSH, the reachable workspace would consist of multiple continuous domains not intersecting each other. It can also be observed that the range of variation ofz is larger for NSH-based tensegrity systems compared to simple tensegrity systems. This would imply that for the same force input to the system, NSH-based tensegrity is more compliant. Therefore, NSH-based tensegrities might have better shock absorption and impact mitigation capabilities compared to traditional tensegrity systems. II) Simulation results for the three-bar tensegrity under cyclic external load are shown in Fig. 7-8. The difference from their corresponding initial values of x, y, z-coordinates of the node n 5 for tensegrity with and without NSH are shown in the Fig. 7a-b. It can be observed that tensegrity without NSH has higher values of node displacements compared to tensegrity with NSH. This result is counter-intuitive because NSH ought to increase system compliance, thus increasing the displacement of tensegrity nodes under the influence of the external forces. The reason for this behavior will become clear after we examine simulations results at other damping values. Simulation results for the three-bar tensegrity with damping term set to d i = 50 Ns/m are shown in the Fig. 7c-d,  FIGURE 9. Results of steady simulations with actuation of strings.
while simulation results at d i = 5000 Ns/m are shown in the Fig. 7e-f. It can be observed that compared to the initial case ( Fig. 7a-b), node displacements in the lower damping term case are larger. This means that decreasing the damping term increases the node displacements as expected. Also, node displacements in the tensegrity without NSH are smaller than in the tensegrity with NSH. In other words, lowering the damping term gives the anticipated result, in which tensegrity with NSH is more compliant than tensegrity without NSH. Increasing the damping term leads to much smaller node displacements compared to the initial case. Also, in the higher damping term case, tensegrity without NSH is more compliant than tensegrity with NSH, which is similar to the results in Fig. 7a-b. In other words, increasing the damping term reduces overall node displacements and decreases the compliance of tensegrities with NSH. The reason for the observed behavior might be also due to the buckling of specific beams and geometric factors (location and direction of the external forces with respect to the node n 5 ). Moreover, node displacements clearly demonstrate the stabilizing effect of damping forces. Under the effect of external periodic loading, tensegrity with high damping terms quickly achieves periodic forced response (all transient effects quickly die out within the first period). While in low damping term case, node displacements have transient effect even at the fourth period. At high damping term values, NSH also displays stabilizing effect: presence of NSH is equivalent to increase of damping term. For example, at d i = 500 Ns/m the steady-state response of the three-bar tensegrity with NSH under cyclic loading is more than two times smaller than the response of tensegrity without NSH. Elastic string force along the string n 1 -n 4 is shown in Fig. 8 for two systems: with and without NSH (d i = 500 Ns/m). String forces in tensegrities with NSH have much smaller magnitudes, and also the variation of string forces due to NSH buckling can be observed. It implies that NSH reduces internal elastic forces in strings and bars and helps to mitigate the possible damaging effect of external forces. These findings again reinforce our conclusion that incorporating NSH into tensegrities enriches their dynamics and improve their capability to withstand external impacts. III) Actuated simulation results are presented in Fig. 9, where relative displacement coordinates ( x = x(t) − x(0) and similarly for y and z) of the last nodes of each of the systems are shown. For clarity, simulation results for the time interval 2-10 s are depicted. It can be observed that in the two-bar system, node coordinate changes are larger for the system with NSH compared to the system without NSH. A similar conclusion cannot be made for the three-bar and six-bar systems. The observed pattern is cyclic, but nonlinear, which shows that tensegrity dynamics is nonlinear. In the six-bar results, oscillations can be observed. This indicates that damping is not sufficient to fully suppress the transient oscillations. Relative displacement coordinates are the largest for the three-bar system. The difference between the behavior of tensegrities with and without NSH are the most pronounced for the two-bar tensegrity. Overall, we can conclude that NSH alters tensegrity dynamics significantly.
IV) In order to better understand the dampening effect of NSH, additional simulations were performed on a 1D simple lumped-parameter mass-spring-damper system with NSH and without NSH, as given in (15). In this case, f ext was a step function with amplitude 100 N, period of 1 s and duty cycle of 50%. Results of external power, which is calculated asżf ext , for this system are shown in Fig. 10. The tensegrity with NSH has more overall absorbed energy     (the area under the power-time curve) than the tensegrity without NSH. A video visualizing these simulations is provided as supplementary material.

V. CONCLUSION
In this work, preliminary results for describing tensegrity dynamics with nonlinear elastic elements are presented. Specifically, the effects of negative stiffness honeycombs on tensegrities are investigated. Simulations are performed for three tensegrity systems: 2D two-bar, 3D three-bar, and 3D six-bar structures. For each system, two types of simulations were performed: one with negative stiffness honeycombs and one without them. Results demonstrate that negative stiffness honeycombs result in interesting steady-state behavior, reduced string forces, increased energy absorption capabilities, and mitigation of the destructive effects of external forces.
Future works include experimental verification of the tensegrity dynamics formulation and consideration of NSH masses in tensegrity dynamics.

APPENDIX
Here, string and bar connectivity matrices are given for the considered three systems. In the tables, 1 denotes beginning node of a string (bar), while −1 denotes their end-node.