Kinematics Analysis and Optimization of a 3-DOF Planar Tensegrity Manipulator under Workspace Constraint

Tensegrity mechanisms (TM) are well-appraised for their compliance and lightweight, making their design theory a hot research topic recently. However, due to unconstrained DOFs, the design and analysis of TMs are generally more complicated than traditional mechanisms composed of rigid links and joints. A compact 3-DOF tensegrity manipulator is introduced and an advanced two-step kinematic modeling method is proposed in this paper. This method is first assumed that bars and cables are rigid to estimate the equilibrium state using the energy-based method. Then, the flexibility of bars and cables is considered, and the force density method equations are solved utilizing the equilibrium state obtained by the previous step as the initial guess of iteration for fast computation. Based on the two-step method, the performances of the TM, such as workspace, manipulability, potential energy, and stiffness, are analyzed. Thereafter, the installation position and size of the manipulator are optimized under the workspace constraint. In the optimization process, discrete points on the prescribed task workspace contour are used to simplify the constraint to improve computational efficiency. Finally, study cases are investigated to validate the proposed method, and the feasibility of the discrete sampling method for constraint simplification is also verified.


Introduction
The word Tensegrity (tensile-integrity) was first proposed by Fuller [1] to describe a structure only containing compression members (struts) and tension members (cables or springs). As a structure characterized by self-stressing and high strength-to-weight ratio, tensegrity quickly gained attention in the fields of architecture [2,3], biology [4][5][6], and art. In the past 20 years, there has been an increasing interest in its applications in deployable mechanisms [7][8][9], manipulators [10], and mobile robots [11][12][13][14]. Generally, tensegrity mechanisms (TM) use springs to achieve a larger workspace on the premise that the mechanism is tensioned and stable. For some applications, TM is advantageous due to its compliance. However, unlike traditional mechanisms, TM has unconstrained DOFs and its pose is determined by solving self-equilibrium equations. Therefore, except for 1-DOF mechanisms [15], it is formidable to obtain an analytical kinematic model of TM, which brings challenges to its design and analysis.
The kinematic modeling method for TMs is being studied persistently and various methods have been developed. The energy-based modeling method determines the equilibrium position through the minimum point of potential energy and can be used for kinematic modeling of TMs. Following this track, Arsenault et al. [16] gave a kinematics model of a 2-DOF tensegrity mechanism using an energy-based approach and analyzed its stiffness and dynamic characteristics. They pointed out that the actual control of TM is more involved than the traditional mechanism due to unconstrained DOFs. Boehler et al. [17,18] used an energy-based approach to calculate the workspace of a planar TM. They further studied the effects of the springs on stiffness, workspace, and actuation requirements. Wenger et al. [19] analyzed the multi-solution problem of a planar TM with three linear springs and classified the solutions according to stability. However, in the energy method, the second partial derivative of the potential energy function is required to determine the stable solution, which is complicated for complex mechanisms. Apart from the energy approach, the principle of virtual work [10], screw method [20], and force density method (FDM) [21] are also used for kinematic modeling.
In addition to the kinematic modeling technique, configuration and actuation design of planar TMs also interests researchers. One typical planar TM structure is obtained by adding actuators to the Snelson Cross (SC), an X-shape planar tensegrity formed by two crossing bars and four cables. Moored et al. [22] designed and analyzed an artificial pectoral fin based on SC, and the multiple cable-routed and strut-routed actuation strategies were adopted. Begey et al. [23] took actuated SC as an example to study the workspace, manipulability, potential energy, and internal forces under different actuation modes. Their research revealed that the actuation mode has significant effects on the behavior of the Snelson Cross. Based on the idea of using SCs as building blocks for assembly, they investigated two planar tensegrity manipulators that meet a remote center of motion constraint [24]. Generally, the performance of the TM is directly determined by installation position, size, and springs. However, it is arduous to obtain the optimal design due to its self-stressing characteristics. In practical applications, a proper workspace of the TM is required, and thus the problem of optimal design for TM under workspace constraint should be investigated.
A 3-DOF planar tensegrity manipulator, considering the flexibility of cables and bars, is introduced in this work. Due to the flexibility of members, the analytical expression of potential energy and its partial derivative are complex, which brings challenges to the form-finding of the energy-based approach. Meanwhile, the stiffness of the bars is generally much greater than that of the springs, making the equilibrium matrix close to singularity when the FDM is used for modeling. The initial guess has a significant effect on the efficiency of the iterative approach in solving the FDM equilibrium equation. In view of the above problems, a two-step kinematic modeling method combining the energy-based approach and FDM is proposed to analyze the workspace, manipulability, and stiffness of the TM. This method first assumes that the bars and cables are rigid to estimate the self-equilibrium state using the energy-based approach. Then, considering the flexibility, the FDM equilibrium equation is established, and the state obtained in the previous step is used as the initial guess of iteration. In this way, determining the initial value of FDM iteration is addressed, and the efficiency of kinematics computation is improved. Besides, a minimum size optimization method under arbitrary workspace constraints is proposed. In order to improve computational efficiency, discrete points on the contour of a prescribed target workspace are used to simplify the constraint in the proposed method.
The rest of this paper is organized as follows. Section 2 presents the 3-DOF planar TM and its configuration and actuation mode. In Section 3, the two-step kinematic model is derived, and the solving procedure is discussed. Section 4 introduces definitions of performance criteria, such as workspace, manipulability, and so on. Section 5 discusses the optimization method under workspace constraints, and study cases are discussed in Section 6. Finally, conclusions are drawn in Section 7.

