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 Rjx t 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 as

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 2MS, where S is the highest direct stiffness of any node relative to adjacent nodes (Barnes, 1999). For quasi-static 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 quasi-static 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 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.

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.

Combining the DRM-based motion tracking method and GA-based optimization algorithm, the flowchart of the algorithm is as shown in Figure 1. The main steps of the algorithm are as follows:

The feasibility of the proposed model and method in solving tensegrity system shape control problems is tested in a double-layered 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⁵ 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₁, e₂, …, e_nc}.

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; Fit and Fit-Δt are the internal forces of the i^th element at times t and t-Δt, respectively; and eit and eit-Δt 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.

The application of the proposed model and method to locomotion control in tensegrity systems is tested in a six-strut 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².

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 (Li et al., 2017). 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 Δucxytf is the horizontal motion distance of the mass center at time t_f; ucxts and ucyts are the x and y coordinates of the mass center at time t_s; and ucxtf and ucytf 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 same computer used in section Shape Control. Two cases, using the TC state and TO state, as the initial state, are considered. The number of stand triangles and triangles adjacent to the stand triangle are shown in 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 double-rolling gait is 10.48 cm; about 1.68 times that generated by the single-rolling gaits.

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 TO-7 (3, 4, 7), TC-8 (4, 7, 9), and finally TO-11 (7, 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.