In this paper we present a control methodology for robots made up of high-elongation linear actuators connected together at universal joints, which we call Linear Actuator Robots (LARs). Such robots can change their shape dramatically through the coordinated actuation of their linear members. The control of robots with a large number of degrees of freedom is important to allow robots to become increasingly flexible to a wide variety of tasks. In the case of a linear actuator robot, the robot can change shape to be better suited for a multitude of tasks including locomotion, manipulation, and matching of 3D target shapes (shape morphing). Such a robot would be valuable in search and rescue missions, where a LAR could flexibly maneuver over uneven terrain and morph into a custom manipulator to clear debris. LARs can also serve as a type of "programmable matter", changing shape to represent 3D objects and responding to a human designer's digital manipulations in real time.
In this work, we present nonlinear optimization techniques to enable a linear actuator robot to locomote. We present a differential kinematic analysis of LARs, relating the velocities of the nodes in the structure to the rate of change of the actuator lengths. This allows us to link concepts from graph rigidity to the control of the robot structure. We use this N. Usevitch Fig. 1: We present two algorithms for Linear Actuator Robots (LARs) to locomote. This figure shows snapshots of an optimized everting gait for a LAR with 10 edges and 5 vertices, computed with one of our algorithms. The LAR shown is a passive mockup made from 10 car antennas, and handpositioned to illustrate the gate. kinematic analysis to derive two on-line planning algorithms for locomotion which are both based on the same underlying nonlinear optimization algorithm tailored to the kinematics and constraints of LARs. A passive mock-up of a Linear Actuator Robot executing one of the optimized locomotion trajectories presented in this paper is shown in Fig. 1. In this case, the robot is composed of 10 actuators (passive car antenna elements), and 5 nodes (with spherical joints formed by magnets attached to steel balls).
Related Work: To best understand the related work, we first present an overview of different implementations and applications of linear actuator robots, and then discuss in detail the specific control methodologies. Early work on TETROBOTs proposed robots composed of linear actuators arranged in repeated graphical motifs of tetrahedrons or octahedrons to facilitate kinematic computations [1]- [3]. Robots of this type have been referred to as variable geometry trusses [4], and have been proposed as manipulators [5], as platforms that allows locomotion over various terrains, and as robots with the ability to change shape to adapt to tasks that may be unknown a priori [6]. Other physical variants of shape changing robots based on linear actuators include the modular linear actuator system presented in [7], an octahedron designed for burrowing tasks made of high extension actuators presented in [8], and an active-surface type device that uses prismatic joints to deform a surface into arbitrary shapes while respecting some constraints [9]. In [10] a user interface is presented that allows a novice user to create a large scale truss structure, and then animate its motion by inserting a few linear actuators. In [11] a 2D structure is built from a collection of triangles with a guide). For this reason, we propose that better quasistatic planning methods are valuable steps towards improved system performance. We also note that whereas the dynamics of tensegrity systems often play a large roll in their response, the linear actuators used in linear actuator robots are often relatively slow, leading to less emphasis on leveraging the dynamics. In this work we follow the precedent of the prior work on linear actuator robots and utilize a kinematic model.
In this paper we present the kinematics of linear actuator robots with arbitrary graphical structure, including overconstrained structures, and utilize an optimization based approach for planning directly over the position of nodes of the robot. This optimization approach allows our method to be customized to different tasks and cost functions. We consider actuator constraints (to enforce min-max elongation), physical constraints (to prevent self-intersection and enforce a minimum angle between connected actuators), and constraints to avoid kinematic singularities in arbitrary robots in designing locomotion algorithms. The key contributions of this work are two algorithms to solve the following problem:
Problem 1: Move the center of mass of the robot in a prescribed direction v cm , or along a prescribed trajectory x cm (τ ), ensuring that the robot is always physically feasible and that the robot does not pass through any singular configurations.
The first algorithm we propose solves an online optimization to minimize an objective function that considers only the current state and motion of the robot while ensuring physical feasibility. This method applies to any robot that is in an infinitesimally rigid configuration. However, this method does not guarantee the persistent feasibility of the robot's motion, meaning it is possible the robot will reach a configuration from which it cannot continue without violating physical constraints (i.e., it might get tangled up). To ensure persistent feasibility, we present another method where we solve an offline optimization that generates periodic motion primitives to move a robot from a starting configuration to an equivalent configuration centered on a new support polygon. This motion primitive is then used by a high-level planner to plan paths from an initial configuration to a goal. We refer to this method as the two-tiered planning approach. This method guarantees persistent feasibility of the trajectory, but requires that the initial configuration of the robot satisfy certain symmetry requirements. The performance of the two algorithms is compared in simulation study in which we find that the two-tiered planning approach gives better performance in terms of cost.
This paper builds upon our past work by the authors presented in [42] on modeling and providing an algorithmic foundation for this class of robots. This work adds singularity constraints and angle constraints to our past work, and the solution method for the online optimization method has been improved. The two-tiered planning algorithm using motion primitives is also new in the present work. We also note that the work presented in the conference paper has served as the foundation of the work in [23,43].
The rest of the paper is organized as follows: Sec. II formalizes a model for LARs and derives the forward and inverse kinematics relating the change in actuator lengths to node positions. Sec. III describes the physical constraints imposed to ensure the robot motion is feasible. Our single step locomotion algorithm is given in Sec. IV, and our twotiered approach is presented in Sec. V. These methods are compared in Sec. VI. In Sec. VII we discuss the performance of the kinematic plan in the presence of dynamic effects, and conclusions are given in Sec. VIII.
Formally, we model a Linear Actuator Robot (LAR) as a framework which consists of a graph G and vertex positions p i ∈ R d . Our methods are applicable to LARs embedded in Euclidean space of arbitrary dimension d, but we focus on the embeddings in 3D (d = 3). The graph is denoted as G = {V, E}, where V = {1, . . . , N } are the vertices of the graph, and E = {. . . , {i, j}, . . .} are the undirected edges of the graph. The geometry of the robot is fully represented by the concatenation of all vertex positions x = [p T 1 , p T 2 , ..., p T n ] T . In this paper we consider a quasistatic model as opposed to a dynamic model, implicitly assuming that the robot's motion is slow enough that inertial effects are negligible, an assumption that we will further address in Sec. VII. We define a length vector L, which is a concatenated vector of the lengths of all edges in the graph
The vector L is of length n L , equal to the number of edges of the graph, and can be directly computed from the framework (G, x). We use the notations L(x) to indicate the length vector induced by a set of node positions x. We note that the relationship in (1) is the constraint on the node positions created by an edge. Note that L(x) represents the "inverse kinematics" for LAR robots since it is a function that maps from the vertex positions (analogous to the end effector position in a serial manipulator) to the lengths of the linear actuators (analogous to the joint positions in a serial manipulator), and it is trivial to obtain (as also noted by [1]).
While it is trivial to obtain the edge lengths from the node positions, our task is to invert this relationship and control the node positions by changing the edge lengths. In a network of linear actuators, each link length imposes one constraint on the node positions as given in (1). Finding the vertex positions from the link lengths means finding node positions that satisfy all of the constraint equations up to translation and rotation of the entire network. Several classes of solutions exist based on the rigidity of the underlying graph. Examples of a few of these classes are shown in Fig. 3, and our analysis of the device kinematics in the following section will depend on the rigidity of the underlying graph of the robot. If the system of equations has infinite solutions the framework is not rigid, as it is possible to move the system relative to itself without violating length constraints as in Fig. 3(i). A framework is rigid if there are a discrete number of solutions to the constraint equations, and all deflections of the system relative to itself violate the length constraints.
Of particular use to our analysis are graphs that are infinitesimally rigid, meaning that all infinitesimal deflections of the system relative to itself violate the length constraints.
If the system is infinitesimally minimally rigid and a minimal set of constraints is applied that is linearly independent of the link constraints, the combined matrix H = [R T C T ] T is full rank and square, and hence invertible, allowing us to write
Note that this is the form of a driftless dynamical system, and that H(x) -1 is the Jacobian matrix relating the motion of the actuators to the motion of the nodes. The vector L describes the rates of change of the linear actuators, and hence is the input to the system. Equation (5) shows that when the H(x) is invertible, each input channel Lk can be commanded independently of the others. The fact that this matrix is invertible means that the input space is all possible length velocities, allowing us to make the following proposition:
Proposition 1: Given an infinitesimally minimally rigid framework with the minimum number of constraints to the environment, the length of each edge can change independently.
This means that it is not necessary to coordinate movements between lengths as long as the graph remains minimally infinitesimally rigid.
If the system is infinitesimally rigid but not minimally rigid, it is over-constrained and some motions of the linear actuators must be coordinated. In this case, the H matrix is skinny, with more rows than columns. Taking the singular value decomposition of the combined H matrix,
The bottom rows of this expression can be expressed as a constraint which encoding how certain lengths must move in a coordinated fashion:
By utilizing this constraint, redundant rows of the H matrix and their corresponding elements in the vector [ LT 0] T can be removed until it is square and full rank, and hence invertible. We call the reduced H matrix and L vector the master group, and we denote them as H m and L m respectively. We refer to the removed rows as the slave group, and denote as H s , and the removed actuator inputs as Ls . We note that the actuators chosen for the master and slave groups are partially up to the user's discretion, and could potentially change based on configuration. As an example procedure, an algorithm could initialize H m = H, and H s as an empty matrix, and then iterate through each row of the H m matrix. If it finds a row linearly dependent on the previous rows from the H m matrix and places it in H s , and removes the corresponding element of L m and places it in L s . This allows us to express the system as follows:
Due to (10) the input space is restricted such that only combinations of link velocities that satisfy the constraint can be physically realized. The master inputs Lm can be picked arbitrarily, but Ls must be chosen to satisfy the constraint equation.
This system can be expressed in the standard form of a linear dynamical system, ẋ = Ax + Bu where A = 0, u = [ LT 0 T ] T , and B = H m (x) -1 . We now make the following proposition:
Proposition 2: A framework that is infinitesimally rigid is fully actuated.
This means that for an infinitesimally rigid system control of every degree of freedom can be achieved given control of the rate of change of the actuator lengths and the motion of the contact points. This has the key advantage of allowing us to plan our motion in terms of node positions, and then use the [R T C T ] T matrix to determine what input to apply to the actuators.
Locomotion requires finding a method to actuate the robot to move while it maintains physical feasibility. We define feasibility as follows:
Definition 1: A framework (G, x) is feasible if it meets three types of physical constraints: (i) the lengths of all actuators fall within a fixed maximum and minimum length range, (ii) the actuators do not physically intersect (except at the endpoints of two connected actuators), and (iii) the angles defined by two actuators connected at a joint remain above a minimum value.
To ensure that all motions of the robot are physically feasible, we detail the form of the constraints and quantify how many of each type of constraint occurs in the optimization based on the characteristics of the underlying graph. In addition to these physical constraints, we also present constraints to prevent the robot from crossing configurations where infinitesimally rigidity is lost, which correspond to the singular configurations of the robot.
For physical feasibility to be preserved, all actuators must be maintained between a maximum and minimum actuator length. The squared length of actuator k that connects nodes {i, j} is quadratic in x, and the constraint that it remain within the set maximum and minimum length can be expressed as
where A k is a matrix where the only nonzero entries are
We note that constraints of the quadratic form x T Qx ≤ c, where c is a positive constant, are convex if and only if Q is positive semi-definite. We note that A k is the Laplacian matrix of a graph that contains only edge k. As the Laplacian matrix is always positive semi-definite, the maximum length constraint is convex in the node positions while the minimum length constraint is not. Thus our algorithms will handle non-convex and nonlinear constraints.
We also enforce the constraint that actuators do not collide physically, except for at the vertices where they are joined. To determine if two actuators cross, the minimum distance between them must be greater than d min , a positive diameter of the actuator assuming that the actuator can be represented as a cylinder. The minimum distance between actuators connecting vertices i, j and k, l is denoted as d kl ij , and can be expressed as follows:
) These links are not in collision if d kl ij > d min . Efficient algorithms for this computation have been explored previously [47]. Checking the pairwise distances between all edges in a graph requires checking
constraints of the type expressed in (12). As we do not compute the distance between actuators that are connected at a node, the number of constraints is reduced by the number of pairwise distances between edges that meet at a node, which for node i is given by
where g i is the degree of the node. Thus the total number of constraints to avoid collisions between actuators is
Another key physical constraint is that the angle between connected actuators remain above a certain value, which is especially important when the actuators have a high elongation ratio. An angle constraint between two edges is a function of 3 vertices. We define p i as the position of the shared node between two edges, and p j and p k as the other vertices of the two edges. The angle constraint is:
The number of angle constraints can also be expressed in terms of the degree of the nodes of the graph:
Our proposed optimization approach is based on the observation that if the robot is infinitesimally rigid, we can directly optimize a path for the node positions and recreate the needed actuator trajectories. For this assumption to remain valid, the robot must maintain its infinitesimal rigidity, meaning the rigidity matrix R must remain of rank 3n -6. Designing controllers that maintain infinitesimal rigidity has been a topic in formation control of multi-agent systems [48,49]. In [48] the rigidity eigenvalue for frameworks in R 3 is defined as the 7th smallest eigenvalue of R(x) T R(x) and the gradient of the rigidity eigenvalue with respect to the node positions is used as part of a controller. In the general case, infinitesimal rigidity can be enforced using the following constraint:
where λ 7 is the 7th smallest eigenvalue of the R(x) T R(x) matrix and λ crit is its minimum allowable value.
One problem with ( 16) is that the magnitude of λ 7 changes quadratically with network size. To provide a constraint that is invariant to network scale, we instead use the worst case rigidity metric, taken directly from [50], and defined as:
It has been noted that if a framework is infinitesimally rigid in one configuration it is infinitesimally rigid almost everywhere, meaning that for a graph with one infinitesimally rigid configuration, the set of non-rigid configurations is a set of zero measure [51]. We make the observation that the configurations where the robot loses rigidity often divide the state space into disconnected regions. We define each of these regions as a rigidity equivalence class as follows:
Definition 2: (Rigidity Equivalence Class) The framework F 1 = (G, X 1 ) and the framework F 2 = (G, X 2 ) are in the same rigidity equivalence class if a continuous path x(t) exists such that x(0) = X 1 , x(T ) = X 2 , and the rigidity matrix R(G, x(t)) is maximal rank for all t ∈ (0, T ).
Analytically characterizing these rigidity equivalence classes for an arbitrary graph has proved challenging. However, we are able to make a statement for the case of graphs that contains 3-simplex (a complete tetrahedron) as a subgraph. Fig. 4: The values of the Worst Case Rigidity index (17) as the position of node E changes linearly between the left and right configurations. Without edge CE (shown in yellow) the worst case rigidity index goes to 0 when E is co-planar with DAB, while with the yellow edge, the worst case rigidity index remains greater than 1.
For each complete tetrahedron, we define its orientation as the sign of the signed volume which is computed as
These preliminaries allow us to make the following statement:
Theorem 1: Let F 1 = (G, X 1 ) and F 2 = (G, X 2 ) be two minimally rigid frameworks in R 3 . If there exists a subgraph of G that is a 3-simplex and F 1 and F 2 contain the simplex with opposite orientation, the two frameworks lie in different equivalence classes.
Proof: Finding a smooth path x(t) for the vertices of a simplex from one orientation to the other requires the signed volume to smoothly change signs, passing a configuration where V = 0. When V = 0 for a simplex, one of the edges of the simplex is a linear combination of the others, meaning there is a redundant edge in the R matrix. For a minimally rigid graph the R matrix has 3n -6 rows, so any linearlydependent edges indicate the matrix is not maximal rank and hence not infinitesimally rigid.
One general question in the design of linear actuator robots is if an over-constrained network is necessary, or if a minimally rigid network is sufficient. We give an example where an overconstrained robot can achieve motion through a configuration that would represent a singularity were the robot minimally rigid (shown in Fig. 4). In this example, we first consider the robot to be only composed of the blue edges (edge CE is not present). In this case both the left and right configurations are infinitesimally rigid and simplex ABED has different orientation in each configuration, meaning that the two configurations lie in different rigidity equivalence classes by theorem 1. If the node positions are linearly interpolated between the two configurations, the rigidity index in (17) goes to 0 when node E is coplanar with nodes ABD. The addition of the yellow edge, which makes the robot over-constrained, allows rigidity to be maintained throughout the transition, as shown by the plot in Fig. 4. We note that with the yellow edge this graph is the fully connected 5-node graph, known as the K5 graph. The K5 graph displays another interesting property:
Theorem 2: The rigidity matrix R(x) for a robot represented by a complete graph of 5 or more nodes only loses rank at configurations where the robot has actuators in collision.
Proof: For a node in a complete graph to have an unconstrained infinitesimal motion, its neighboring edges must not span R 3 , meaning that all nodes must lie in the plane. Complete graphs with 5 or more nodes do not have planar non-crossing embeddings.
This result means enforcing the constraint that no actuators collide for the K5 graph naturally enforces the graph rigidity constraint. Whenever we evaluate a K5 robot in this paper, we leverage this result and do not enforce the rigidity maintenance constraint.
Our approach to finding a trajectory for a linear actuator robot is to use an optimization to solve for a discretized trajectory containing N conf ig configurations we denote as x j where j = 1, 2, ...N conf ig . The optimization solution guarantees that the configurations x j satisfy the constraints defined above which we now succinctly express as f (x j ) ≤ 0. However, the nonconvex nature of the constraints means that it is possible that the intermediate configurations (the configurations between x j and x j+1 ) may violate the constraints. To address this, we enforce a constraint that two sequential configurations must be close together in terms of the distance each node travels. We define this constraint as
We assume that the intermediate configurations between x j i and x j-1 i are given by linear interpolation. From work on sampling-based motion planning [52], the maximum violation of a constraint between two configurations can be bounded by using the Lipschitz constant, K of the constraint as follows
Given the Lipschitz constant for each constraint function, it is possible to augment the constraints with a buffer such that satisfying the buffered constraints and the constraint in 19 ensures satisfaction of the true constraint. In our case, we assume that the constraints already include this buffer. In practice we choose d move to ensure that two edges can not jump over each other without violating the collision constraint by picking 2d move ≤ d min .
Our first approach to solving the locomotion problem involves solving an online optimization to move the center of mass to a desired position for one time step. It acts greedily to minimize an objective function for a time step, and does not account for making and breaking contact with the ground. Our second method, presented in Sec. V and referred to as a two-tiered planning approach, extends this single step computation to an optimization over multiple steps. The twotiered approach directly accounts for the rolling behavior in the computation, but imposes restrictions that the robot must satisfy certain symmetry requirements.
The position of the center of mass is defined in terms of the mass matrix of the system, M ∈ R 3×3n . Without loss of generality, the quasistatic model allows us to assume that all mass is concentrated at the nodes of the system. In our case, we assume that all actuators are of uniform, evenly distributed mass, and thus half of the mass is assigned to each end of the actuator. The position of the center of mass is given by
where m vec,i is the sum of all of the partial masses assigned to node i. In the uniformly distributed case, m vec,i = di 2N L . With this mass matrix, we can express the velocity of the center of mass as a function of the actuator velocities
We can now pick any L that achieves a desired motion of the center of mass. The maximum rank of M is d, so for a system with many vertices M H -1 will have more columns than rows, and there is freedom in which ẋ is selected to move the center of mass. We define an optimization problem to pick a value of L that minimizes an objective function.
The kinematic relationships derived in the Sec. II apply to a continuous time system. To optimize the trajectory we work in discrete time, denoting the configuration of the robot with the superscript x j . In practice we determine the node velocities by linearly interpolate from the current configuration to the configuration that is the result of the optimization, and determine the necessary actuator velocities using the kinematic relationships. The optimization procedure takes as an input the configuration at x j-1 and optimize the next desired configuration x j . We seek to find a trajectory that maximizes some cost function J(x) while satisfying constraints. As this optimization minimizes the cost function over only one step, we refer to this optimization approach as the greedy method. The complete optimization problem is given as follows:
subject to
The choice of cost function J(x) will be discussed in the following section. Eq. ( 24) fixes the contact points and is the discrete time version of the ground constraint, where b is a vector of locations of the vertices in the support polygon. In the locomotion optimization we also enforce the linear constraint that no nodes pass through the ground, Gx > 0, where G = I n ⊗ diag([0 0 1]). We denote all of the feasibility constraints, including maximum and minimum actuator length (11), actuator collision constraints (12), angle constraints (14), and singularity avoidance constraints (17) as f (x j ) ≤ 0 as given in (26).
By defining this problem as an optimization problem, the system will take the action that instantaneously optimizes some objective, J(x). One intuitive choice for the cost function is J(x) = L(x) 2 = R(x) ẋ 2 , which penalizes large actuator velocities. In discrete time, we approximate this cost as:
As mentioned previously, one potential issue with a singlestep method is that persistent feasibility is not guaranteed. One heuristic to prevent the robot from getting tangled up in an unfavorable configuration is to try and keep the network as close as possible to a fixed operating point, such as attempting to keep all actuators close to a nominal length l N . This can be encoded with an objective function
We will quantitatively compare the results of using both of these cost functions in Section VI.
We find a feasible solution to the optimization problem using the sequential quadratic programming algorithm available in the matlab fmincon toolbox. The most computationally expensive part of this algorithm is repeatedly checking to see if the nonlinear constraints are violated, a process that could be parallelized in a future implementation. To speed computation, ∂f (x) ∂x is computed analytically before operation. We demonstrate the character of the solutions that result from this optimization we will show the types of trajectories generated when it is applied to different robots in the following sections.
1) Randomly Generated Robot: To show the generality of the algorithm to a wide variety of robots, the locomotion of a randomly generated minimally rigid 7 node robot is shown in Fig. 5. The initial configuration of the robot is obtained by starting with a triangular base and iteratively adding one node and connecting it with 3 randomly selected existing nodes. For each node, the node position is randomly regenerated until all constraints are satisfied. The objective function presented in ( 28) is used. For these simulations (and for simulations throughout this paper) the actuator lengths were constrained to remain between 0.5 and 4 units, the minimum angle between connected actuators was 10 degrees, and the minimum distance between actuators was set at 0.15 units. The minimum value of the worst case rigidity index was set at 0.005. The robot has an emergent, almost amoeba-like gait as it moves. Videos of this motion are available in the supplementary materials.
2) Trajectory Tracking: In order to demonstrate the ability of the system to follow a trajectory the K5 robot was controlled to move its center of mass towards waypoints that make up the corners of a predefined trajectory. The resulting trajectories when both ( 28) and ( 29) are used as the objective are shown in Fig. 6. We use the same values for the physical constraints as for the random robot, but do not enforce the rigidity maintenance constraint for the K5 robot due to Theorem 2. The variance from the prescribed trajectory occurs because of the rolling motion when the center of mass leaves the support polygon. In order to illustrate the effectiveness of this method in preventing constraints from being violated, Fig. 7 shows that during the trajectories shown in Fig. 6, the various physical constraints on the robot are often active but are not violated.
This trajectory tracking test also gives a sense of the robustness of the control algorithm. A downside to the approach of repeatedly solving the optimization is that persistent feasibility is not guaranteed, meaning it is possible that the network reaches a configuration where it cannot continue without violating some constraint. In the case where the objective function was (28), a configuration was reached where the device could not continue to match the desired center of mass motion without violating constraints (the blue trajectory in Fig. ). With the objective presented in (29), the planner finds a feasible path that completes the trajectory (the red trajectory in Fig. 6). Simulation results on a variety of trajectories show that a common failure mode of the system is if the robot rolls onto a very large support polygon, it may not have the ability to extend its center of mass and roll again. This failure mode may become less significant if future work included a frictional model of the ground and allowed the support nodes to slide along the ground. Interestingly, relaxing constraints does not necessarily guarantee the robot will be able to travel further before reaching a configuration with no feasible solution. Often, relaxed constraints such as a higher upper limit on actuator length lead to failure sooner, as the robot tends to reach more jumbled configurations early in the trajectory. Computation for completing the "S" trajectory involved solving the one-step optimization 1893 times, which took approximately 120 seconds on a laptop computer (Intel Core i7 Processor, 4 cores, 2.80 GHz, 16GB RAM). The average time to solve each optimization was 63 ms, with a standard deviation of 10 ms and a maximum time of 153 ms.
Note that at each time step these methods instantaneously minimize an objective while a desired motion of the center of mass is obtained. The algorithm can be thought of as greedily trying to move the center of mass. However, motion is not optimal for the entirety of the trajectory. The algorithm does not explicitly take into account making and breaking of contact with the surface, which would be required to discuss the optimality of an entire trajectory. To enable discussion of optimality over several steps as well as to directly consider the rolling behavior of the robot, we extend this method to a tiered planning approach.
In this section we extend the one-step optimization of the previous section to an optimization over many configurations of the robot. This multi-step optimization directly accounts for the rolling behavior, whereas the previous method moved the center of mass without consideration for the rolling motion. We use an offline optimization to compute trajectories from a predefined configuration centered on one support polygon to the same predefined configuration centered at the next support polygon, but with the node correspondences changed. This precomputed trajectory serves as a motion primitive for a high level planner that computes a series of support polygons that lead from the robot's initial position to a goal region. When deviations occur from the preplanned trajectory, the only Fig. 6: The path of the center of mass as it travels to each waypoint of an "S". Note that when using the objective in (28) the LAR reaches a point from which it cannot continue. With (29) as the objective the robot completes the trajectory.
computation that occurs online is using the high level planner to adjust the path of support polygons. As the final path of the robot is composed entirely of feasible motion primitives, the resulting path is guaranteed to be feasible. In this section we discuss the necessary symmetry requirements for such a trajectory to exist, present the optimization setup to solve for the motion primitive and our use of a high level planner to combine the motion primitives. We then discuss methods of smoothing the motion primitives during trajectories over a series of support polygons.
We now detail the symmetry requirements that allow a motion primitive to be optimized offline and then stitched together online into long trajectories. This means that the robot must finish a motion primitive in a configuration identical to the starting configuration, but with different node correspondences. If the robot begins in an arbitrary configuration, a path between the initial configuration and the symmetric configuration must be computed and executed. For the symmetric configuration, we restrict the shape of the support polygon of the robot to be an equilateral triangle, meaning that the robot's motion will be over a grid of equilateral triangles. The symmetry requirements are illustrated in Fig. 8. To enable the same primitive to be reused repeatedly, we constrain the starting configuration to have mirror symmetry about the three lines that originate at the vertices of the support polygon and bisect the opposite edge of the triangle, as shown by the red Fig. 9: The resulting motions obtained from the optimization given in (30). The top row shows a rolling gait of the K5 graph, the middle row an everting gait of the K5 graph, and the bottom row shows a rolling gait for an octahedral robot.
subject to All Configurations:
) Non-Tipping Configurations:
The objective function ( 30) is a multi-step extension of (28), where L(x j ) is the vector of all of the edge lengths of the graph with node locations given by x j . This objective penalizes sudden and large changes in the lengths of the actuators, and hence favors trajectories that require small changes in actuator lengths. The linear equality constraint in (31) constrains the location of the contact points for each configuration. Note that C 1 = C k , ∀k ≤ j * , and
as only two support polygons are used throughout the optimization. In (32) three linear inequality constraints keep x com within the support polygon at each configuration to prevent premature rolling. The variable g j i and h j i describe the parameters of the line along edge i of the support triangle. For each configuration, we denote f (x j ) ≤ 0 to represent all of the physical feasibility constraints for one configuration. We write the constraint x Nsteps = x f to ensure the proper final configuration.
We note that for the transition between the tipping configuration and the next configuration (x j * and x j * +1 ) large motions of the node positions are possible due to the rolling, even though change in the edge lengths may be small. For the rolling step alone, we constrain the change in edge lengths to be below a fixed threshold as shown in (37) to prevent the robot from jumping over a physical constraint, such as an actuator collision constraint, between configurations.
1) Tipping Constraints: In addition to the constraints that ensure that each configuration is feasible, we also impose additional constraints that ensure that the robot tips at the predefined tipping configuration. The linear equality constraint in (34) ensures that at the tipping configuration the center of mass lie on the tipping edge of the support polygon. We denote the new node in the support polygon after the tip as node s. The two quadratic equality constraints in (35) and (36) ensure that the positions of node s, denoted x j * s is the proper distance away from each node on the rolling edge, denoted x j * c1 and x j * c2 . Note that these constraints do not fix the height from which the robot tips onto the next support polygon. Were the height and position of the next point specified exactly, the two quadratic constraints would be replaced by three linear constraints, but the optimization would lose the ability to change the tipping height.
2) Optimization Results: The motion primitives produced by solving the optimization problem are shown in Fig. 9. For these results, we use N steps = 40, with the tipping configuration j * = 20. We initialize the optimization such that every configuration or after the tipping configuration is exactly the starting or ending configuration respectively. We initialize the tipping configuration with, all nodes of the initial and final support polygon in place, and the nodes that are not part of the support polygon positioned such that the center of mass is on the tipping edge. The optimization was solved using the fmincon solver available with Matlab.
The top two rows of Fig. 9 correspond to the two different node correspondences for the K5 graph discussed previously. The first row is the resulting motion primitives when using the correspondence in Fig. 8B. Here the center node of the robot before rolling remains the center node after rolling. The second row is the resulting motion when the correspondence in Fig. 8C is used. In this case the node initially in the center of the robot becomes the new node in the support polygon, and the top node of the robot remains the same both before and after rolling. We will refer to the gait with a constant internal node as the rolling gait, and the gait with the switching center node as an everting gait, as the robot seems to be everting its inside and outside as it moves. From a practical perspective, the fact that the top node of the everting gait always remains off the ground could allow it to house cameras or other components. We note that this everting gait requires an over-constrained network, as it requires that a simplex present in the initial graph switch its orientation as demonstrated in Fig. 4. The resulting motion primitive for the octahedron is shown in the third row of Fig. 9. The computation time for these offline primitives was 127, 134, and 166 seconds for the K5 inverting gait, K5 rolling gait, and octahedron gait respectively.
We can also compare the resulting motions in terms of cost. As computed by (30), the optimized motion primitive for the everting K5 has a cost of 0.198, the rolling K5 primitive has a cost of 0.062, and the octahedron has a cost of 0.025. Another interesting comparison between the motion primitives is the ratio of the maximum and minimum actuator length. The ratio of the overall longest actuator to the overall shortest actuator is 3.11 for the everting gait, 2.85 for the rolling gait, and 1.58 for the octahedron. These results demonstrate the need for high elongation actuators.
Given a motion primitive developed by the optimization, we need a planner to specify a series of support polygons from the initial configuration to the goal. This task corresponds to planning a path on a triangular grid of candidate support polygons. We use the A* algorithm for this task, but note that any discrete planning algorithm could be used for this task. An example of A* finding a path through obstacles is shown in Fig. 10. If a feasible trajectory is found that solves the optimization and a feasible path is found between the starting configuration and the goal region, a feasible path that satisfies all constraints is possible from the start to the goal region. Note that in the case of an environment with obstacles, the collision checking performed as part of the A* algorithm depends in part on the robot gait. The maximum extent of the computational gait must be used by the planner to ensure that it is possible to move from one support polygon to another. If a path of collision free support polygons that leads from the start to the goal exists, the A* algorithm is guaranteed to find it. However, it is possible that if A* fails to find a path, the robot could pass through the environment by using a different motion primitive.
In this section we leverage symmetry in the robot and the triangular grid of candidate support polygons to consider motion primitives for moving between several support polygons, as opposed to just moving from one support polygon to its neighbor. We present two approaches: one where we relax the requirement to return to the symmetric configuration between every step to a requirement to return to the configuration at larger numbers of intermediate steps, and a second approach where we optimize a trajectory that enables a robot to continue in a straight line indefinitely without returning to the symmetric configuration.
1) Smoothing over Multiple Support Polygons: In the extreme, we could optimize directly over the entire trajectory from beginning to end, but such a procedure may be expensive to compute online. Instead, we quantify the marginal gain of optimizing over trajectories of increasing length, but while maintaining the same support polygons. We note that for a robot traveling through a triangular grid, if we eliminate the option to move backwards at every step the robot can choose to roll over the left or right edge. Shown in Fig. 11 is a partial triangular grid that gives the sequence of turns to arrive at each cell, assuming the initial motion is from the "start" to the "1" cell. Each path can be represented by a p -2 digit binary word, where p is the number of transitions between support Fig. 11: The number of possible trajectories available when the robot rolls over one edge. Note that the trajectories denoted with a prime are the same as their counterparts, but with each letter switched. The marginal gain of adding more steps seems to be decreasing. polygons or rolling events. By symmetry of the robot and the grid of support polygons, switching all entries in the binary word results in a mirrored trajectory. This means that for p (where p ≥ 2) steps there are 2 p-2 possible paths to compute. This means that for paths with two rolling events there is only a single motion primitive possible, meaning there is no loss of generality for optimizing the trajectory over two steps as opposed to a single step.
To understand the cost savings of optimizing over multiple support polygons, we compute the cost of moving 1 to 4 steps along the pattern of support polygons shown in Fig. 12, along with the direct comparison of the center of mass path when both 1 and 4 steps were used. The cost to complete a single roll is shown in Fig. 13. We note initial improvement in the cost when moving from one step to two steps, but observe diminishing returns by using longer and longer primitives. Qualitatively, the motion of the center of mass in the smoothed and unsmoothed trajectories is shown in Fig. 12. For the K5 graph a reasonable compromise appears to be to always use the two step motion primitives unless the robot is within one step of the goal. We illustrate combined behavior of the A* planner and the smoothed primitives to navigate between the waypoints of the "S" trajectory shown in Fig. 14.
2) Gaits with no return to the nominal configuration: The smoothing methods presented previously in this section relaxed the requirement of returning to the symmetric configuration from every step to every N steps, where N is some integer number of steps. An equivalent optimization approach could also be used to develop primitives that enable moving between different intermediate configurations without passing through the symmetric configuration. A key question with this approach is how to define the best intermediate configurations. One option is to include the shape of the intermediate configuration as part of the optimization itself. As a demonstration, we develop primitives that allow the octahedron and K5 robots to locomote along an arbitrarily long straight path of support polygons, similar to the motion shown in Fig. 12. Whereas we previously optimized a trajectory that starting at a given symmetric configuration and ending at equivalent symmetric configuration, we now optimize a trajectory that starts at a tipping configuration and ends at an equivalent tipping configuration for the next support polygon, where the shape of the tipping configuration itself is part of the optimization. We encode the symmetry between the first and last configurations as follows
Where A and b define a linear transform and necessary assignment of node correspondences to ensure that the initial and final configurations are equivalent. We repeat the optimization presented in ( 30)-( 38), but including the initial configuration as one of the optimization variables, and replacing (38) with (39). The resulting gaits are demonstrated in the supplementary video. The cost of this smoothed gait for the octahedron, K5 inverting gait, and K5 rolling gait is 76%, 51% and 46% the cost of the repeatedly using the one step trajectory that starts and ends in the symmetric configuration, representing a substantial savings. Utilizing these gaits online in the robot requires storing the repeating gait as well as the trajectory to move to and from the symmetric configuration to be used at the beginning and end of the straight line trajectory. This means that the memory required to store this gait is equivalent to the memory required to store a two-roll primitive. Due to the cost to move from the initial configuration to the rolling configuration, the multi-step primitives such as those shown in Fig. 12 are superior for short sequences of support polygons. However, as the length of the trajectory increases, the cost of using the repeated gaits approaches the cost obtained by optimizing over the entire trajectory, but requires a smaller amount of memory to store. Future work could seek to define intermediate gaits and other shapes that enable different behaviors such as turning. Fig. 12: A series of one step trajectories stitched together (red) compared with the smoothed optimization over three steps (black) for the rolling gait of the K5 network. Note that the path of the center of mass is more direct for the smoothed primitive than the compilation of single step primitives.
We now compare the behaviors of the greedy, roll-unaware planning method presented in Section IV with the two-tiered planning method presented in V. We find that, on average, the two-tiered planning method finds more efficient trajectories than the greedy approach. In addition, the two-tiered planning approach always finds a successful trajectory if a sequence of support polygons exists that leads to the goal, while the greedy approach is often unable to find a successful trajectory.
However, we note that the one-step planning method applies to every infinitesimally rigid robot, while the two-tiered planning approach applies only to robots of a restrictive symmetry class.
Conceptually, we can compare the behavior of the two planners by comparing the resulting center of mass trajectories in Fig. 6 and Fig. 14. With the two tiered planning approach, the trajectory of the center of mass takes a less direct path between waypoints, as the constraint to move the support polygon along the triangular grid ensures that the center of mass does not move in a straight line. Despite the apparent inefficiency of the trajectory from the two-tiered planner, we find that it results in lower cost trajectories. We hypothesize that this occurs because the robot remains in a better conditioned state. The two-tiered planner generates trajectories with consistent motion of all of the free nodes, while in the trajectories of the greedy planner free nodes seem to be flailing about a relatively steady center of mass trajectory.
For a quantitative comparison of the performance of the planners, we generated 100 sequences of 5 random waypoints in a 5 unit by 5 unit region and use the proposed planning methods to find a trajectory to visit the waypoints sequentially. As the output of the optimization is a kinematic trajectory, we scale the trajectories such that completing the entire trajectory takes 1 unit of time, and convert the trajectory to continuous time by linearly interpolating the node positions between the discrete configurations returned by the optimization. This rescaling ensure an equivalent average velocity between the different experiments, and allows a direct comparison in terms of the cost. To evaluate the cost of the trajectories we use the following cost function:
where d is the sum of the straight line distances between the waypoints. This cost is a continuous time version of (28), divided by the path length to give an efficiency metric as the average cost to move a unit distance. Using both the K5 and the octahedron robot, Fig. 15 compares the efficiency of the paths resulting from the two-tiered planning method (using both the rolling and everting primitive for the K5 graph), and the greedy method using both ( 28) and ( 29) as the objective.
For both the octahedron and the K5 graph, the two-tiered planning approach attains lower cost and lower variance than the greedy approach for the same robot. Interestingly, for both the octahedron and the K5 graph, a lower overall cost is obtained by using (29), the cost function that penalizes actuator for deviating from a nominal length, as the cost as opposed to (28), which penalizes changes in actuator length at each time step and is the single-step version of (40). This seems to indicate that long term efficiency is achieved by keeping the robot in a relatively well-conditioned state.
In addition to cost, the other key criteria by which to evaluate the planners is their ability to find a complete path without violating constraints. For the 100 randomly generated trajectories and using the K5 robot, the one-step optimization method successfully found a path between all waypoints for 12% of the trials when using (28) as the objective, and for 70% of the trials when using the formation-control based objective given in (29). We repeated the experiments with Average Error (cm) Octahedron K5 Rolling K5 Everting Fig. 16: A comparison of the average error between the planned quasistatic trajectory and the position of the nodes in a dynamic simulation where the quasistatic input serves as an input to a PID controller. The error is compared against the overall time that it takes to complete the trajectory. The average error decreases as the movement proceeds more slowly.
overall speed of the robot decreases. These demonstrations also illustrate that the quasistatic trajectories lead to useful behaviors in a fully dynamic system.
In this paper we have derived the differential kinematics for networks of linear actuators connected at universal joints and have shown that if the embedded graph describing the robot is infinitesimally rigid, any desired motion of the nodes can be achieved through some motion of the edges. We then frame the locomotion problem as a nonlinear optimization over the node positions, while enforcing constraints that guarantee the feasibility of the robot. We also discuss that constraints to maintain the infinitesimal rigidity of the robot tend to divide the state space of the robot into separated regions, even though the singular configurations themselves make up a set of zero measure. We discuss the control of both minimally rigid graphs and over-constrained graphs, and demonstrate that over-constrained graphs can achieve some behaviors that minimally rigid graphs cannot, such as the everting locomotion gait of the K5 graph. We present two planning schemes: one where we solve a single step nonlinear optimization online to achieve a desired instantaneous motion of the center of mass, and another where we optimize over many configurations that compose a motion primitive, including in the optimization direct consideration of the rolling behavior. The single-step approach is applicable to robots of arbitrary configuration, but there is the possibility that the robot will reach a state from which it cannot continue in the desired direction without violating physical constraints. While these is no guarantee of persistent feasibility with this approach, we have found that long trajectories can be achieved based on the choice of the cost function. The two-tiered approach ensures persistent feasibility, but requires the robot to satisfy certain symmetry properties.
In future work we will attempt to blend the properties of these two control approaches, namely by finding ways to guarantee persistent feasibility for robots composed of linear actuators in arbitrary infinitesimally rigid configurations. We will also extend our planning approach to directly consider dynamic effects, explicitly considering forces in the members and incorporating the inertia properties of the system. One possibility is to use these kinematic trajectories as starting points for a dynamic model, similar to the methods in [31]. We will also explore taking the centralized controllers presented in this paper and finding a distributed version of a similar controller, where computation is performed locally at the actuators in the system.
In future work we will demonstrate this system with novel robotic hardware of the type presented in [17]. We have noted that the current optimization procedure does not have any guarantees on global optimality. We will work to find convex relaxations of the non-convex constraint such that we are able to find suboptimality guarantees for certain frameworks.