Configuration and Actuation Mode
The 3-DOF compliant tensegrity manipulator discussed in this paper is constructed by stacking two SCs, as shown in Figure 1. These two SCs share the member l 26 for compactness and high stiffness, which is different from the manipulator introduced in Ref. [24]. The proposed mechanism consists of seven passive linear springs, two passive bars, two actuated cables, and two actuated bars. Linear springs are essential tensile elements in TM that provide compliance. The length of passive bars l 25 and l 36 are equal to have symmetric structure. Node A 0 is fixed at the origin of coordinate system A 0 -xy, and nodes A 1 and A 7 are connected with the base by prismatic joints. Analogously, end-effector g is connected to nodes A 3 , A 4 , and A 5 . Cables and bars are generally considered rigid in tensegrity mechanisms to simplify the modeling due to their relatively high axial stiffness compared to the springs. However, in practical implementations where non-metallic material (e.g., plastic) is used to reduce weight, the deformation of cables and bars cannot be ignored. Therefore, the compliance of cables and bars is considered in this work. Besides, all joints are frictionless and the effect of gravity is ignored. In order to obtain determined planar motion and keep all the springs tensioned, four actuators are used. However, in the inverse kinematics analysis, the end-effector g = [x 4 , y 4 , θ] T only provides three constraining equations, making the mechanism actuationredundant. Therefore cables l 23 and l 56 are designed to actuate synchronously so that the mechanism can have determined motion. The actuation strategy for l 23 and l 56 is given as where l c is constant. Using the proposed actuation strategy, the forward and inverse kinematics can have smooth motion and stable solution, respectively.

Step-1: Energy-Based Preliminary Modeling
For traditional mechanisms, the number of drives needs to be greater than or equal to the DOFs of the mechanism to ensure the unique motion. However, for the TM with springs, the spring can not be regarded as a constraint. The mechanism maintains equilibrium under springs and can deviate from the equilibrium position under external force. This is also the reason why such mechanisms are compliant. When cables and bars are assumed to be rigid, the mechanism presents three unconstrained DOFs if the actuation variables q a = l 16 l 27 l 23 T are given. The generalized coordinates of unconstrained DOFs are selected as q u = θ 1 θ 2 l 17 T , as shown in Figure 2. Based on geometric conditions, the coordinates of the nodes can be obtained as According to the Cosine theorem in triangles A 2 A 5 A 6 , the angle θ 4 can be written as Then the coordinates of A 5 is expressed as with In the same way, it is easy to calculate the coordinates of A 3 .
As the coordinates of all nodes are expressed, the potential energy stored in springs can hence be written as where k ij is the stiffness of spring A i A j , and l 0,ij represents the free-length of the spring. In the energy method, the equilibrium state of the tensegrity mechanism is the minimum potential energy state. Accordingly, it can be obtained by solving the following equations Equations (3) and (5) require that A 2 and A 6 cannot coincide, otherwise, the mechanism will have a singular configuration. In addition, when the two groups of rods are collinear, the second derivative of potential energy is less than 0, which is an unstable equilibrium state, resulting in a singularity of the mechanism too. Without considering the symmetry factor, the typical singular configuration is shown in Figure 3, and the singularity needs to be judged when solving.

