A General Model for Both Shape Control and Locomotion Control of Tensegrity Systems

Tensegrity systems composed of tension and compression elements have the potential for use in configurable structures and locomotive robots. In this work, we propose a general mathematical model for controllable tensegrity structures. Additionally, a method combining a genetic algorithm (GA) and dynamic relaxation method (DRM) is developed to solve the model. Our proposed model and method are applied to a typical shape controlled tensegrity and a typical locomotive tensegrity system. Firstly, the shape control of a two-stage tri-prism tensegrity is considered, and a collision-free path with minimum energy consumption is identified by using our approach. Secondly, gait design and path planning of a six-strut tensegrity is considered, and optimal gaits and motion paths are obtained by using our approach. The generality and feasibility of the proposed approach is conceptually verified in these implementations.


INTRODUCTION
Tensegrity systems are a class of special prestressed pin-jointed bar assembly, composed of compression elements and tension elements. A distinguishing feature of tensegrity structures is that their shape and mechanical properties can be actively controlled by prestressing their structural elements, making them good candidates for structural systems requiring controllable shapes and mechanical properties, such as smart structures (Shea et al., 2002;Fest et al., 2003;Motro, 2003;Al Sabouni-Zawadzka, 2014), deployable structures (Fazli and Abedian, 2011;Veuve et al., 2015;Kan et al., 2018;Oh et al., 2018), tunable metamaterials (Fraternali et al., 2014;Amendola et al., 2018;Liu et al., 2019), and robots (Paul et al., 2006;Mirletz et al., 2015;Kim et al., 2017;Cera and Agogino, 2018;Park et al., 2019;Wang et al., 2019). The actuators of a controllable tensegrity structure can usually be idealized as active members with variable rest lengths. The main problem to be solved in the control of a tensegrity system given an actuation configuration is determining the actuations (i.e., rest length changes of the active members) required to drive the structural system from its initial state to the target state.
If the actuations are imposed on the structural system very slowly and the structural system remains in static equilibrium throughout the control process, it is deemed a quasi-static system and the control problem can be interpreted as that of finding a static equilibrated path to connect the initial and target states. Most recent studies on shape control of tensegrity structures have treated it as a quasi-static process (Shea et al., 2002;Sultan et al., 2002;Xu et al., 2014). In particular, the prestressable equilibrium manifold of symmetrical prism tensegrity structures can be identified analytically, and then feasible control paths can be determined based on the equilibrium manifold using a search strategy (Sultan et al., 2002;Sultan and Skelton, 2003). For general cases in which the equilibrium manifold cannot easily be determined, an algorithm based on rapidlyexploring random trees has been proposed to find feasible actuation sets for the shape control of tensegrity systems (Xu et al., 2014). Form-finding methods, such as dynamic relaxation method (Fest et al., 2003) and non-linear force method (Xu and Luo, 2009) have been used to track the quasi-static motion of shape controlled tensegrities. Meanwhile, the dynamic effect have to be considered in the studies of locomotive tensegrity systems, due to the stronger actuations and environmental interactions involved in such systems. The potential for using tensegrity structures as locomotive systems has recently attracted considerable attention. A tensegrity swimmer was developed to achieve propulsive performance with closed-loop control (Bliss et al., 2013). Omer et al. (2011) proposed a 2D tensegrity robot to mimic caterpillar locomotion. Böhm and Zimmermann (2013) proposed a vibration-driven mobile tensegrity robot. The DuCTT (Duct Climbing Tetrahedral Tensegrity) with the ability to traverse complex duct systems was presented and demonstrated (Friesen et al., 2014(Friesen et al., , 2016. Spherical tensegrity robots with potential application in planetary explorers (SunSpiral et al., 2013;Sabelhaus et al., 2015) have been most intensively studied (Khazanov et al., 2013;Kim, 2016;Luo and Liu, 2017;Lu et al., 2019). To track the dynamic motion of locomotive tensegrities, Runge-Kutta method (Rovira and Tur, 2009), multi-body kinematic and dynamic simulation (Lin et al., 2016), and commercial physical engine (Zhao et al., 2017) have been used. Different descriptions and formulations are usually used for the shape control and the locomotion control problems of tensegrity systems, due to the difference in their application scenarios and the difference in the academic background of the researchers. In this paper, both the shape controlled tensegrity systems and the locomotive tensegrity systems will be modeled by the same mathematical formulations and a DRM-based motion tracking algorithm applicable for both quasi-static and dynamic motions will be adopted.
This paper is laid out as follows: section Model for Controllable Tensegrity Systems presents the general mathematical model for controllable tensegrity systems. In section Control Optimization Method, a method incorporating a GA-based optimization scheme with a DRM-based motion tracking algorithm is developed to solve the model. The proposed model and method are implemented in a tensegrity system shape control application in section Shape Control. Section Locomotion Control presents the application of the proposed model and method to locomotive control in tensegrity systems. Finally, section Discussion concludes the study.

MODEL FOR CONTROLLABLE TENSEGRITY SYSTEMS
For a tensegrity system composed of n nodes and n c elements, an element's type is given by the element type vector D = [D 1 , D 2 , . . . , D n c ] where D i is defined as: For a controllable tensegrity system with n a (n a ≤ n c ) active elements, the locations of the active elements are described by an active element location vector in which "active element" means the element whose length can be actively changed by actuators, and "passive element" means the element whose length cannot be actively changed. The active element location vector satisfies that sum (D A ) = n a . The control strategy of a controllable tensegrity system can be described by an actuating function . During a given time period t s , t f , the actuating function of the i th member is defined as e i = f i (t), where f i is a continuous function of time t. Then, the actuating function is written as = e 1 , e 2 , . . . , e n c where e i ∈ [−e s i , e l i ], in which e s i and e l i are the allowable shortening and elongation for the i th member, respectively. Note that e s i = e l i = 0 for passive members (D Ai = 0). The rest length change ranges of the elements can be expressed as E = −e s 1 , e l 1 , −e s 2 , e l 2 , . . . , −e s n c , e l n c . Hence, at any time, the rest length change function of the controllable tensegrity system is required to satisfy the condition that ∈ E.
The rest length of members L t R at time t satisfies The internal force F t i of the i th member is determined by where E i , A i , L t Gi , and L t Ri are the elasticity modulus, crosssectional area, deformed length, and rest length of the i th member, respectively. Meanwhile, the internal force F t i of the i th member is required to be within the strength limitations, i.e., Frontiers in Built Environment | www.frontiersin.org where F l i and F u i are the lower bound and upper bound of the strength of the i th member, respectively. Specifically, F l i and F u i are given by where F ci , F bi , and F ti are the compression strength, buckling strength, and tension strength of the i th member, respectively.
Having defined f = F l 1 , F u 1 , F l 2 , F u 2 , . . . , F l n c , F u n c , the strength constraint on the members of the controllable tensegrity system can be written as F t ∈ f .
The equilibrium equation of the tensegrity system at time t is written as where A t (n r × n c ) is the equilibrium matrix, F t e (n r × 1) is the external nodal load vector, and n r is the number of degrees of freedom. Hence, the unbalanced nodal force R t of the system at time t can be expressed as Then, according to Newton's second law: where M is the nodal mass matrix of the system, C is the viscous damping of the system, and U is the matrix of nodal coordinates of the system. Substituting Equation (9) into Equation (10) yields The control objective of the system is case-dependent. It could be a special requirement of the nodal displacement, locomotion distance, or the energy cost to generate the required shape adjustment or motion. Since it is usually a function of the control strategy, it can be conceptually represented as G( ). Based on the above definitions of the control strategy, constraints, and control objective, a general model for the controllable tensegrity system can be written as where s represents free space in an environment where the controllable system does not interfere with boundaries or obstacles, and C( ) represents the additional constraints on the system in specific situations.

Motion Tracking
An incremental procedure based on the dynamic relaxation method (Xu and Luo, 2013) is adopted to simulate the motion of the controllable tensegrity system under a given actuation. The duration of the actuation is discretized by a given time increment of t. The state of the system is tracked by a DRM procedure that starts from a known state at time t and ends with a new state at time t + t after a time increment of t.
Step by step, the motion of the system under the given actuation is discretely depicted. Details of the DRM-based motion tracking procedure are given in the following paragraph.
Under a given actuation t , the residual force R t of the system is determined by Equation (9). According to Equation (10), for a node j in the direction x at time t, we have where R t jx represents the residual force of node j in the direction x, and is a component of R t ; M j represents the equivalent nodal mass of node j; and C jx represents the damping of node j in the direction x. The nodal acceleration can be approximated in centered finite difference form as Substituting Equation (14) into Equation (13), the nodal velocity at t + t/2 can be expressed aṡ Then, the nodal coordinate at t + t can be expressed as By setting t = t + t, the actuation t , equilibrium matrix A t , internal force F t , external force F t e , and residual force R t can be updated accordingly. The above central difference process repeats until t = t f . The numerical stability of this central difference process is guaranteed by using a time increment t no larger than 2M S , where S is the highest direct stiffness of any node relative to adjacent nodes (Barnes, 1999). For quasistatic problems, the motion tracking can be interpreted as a series of form-finding processes (Xu and Luo, 2013). In each form-finding process, a sub-iteration scheme is implemented whereby each time step is iterated so that convergence is achieved for each time step (Senatore and Piker, 2015). Fictitious values for mass, stiffness, and damping are used to improve the performance of the algorithm for convergence to the static solution (Zhang et al., 2006). In particular, a technique called "kinetic damping" is usually adopted to expedite the convergence of algorithm (Barnes, 1999;Zhang et al., 2006). Note that the real value for stiffness should be substituted into the each form-found system to evaluate whether the stress and stability limits of the structural elements are met. For dynamic problems, real values for mass, stiffness, and damping are used by the algorithm. The effectiveness of it in dealing with dynamic problems has been demonstrated by a similar procedure called "finite particle method" (Yu and Luo, 2009;Yu et al., 2011). Note that a DRM-based scheme for both quasistatic analysis and dynamic analysis also has been proposed by Senatore and Piker (2015).
An internal collision-detecting strategy proposed in the literature (Xu et al., 2014) is used here to check whether there is internal collision between any pair of elements. This must be avoided, since any internal collision may cause the moving system to become stuck. As a result, the motion tracking algorithm stops if an internal collision is detected, and the path and corresponding actuation are deemed infeasible. On the other hand, collision between the structural system and its environment is allowed and must be accounted for. The penalty FIGURE 1 | Flowchart of the solving method.
Frontiers in Built Environment | www.frontiersin.org function method (Oden and Kikuchi, 1982;Papadopoulos and Taylor, 1992) is used to consider the interaction between the structural system and its environment. An obstacle in the environment is assumed to be a series of surfaces covered with springs, and the nodes of the tensegrity structure are allowed to penetrate the surface. The normal contact force F n between the node and the surface is evaluated by where δ is the penetration depth; and k n is the virtual normal stiffness of the surface. The friction force F f between the node and the surface is determined by where µ is the friction coefficient.

Optimization Algorithm
To find a control strategy that minimizes the given objective, an optimization algorithm is required. The control problem modeled by Equation (12) is essentially a problem of combinatorial optimization, and the search domain is usually very large. Direct search algorithms such as evolutional algorithms, simulated annealing algorithms, and genetic algorithms are considered suitable for this problem. Without loss of generality, a genetic algorithm is adopted here.

Solution Flowchart
Combining the DRM-based motion tracking method and GAbased optimization algorithm, the flowchart of the algorithm is as shown in Figure 1. The main steps of the algorithm are as follows: (1) Define the genetic representation and the fitness function; initialize a population.
(3) For each individual, build a numerical model of the corresponding tensegrity system for DRM; initialize the parameters of the numerical tensegrity model, and set t = t s ; perform Step (4) to Step (6). (4) Apply a control strategy to the numerical tensegrity model through the actuating function t . (5) Calculate the residual nodal force R t ; solve Equation (11) by the central difference method; update the nodal accelerationÜ t , nodal velocityU t , and nodal coordinate U t ; apply constraints. (6) If t = t f , go to Step (7); otherwise, set t = t + t, and go to Step (4). (7) Evaluate the fitness of the population. If the population satisfies the convergence condition or the number of generations reaches the maximum permissible, go to Step (9); otherwise, go to Step (8). (8) Apply the selection, crossover, and mutation operators, and create a new generation. Then, go to Step (2). (9) Output the solutions.

Shape-Controlled Tensegrity System
The feasibility of the proposed model and method in solving tensegrity system shape control problems is tested in a doublelayered tri-prism tensegrity structure, which is used as a typical example of a shape-controlled tensegrity system in a previous study (Xu et al., 2014). As shown in Figure 2, the shape-controlled tensegrity consists of 12 nodes, 6 compression members, and 18 tension members. The initial locations and target locations of the nodes are given in Table 1. The tension members comprise six vertical cables, six diagonal cables, and six saddle cables. The connectivities of the elements are also shown in Figure 2. The triangle bottom and top of the tensegrity structure are assumed to be rigid, and the three bottom nodes is fixed to the ground. There is a ceiling at z = 2.830 m in the space. Properties of the members used for the shape-controlled tensegrity structure are given in Table 2. A fictitious nodal mass of 1.00 is assumed because the dynamics of the system was not considered in this shape control problem. A linear stiffness of 2 × 10 5 N/m for struts and a linear stiffness of 200 N/m for cables are assumed. The "kinetic damping" is used to expedite the convergence of the DRM. The GA's parameters are set as follows: The population size is 50; the selection type is stochastic universal sampling; the crossover type is two-point; the crossover probability is 1; the mutation probability is 0.7 divided by bits of coding; the maximum number of generations is 200. The algorithm stops when the number of generations reach the maximum number specified.
Vertical and diagonal cables are used as active elements whose rest lengths can actively change in the ranges of 0.40-3.20 m and 0.05-4.00 m, respectively. The control strategy represented by the actuating function is simplified using the following approach: The actuating function for an active element is divided into a series of actuating steps; in each step the rest length of the active element is assumed to be linearly changed. Given the initial rest length, initial time, end time, and the rest length at the end of each step, the linear actuating function of the element in each step is easily determined. Denoting the rest length change of the element i in the k th step as e ik , the actuating function of an element with q steps can be expressed as e i (e i1 , e i2 , . . . , e iq ). Further assuming that all the active elements possess the same division of actuating steps, i.e., all have q steps and the start time and end time of each step are the same for all active elements, the actuating function of the system can be expressed as {e 1 , e 2 , . . . , e nc }.

Path Optimization
The requirement to be collision-free is considered as an additional constraint, expressed as S∈ c , in which S represents the state of the system and c represents the states without internal collision. The target state of the control is imposed as a further constraint that U tf = U target . The target state, i.e., the target location of nodes, is given in Table 1. This corresponds to change the center of the top triangle of the double-layered tensegrity from the initial position (0.000, 0.000, and 2.534) to (1.857, 0.249, and 2.065). The energy cost of the control strategy is used as the objective which can be expressed as: where n c is the number of elements; F t i and F t−∆t i are the internal forces of the i th element at times t and t-∆t, respectively; and e t i and e t−∆t i are the rest length changes of the i th element at times t and t-∆t, respectively. The control model is rewritten as The reciprocal of the objective, i.e., 1/G( ), is selected as the fitness for the GA. The number of actuating steps is assumed to be 2, i.e., q = 2. The computation is carried out on a personal computer with a Intel Core i7-4790k CPU @ 4.0 GHz and 16 GB RAM. It has taken 293 min to carry out the 200-generation revolution. It is observed that the total energy assumption decreases rapidly in the first 46 generations and then converges slowly on 258.17 N·m (Figure 3). The optimal control strategy, corresponding to the minimum energy cost, is shown in Table 3, and the equilibrium configuration at the end of each actuating step of the optimal control strategy is shown in Figure 4. The motion trajectory of the center of upper triangle during the control is shown in Figure 5, in which the thick blue line represents the optimal path generated by the optimal control strategy, and the thin red lines represent 10 feasible paths generated by 10 feasible control strategies found in the evolution process. It is observed that the motion path generated by the optimal control strategy is much shorter than those generated by the unoptimized control strategies.

Locomotive Tensegrity System
The application of the proposed model and method to locomotion control in tensegrity systems is tested in a sixstrut locomotive tensegrity structure. As shown in Figure 6, the tensegrity system comprises 12 nodes, 6 compression elements, and 24 tension elements. The connectivities and numbers of the   elements are also shown in Figure 6. The 24 tension elements generate a pattern of 20 triangles on the outer surface of the structure. Eight of the surface triangles, each of which is closed by elements, are called closed triangles (TC). The other triangles are called opened triangles (TO) because each has only two tension elements and one edge of the triangle is open. For instance, as shown in Figure 6, the triangle with nodes 3, 7, and 12 is a TC, and the triangle with nodes 3, 11, and 12 is a TO. In the ideal state without external load and gravity, the compression elements and the tension elements have the same geometrical length. As a result, in the ideal state all the TCs and TOs are identical. The six compression elements are used as active elements, with an actuating range of [−5, 5] cm. The actuating speed is assumed to be 10 mm/s. The rest lengths of elements at the initial state and the mechanical properties of the elements are as shown in Table 4. The initial state of the system is obtained by applying gravity to the ideal state. A damping coefficient of 0.01 is assumed for the elements. Environmental parameters are set as follows: A stiffness of 1,000 N/m, and a friction coefficient of 0.5 are assumed for the ground; the gravitational acceleration is −9.8 m/s 2 .

Gait Design
The gait of the locomotive tensegrity structure is designed to satisfy the conditions that (1) the system acquires a certain amount of motion after the gait; and (2) the rest lengths of active elements are not changed before and after the gait, i.e., The motion of the system can be expressed by the displacement of the mass center as u t c = ||u t c -u ts c ||, in which u c = [u cx , u cy , u cz ] is the coordinate of the mass center.
Regarding the symmetry of the tensegrity system, there are two basic stand states in the locomotive system: the TC state and the TO state. As indicated by the names, the tensegrity structure stands on the ground with a TC in the TC state, and with a TO in the TO state.
There are two typical types of gaits for the six-strut locomotive tensegrity system, namely crawling gaits, and rolling gaits . It can be easily recognized that rolling gaits usually generate a larger motion than crawling gaits. To obtain rolling gaits rather than crawling gaits, the objective function for gait design is defined to maximize the horizontal motion of the mass center of the system, i.e., where u t f cxy is the horizontal motion distance of the mass center at time t f ; u t s cx and u t s cy are the x and y coordinates of the mass center at time t s ; and u t f cx and u t f cy are the x and y coordinates of the mass center at time t f . An additional constraint on the number of active elements used in the gait is adopted to impose control on the number of active elements used in the motion. The additional constraint is written as where n a is the number of active elements used in the gait, and N a is an integer constant no larger than 6, specified by the designer. Substituting the above objective in Equation (22) and additional constraints Equations (21) and (23) into Equation (12), an optimization model for the gait design can be obtained: The method proposed in this paper is used to solve this model. The GA used here adopts a maximum number of generations of 100, and the other parameters of it are same as those used in section Shape Control. The computations are carried out with the      Figure 7. In each case, six computations, corresponding to N a = 1, 2, 3, 4, 5, and 6, are performed. The evolution process of the mass center's horizontal motion distance in the computations are shown in Figure 8. In the case using the TC state as the initial state, three typical gaits that, respectively, drive the system rolling from the current stand triangle to the three adjacent triangles are found ( Table 5). The motion distance of the mass center generated by these gaits is in the range 6.40-6.50 cm. In the case using the TO state as the initial state, there are also three gait types. Two respectively, drive the system rolling from the initial stand triangle to two of the three adjacent triangles, whereas the third drives the system rolling twice, from the initial stand triangle to an adjacent triangle, and then to an adjacent triangle again ( Table 5). The motion distance of the mass center generated by the single-rolling gaits is 6.24 cm. The motion distance of the mass center generated by the doublerolling gait is 10.48 cm; about 1.68 times that generated by the single-rolling gaits.

Motion Path Planning
For a given initial and target location pairing, there are generally two approaches to realizing motion from the initial to the target location. One is connecting the initial location to the target location using several basic gaits, as given in Table 5. This approach is suitable for regular movement on flat terrains with given maps, and a geometry-based algorithm has previously been verified for finding the path for locomotion of a six-strut tensegrity robot (Lu et al., 2019). The other approach is a direct search for a control strategy able to drive movement of the system from the initial location to the target location. This approach is considered more flexible than the first, but normally carries greater computational cost because it has to numerically track motion within the path planning scheme. In the current work, the second approach is adopted. The locations A(x A , y A ) and B(x B , y B ) are defined as the initial location and the target location. The objective function is defined as: where d is the distance between the current location and the target location, and ε is a positive small constant. The aforementioned additional constraint of no internal collision is included. The motion path planning model can thus be written as The initial location is set as A (0, 0) and the target location as B (0.1, 0.1). The starting time is set as t s = 0, the end time as t f = 150 s, and the positive small constant ε as 0.0005 m. Time is discretized by an interval of 1 s. The GA used here adopts the same set of parameters as used in section Gait Design. It has taken 1,750 min to complete 75 generations on the same computer as that used for the motion path planning optimization. The evolution of the average and minimum objectives in the computation is shown in Figure 9. The minimum distance between the current and the target location decreases to 0.0004 m in the 69th generation. The individual corresponding to the minimum objective in the last generation is the optimal control strategy, represented by the rest length change spectra of the six active elements, as shown in Figure 10A. The motion trajectory of the system's mass center under the optimal control strategy is shown in Figure 10B. Further investigation into the motion reveals that the system experiences three complete rolls and the stand triangle shifts from TC-6 (3, 7, 12) to 4,7), 7,9), and finally 8,9). Besides the rolling movements, some crawls and shape adjustments occurs, resulting in small-scale movements, as indicated by the red dashed-line rectangles in Figure 10B.

DISCUSSION
In this paper, the generality of the proposed model is ensured by the following factors: It does not deliberately distinguish between the shape-controlled tensegrity structure and the locomotive tensegrity structure. A unified dynamic formulation is used for the controllable tensegrity structure, though the shape control process is usually deemed to be quasi-static. The DRM-based algorithm using the central difference explicit scheme can track both the deformations and global rigid body motions of tensegrity systems. Conceptual expressions of the subjective and constraints were adopted. The control problem for a controllable tensegrity structure was interpreted as a combinational optimization problem. A similar optimization model was used for shape-controlled tensegrity structures in previous studies (Xu and Luo, 2009;Xu et al., 2014). However, in the current work, quasi-static behavior was assumed for the structural system, and the global rigid body motion of the structural system was eliminated by proper restraints. To some extent, the proposed model may be considered an extension of the previous model.
The global optimization algorithm used in the paper, i.e., the GA, worked well in the shape control example in section Shape Control and the locomotion control example in section Locomotion Control. In both examples, the number of control steps was provided in advance. Situations in which Frontiers in Built Environment | www.frontiersin.org the number of control steps cannot be determined in advance, as in the path planning situations considered in the literature (Xu et al., 2014;Lu et al., 2019), are difficult to handle with a GA. Further considering that there is no constraint on the number of control steps in the mathematical model given in Equation (12), the generality of the GA does not match that of the model. Various other global optimization algorithms besides a GA could be used, based on the specific application.

CONCLUSION
In this paper, a general model for both shape control and locomotion control of tensegrity systems is proposed. A method combining GA and DRM was adopted to solve the model. The proposed model and method were applied to typical tensegrity system shape and locomotion control problems. The results demonstrate that both the shape control problem of a doublelayered tri-prism tensegrity structure, and the locomotion control problem of a six-strut tensegrity structure can be modeled by the proposed model, and solved by the proposed method. We believe the proposed model is also suitable for control problems in other types of tensegrity system. The applicability of the proposed GA-based method to control problems of higher complexity and larger scale than those presented in this work could be subject of future investigation.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.