Step-2: FDM Equilibrium Equations
The equilibrium equations of FDM are essentially the equilibrium equations of node forces. This method originates from the form-finding problem of cable-nets and can also be applied to the tensegrity mechanism. There are 8 nodes and 13 members in the proposed mechanism, and a 13 × 8 connectivity matrix C can describe the topology information of the TM. Supposing that member m(m = 1, 2, . . . , 13) connects node A i and node A j (i, j = 0, 1, . . . , 7). The components of the m-th row of C are defined as Let x, y ∈ R 8 denote the vectors of nodal coordinates in x-, y-directions, respectively. The coordinate difference vectors (member vectors) u, v ∈ R 13 expressed in a matrix-vector form as follow Based on u and v, the lengths of members can be obtained as Springs, bars, and cables are assumed to deform linearly. Therefore, given the free length l 0,m and axial stiffness k m , the force density of member m can be calculated as Then, the 8 × 8 force density matrix E is defined as [21] where the diagnosed matrix Q is expressed as If there is no external force acting on the mechanism, the equilibrium equations are Note that E can be expressed by x and y. Equation (14) is a set of non-linear equations with respect to nodal coordinates. The geometric constraints considering the fixed constraint and prismatic joints can be written as Meanwhile, in order to avoid the mirrored configuration, the coordinate range should be constrained by

Kinematics Solving Procedure
As the dimension of q a and the end-effector g are all three, there is only one unique stable solution within a reasonable range after excluding symmetrical results. According to the two-step method, both forward and inverse kinematics require solving nonlinear equations twice, as shown in Figure 4. Iterative methods, such as the Levenberg-Marquardt algorithm, are suitable for solving these equations.

Stability
Eigenvalues of the stiffness matrix can indicate the stability of a tensegrity mechanism. Denoting X = [x T , y T ] T as the generalized coordinates vector and the 16 × 16 global stiffness matrix K can be computed as [21] with where I 2 is a 2 × 2 identity matrix and ⊗ denotes the tensor product. This formulation considers the stressed equilibrium state as the reference state. The stiffness matrix in Equation (17) is a global matrix with respect to all coordinates. Removing rows and columns related to fixed coordinates (x 0 , y 0 , y 1 , and y 7 ), the 12 × 12 stiffness matrix K 1 of all moveable nodal coordinates can be obtained. Let λ = [λ 1 , λ 2 , . . . , λ 12 ] T denote the eigenvalues of K 1 , and the TM is considered stable when all the eigenvalues are positive, as Stability is the basic characteristic of tensegrity. If there are non-positive eigenvalues, it indicates that there is singularity caused by the collinearity of members.

End-Effector Workspace
Except for the limitation of the actuators, spring length is another factor that limits the workspace of the end-effector. The springs need to be elongated to ensure that the TM is tensioned and the length cannot exceed the elastic elongation limitation. At present, the conventional method of obtaining the workspace is discretizing q a or g and then solving the kinematics equations. The workspace of the end-effector considering geometric constraints, stability, and spring lengths, is defined as where Φ 1 and Φ 2 are geometric constraints, λ denote the eigenvalues of the stiffness matrix K 1 , and L S denotes the lengths of springs constrained within L low S , L up S . Cable stress, failure, and buckling of bars are not considered here.

Manipulability
The 3 × 3 kinematic Jacobian matrix maps the relationship between the actuation velocities and the end-effector velocities, which is a valuable means to evaluate the manipulability. It can also be used to determine whether there are singularity states.
The centered finite differences method is used to compute J, and the manipulability measure is expressed as

Potential Energy and Stiffness of End-Effector
Generally, a high level of potential energy is expected to resist deformation when relatively large external loads are applied. However, a low level of potential energy suggests that the manipulator is more mobile and safer in specific situations, e.g., using it as a surgical tool. In practical applications, the proportional adjustment of spring stiffness can directly change the potential energy.
In order to quantify the stiffness intuitively, the translational stiffness (normal stiffness, tangential stiffness) and rotational stiffness of the end-effector can be calculated based on the stiffness matrix. The tangential stiffness of the end-effector can be obtained by applying unit force along the tangential direction (l 45 ) on node A 4 , solving the node displacements using K 1 , and calculating stiffness. The calculation of normal and rotational stiffness is similar to that of tangential stiffness.

Optimization Model
For a designated task workspace, it is not easy to determine the installation location and design member lengths. In some situations, we need to design the TM as compact as possible. Therefore, we develop an optimization method to achieve the most compact design of the mechanism under workspace constraints.
The optimization model is defined as where [x 0 , y 0 ] is the installation location of A 0 , l p is the free-length of two passive bars, l c is the sum length of cables, l s is the free-length of springs (l 12 , l 17 , l 26 , l 35 , and l 67 ), l low a , l up a are the limitation of actuated bars, W t is the task workspace, which is determined by the actual requirements of the end-effector, W(x D ) is the workspace defined by Equation (20). l 23 change in [0.2 l c , 0.8 l c ]. The potential energy and stiffness can be changed without affecting the workspace by changing the spring stiffness in equal proportion. Therefore, they are not considered in the optimized model.

Constraint Simplification
The above problem is a seven-dimensional single-objective constrained optimization problem. As the objective function is simple, the core of this problem is how to transform complex workspace constraints (set constraints) into relatively simple numerical constraints. Using the discretization method to calculate the workspace defined Equation (20) is timeconsuming and unnecessary in the optimization process.
For the TM in this paper, its workspace is a three-dimensional simply connected space. Based on this, as long as the outer contour of the task workspace W t B is enclosed by W(x D ), it is safe to say that the entire task workspace is a subset of W(x D ), as Then, we discretize the contour W t B to obtain finite poses and reduce the computational cost.
Therefore, the constraint function can be defined as with In this way, the workspace constraint is expressed by a numerical constraint function, as illustrated in Figure 5. It is noteworthy that the number of discrete poses on the task workspace contour may influence constraint modeling accuracy. The fewer discrete poses we use, the faster the optimization algorithm runs, but it has the risk that the designed TM cannot reach all the poses, especially for a target workspace of irregular shape.

Optimization Procedure
As shown in Figure 6, after the task workspace is determined, the contour is extracted and discretized to get n poses. The constraint function and objective function value under given design variables can be obtained by solving n sets inverse kinematics equations. Considering that the constraint function of this problem is not expressed analytically, gradient-based optimization algorithms are not suitable, and these methods rely too much on initial values. This paper uses the adaptive simulated annealing (ASA), a global optimization algorithm, to solve this problem.

Workspace Analysis
In this section, we use a numerical example to analyze the workspace form. Two passive bars of the 3-DOF TM are set to 100 mm, and the length of the actuating bars varies in the interval [40, 160] mm. The sum length of two cables (l 23 + l 56 ) is set equal to 140 mm.
All bars and cables have the same axial stiffness of 1.0 × 10 2 N/mm. Three different springs are used in the tensegrity mechanism, as presented in Table 1. In order to better visualize the workspace boundary, we slice the workspace with respect to rotation angle θ, as shown in Figure 7.   Figure 7c,d, when θ = 0 • , the workspace is symmetrical with respect to the line x = 0 mm, and when θ = 30 • , the workspace is symmetrical about the line x = −17 mm. The original intention of the design is to use the cable-actuated SC(A 2 A 3 A 5 A 6 ) to produce angular (θ) and the bar-actuated SC(A 1 A 2 A 6 A 7 ) to produce displacement(x, y). However, the results show that there is a coupling between the two SCs. When the cable-actuated SC generates a specific θ, there is an x-direction displacement, which is the main reason for the offset.
The manipulability measure µ is distributed in the range [13.6, 26.9] with a mean value of 19.8, and there is no singularity in the whole workspace. On the whole, the distribution of µ has two trends. The first is that with the increase of |θ|, µ has a downward trend, significantly when θ = ±60 • , µ drops to near 14. The second trend is that under the same θ, the farther g is away from the origin, the larger µ is, as shown in Figure 7c.

Potential Energy and Stiffness of End-Effector
As shown in Figure 8a, the maximum and minimum potential energy are 1.0373 J and 0.2378 J, respectively, suggesting that the proposed TM requires at most 0.7995 J actuating energy to move even if there is no external load. At different θ, the potential energy distribution is similar, and it is mainly affected by the translation of the end-effector. As expected, when the actuated bars are both at the maximum length and l 23 = 0.5l c , the TM own maximum potential energy, as shown in Figure 8b.   Figure 8d. The tangential stiffness varies in [0.0134, 0.0441] N/mm, and its mean value is 0.0193 N/mm, as shown in Figure 8e. These results show that the tangential stiffness is significantly lower than the normal stiffness. One reason for this phenomenon is that the connections at A 0 and A 4 are equivalent to a series connection springs. The unevenness of stiffness in different directions is a drawback of this mechanism, and excessive tangential force must be avoided during application.

Effect of the Flexibility of Bars and Cables on Kinematics
By calculating the inverse kinematics under different bar/cable stiffness, we can understand the influence of bar/cable flexibility on kinematics and further evaluate the necessity of the proposed two-step method. The mechanism configuration, size, and springs are consistent with Section 6.1, and all bars and cables have the same axial stiffness. We can define the stiffness ratio as The effects of the flexibility of bars/cables can be quantitatively expressed as where q flex a is inverse kinematics solution of the two-step method, q rigid a is the inverse kinematics solution assuming bars and cables are rigid.
When θ = 30 • , the mean values of ∆ are shown in Figure 9. As expected, the larger the stiffness ratio of the bar to the spring, the smaller difference of the q a between the rigid and flexible models. The relationship between ∆ and η is close to linear in logarithmic coordinate. When η = 10, the average ∆ is 8.687 mm, and when η = 10 5 , the value falls into 8.132 × 10 −4 . Meanwhile, the distribution of ∆ is also related to g. Take η = 10 and θ = 30 • as an example. The maximum ∆ is 21.462 mm at g = [−26, 104, 30] T and minimum ∆ is 4.044 mm at g = [−26, 200, 30] T . Based on the comparison results, if the ratio of the bar stiffness and spring stiffness is less than 10 2 , the flexibility of the bars and cables have a relatively noticeable effect on the kinematics. It is necessary to use FDM for accurate modeling. When the stiffness ratio is more than 10 5 , it will not produce an obvious error if the bars/cables are assumed rigid.

Optimization Examples
Two examples are used here to validate the proposed optimization method, and the initial design is presented in Section 6.1 with F obj = 640 mm.
The first example is to design a TM that can rotate around a fixed point, and this remote-center-rotation constraint is modeled as This remote-center-rotation manipulator is usually used in medical applications [25,26]. In this situation, the task workspace is a smooth curve, and we have W t B = W t . Using the proposed contour sampling method in Section 5.2, we use eight sample poses evenly distributed on the ask workspace to model the constraint.
After generating 4947 designs by ASA, the objective F obj declined by 62.3%, from 640 mm to 241.3 mm, as shown in Figure 10a. The optimized workspace(θ ∈ [0, 30] • ) and task workspace are displayed in Figure 10b,c. Two endpoints of the curve almost fall on the contour of the optimized workspace, which means that a relatively optimal design has been found. Figure 11 are manipulability and potential energy plots of the optimized configuration. Compared with the initial configuration, the maximum manipulability increased significantly, but mainly concentrated in the edge region. The potential energy distribution is very similar to the initial design. At the same time, the decline of the overall potential energy level indicates that the elastic potential energy stored by the mechanism itself will also decrease after the overall size is reduced. Another example is to design a TM with workspace expressed as The workspace is an elliptical cylinder region in shape. Analogously, 48 poses on the contour are used to model the constraint. After generating 4997 designs by ASA, the objective F obj declined by 69.4 %, from 640 mm to 195.7 mm, as shown in Figure 12a. As shown in Figure 12b

Conclusions
This paper introduces a 3-DOF planar tensegrity manipulator with a compact structure constructed by two SCs. An advanced two-step kinematic modeling method is proposed to analyze its kinematics characteristics. The first step of this method is to estimate the equilibrium state of TM using the energy method, assuming that the bars and cables are rigid. Then, considering the flexibility of bars and cables, the second step is to obtain the accurate pose by solving the FDM equations using the state of the previous step as the initial design of iteration. After that, a sampling-based optimization method under workspace constraint is proposed. In this method, we use discrete points on the contour of a given target workspace to model constraints numerically. The following conclusions are obtained:

1.
The numerical example suggests that the TM has a relatively sizeable 3-DOF workspace and no singularity under the constraints of driving variables and spring elongation limitation; 2.
The stiffness of the end-effector in the tangential direction is significantly less than that in the normal direction, and excessive tangential force must be avoided during application; 3.
The flexibility of the bars and cables have a relatively noticeable effect on the kinematics when the bars-to-spring stiffness ratio is less than 10 3 , and it is necessary to use the two-step method for accurate modeling. When the ratio is more extensive than 10 5 , flexibility can be ignored.

4.
Two study cases are used to validate the effectiveness of the optimization method. After optimization, the objectives are reduced by 62.3% and 69.4% , respectively. Moreover, the target poses are located on and tangent to the contour of the workspace, suggesting that optimal designs considering compactness are obtained.