A robot made of robots: Emergent transport and control of a smarticle ensemble

See allHide authors and affiliations

Science Robotics  18 Sep 2019:
Vol. 4, Issue 34, eaax4316
DOI: 10.1126/scirobotics.aax4316


Robot locomotion is typically generated by coordinated integration of single-purpose components, like actuators, sensors, body segments, and limbs. We posit that certain future robots could self-propel using systems in which a delineation of components and their interactions is not so clear, becoming robust and flexible entities composed of functional components that are redundant and generic and can interact stochastically. Control of such a collective becomes a challenge because synthesis techniques typically assume known input-output relationships. To discover principles by which such future robots can be built and controlled, we study a model robophysical system: planar ensembles of periodically deforming smart, active particles—smarticles. When enclosed, these individually immotile robots could collectively diffuse via stochastic mechanical interactions. We show experimentally and theoretically that directed drift of such a supersmarticle could be achieved via inactivation of individual smarticles and used this phenomenon to generate endogenous phototaxis. By numerically modeling the relationship between smarticle activity and transport, we elucidated the role of smarticle deactivation on supersmarticle dynamics from little data—a single experimental trial. From this mapping, we demonstrate that the supersmarticle could be exogenously steered anywhere in the plane, expanding supersmarticle capabilities while simultaneously enabling decentralized closed-loop control. We suggest that the smarticle model system may aid discovery of principles by which a class of future “stochastic” robots can rely on collective internal mechanical interactions to perform tasks.


Self-propulsion (1) is a feature of living and artificial systems across scales—from crawling cells to swimming spermatozoa (2), micro- (3) and nanoswimmers (4), running cockroaches (5), and robots (6, 7). It is generally assumed that self-propelling systems require carefully orchestrated integration of many diverse components to perform the seemingly simple behavior of spatial translation. Thus, artificial locomoting systems typically consist of a central controller, a set of actuators and sensors to perform feedback control, and an objective function written in terms of individual system states; such designs have led to progress in machines that robustly and nearly autonomously roll (8), fly (9), and walk (10, 11) in relatively predictable environments.

In contrast to such “deterministically” designed robots, future more “stochastically” designed robots could generate self-propulsion using systems in which a delineation of components is not so clear, such that many redundant and generic elements fluidly interact and collaborate to achieve complex tasks (Fig. 1). Although such designs are potentially advantageous due to wide system reconfigurability and robustness to component damage, it is not yet clear how to build such a system to operate in natural environments. There are several reasons for this, some of which have been anticipated by insights from modular and swarm robotics (1215), physics of active matter (13, 1621), amorphous computing (22), and engineering of reliable systems from unreliable components (23).

Fig. 1 Stochastic robotic collectives.

Future robots may be composed of components whose delineation is neither clear nor deterministic, yet are capable of self-propulsion via the expression of ensemble-level behaviors leading to collective locomotion. In such a robot, groups of largely generic agents may be able to achieve complex goals, as routinely observed in biological collectives.

For one, future stochastically designed robots (and collectives/swarms) may contain so many components (members) (24, 25) that it might be infeasible to carefully arrange and couple the elements to generate coordinated translation or rotation. Further, like in crawling cells where locomotion is generated through cytoskeletal reconfigurations via shape-changing proteins, individual elements may be task-incapable [e.g., unable to move on their own, unlike in collective robot locomotion via mechanical rectification of individual bristlebots in (26, 27)]. In such situations, the robot’s objective should not depend on deterministic interactions between components but instead on emergent ensemble-level behaviors (25, 28). Thus, it becomes a question of leveraging or mitigating the inherent uncertainty of internal component interactions to develop reliable control schemes of the ensemble.

Traditional control synthesis techniques determine which inputs best, and most robustly, enable a system to achieve an objective, such as self-propulsion. It can be challenging enough to find control inputs that realize a well-defined objective in a deterministic system. In the case of a robot composed of robots with highly complex interactions between the system and the environment, and no single configuration of individual components necessary for the robot to achieve locomotion, control synthesis using traditional methods is infeasible. The notion that an ensemble may be able to accomplish a goal independently of the specification of its individual states is incompatible with typical theories of control that assume a central control architecture with full state information.

Here, we sought to discover principles by which a collective can overcome individual locomotor limitations via opportunistic but stochastic mechanical interactions among individuals. Specifically, we studied a simplified robophysical (29) model of controllable, smart, active particles—smarticles—that are immotile but have mutable shape. An enclosed smarticle ensemble—a supersmarticle—however, can self-propel diffusively using interactions arising from the shape modulation of smarticles. Despite stochastic interactions between elements, a supersmarticle is capable of directed motion by selectively inactivating its constituents, which we demonstrate by achieving endogenously steered phototaxis. To understand the supersmarticle diffusion and its dependence on internal mechanical interactions, we developed a model based on kinetic theory. To further explore the ensemble’s abilities, we introduce a data-driven algorithm that enables decentralized control synthesis with respect to ensemble properties. Using this algorithm, we modeled supersmarticle dynamics and demonstrate that the ensemble is capable of rich locomotion by taking advantage of more complex control strategies. We validated our algorithm by leading the supersmarticle on a simple path, demonstrating that collective locomotors may be reliably controlled through their ensemble properties despite being composed of stochastically interacting unreliable elements.


Smarticle dynamics

The smarticle’s form (Figs. 1 and 2A and Materials and Methods) was inspired by insights from a previous study of rigid, non-active “u-particles” (30), which demonstrated how material properties of an entangling collective could vary via changes in the shape of its constituents. Smarticles, when active, perform a “square gait,” inspired by the dynamics of Purcell’s three-link swimmer (31, 32), and depicted in Fig. 2B. Outside of a frictional medium [e.g., (31, 32)], when resting in an orientation where the links’ axis of rotation is parallel to the normal of the surface it rests on, smarticles are incapable of translation or rotation (Fig. 2C) over hundreds of oscillation cycles. The moving links rest above the central link and never interact with the surface (movie S1).

Fig. 2 Smarticle robot dynamics.

(A) Top view schematic w = 5.3 cm and l = 4.9 cm. (B) Clockwise (CW) square gait, with key configurations enumerated. (C) Drift of a single smarticle on a flat surface, executing a square gait over 38τ. (D) Tracked trajectory of a smarticle within an ensemble of other self-deforming smarticles; color gradient (blue to red) represents passage of time 47τ, with τ = 1.6s.

Despite their inability to significantly self-propel, an individual robotic smarticle’s position and orientation can change as a result of a collision, as shown in Fig. 2D. When viewed as an ensemble, a “cloud” of self-deforming smarticles may display weak cohesion on short time scales, forming a rudimentary collective flocking unit (Fig. 3, A and B, and movie S2). That is, unlike single smarticle experiments, we found that the center of mass (CoM) of the cloud could diffuse over scales comparable with the size of a smarticle (see fig. S1).

Fig. 3 Smarticle cloud dynamics.

(A) Snapshot of experimental trial, with the dashed line indicating the boundary of the convex hull area AC. The cloud’s CoM trajectory is illustrated in red, beginning at the black dot and ending at the red dot. Experiment ran for 113τ. (B) Center link trajectory of geometrically repulsive (top) and attractive interactions (bottom). (C) Evolution of φ averaged over 20 trials (black, with gray shaded region representing a single standard deviation); four individual trials are shown in blue, red, green, and brown lines. (D) 〈V2〉 averaged over 20 cloud trials. Raw data are in black; the blue line is moving mean with a window size of 1τ. Red line and area surrounding it represent mean value and single standard deviation of 〈V2〉 noise of an experiment lasting 10τ with seven moving, but non-interacting, smarticles. Here, gait period τ = 1.6 s.

Because of interactions between smarticles, the area fraction φ typically decreased over time as in Fig. 3C. Here, φ = nAp/Ac, where Ac is the area of convex hull of the smarticles (bodies and arms) in the cloud (Fig. 3A), n is the number of smarticles in the system, and Ap is the area of a single smarticle. The decrease in φ was not always monotonic; in certain trials, increases in φ occurred (Fig. 3C). Despite purely repulsive interactions at surfaces, smarticles could both repel and attract their neighbors (see Fig. 3B). This emerges from the particle geometry: Collisions between particles in concave configurations can generate attraction via arm entanglement (30).

After sufficient time, the cloud’s mobility slowed as smarticles separated and no longer interacted strongly. We quantified collective mobility using the cloud’s “granular temperature,” defined as V2=1/3v2nvn2N, where v=ẋ2+ẏ2+(2l+w)θ̇2 sums the translational plus rotational velocity of n smarticles of length l and width w and averages over N experiments (33, 34). On long time scales, a single experiment’s V2 may approach the noise floor (seen in Fig. 3D) (35), thereby limiting the flocking ability. For this study, we determined the noise floor empirically by measuring the granular temperature of non-interacting smarticles.

Supersmarticle dynamics

Given the correlation between φ (Fig. 3C) and 〈V2〉 (Fig. 3D), we hypothesized that we could sustain locomotion on longer time scales by constraining φ of the collective. To achieve this, we confined five smarticles within a ring, creating what we call a supersmarticle. Each smarticle in the supersmarticle starts at a random phase in the square gait and continuously performs a square gait inside an unanchored rigid ring of radius R = 9.6 cm and variable mass m ∈ [9.8 g, 207 g] (Fig. 4A). It takes t = 225τ (where τ = 1.6 s) before two smarticles are >π out of phase. The ring diameter was chosen such that φ and 〈V2〉 remained high, yet there was enough area that jamming was rare and self-resolvable.

Fig. 4 Collective confined diffusion.

(A) Supersmarticle top view; ring inner radius is 9.6 cm. The four gray spheres were used to track the motion of the ring. (B) Granular temperature of five active smarticles confined in a ring; black line is raw data over 10 trials, and blue is a moving window mean with a window size of 1τ. (C) Trajectories, from an experiment, of a smarticle inside the ring (purple), and the ring’s center of geometry (blue). (D) Experimental tracks of ring trajectory for 50 trials; mring = 68 g. The black circle represents the size and initial position of the ring. (E) MSD averaged over 50 and 80 trials, for the active and inactive systems, respectively, all lasting 75τ. The inset shows the average change of γ for active (black) and inactive (blue) systems. The oscillation seen in both the MSD and γ is related to the gait period τ (where τ = 1.6 s).

The ring confinement maintained φ at approximately the value observed at the initiation of the cloud trials (see Fig. 3A). Similarly, 〈V2〉 of the supersmarticle system (Fig. 4B) remained at approximately the value found at the highest φ in cloud trials (Fig. 3, C and D). This led to persistent diffusive transport of the supersmarticle (Fig. 4B). Within the ring, individual smarticles displayed complex interactions, often displacing an amount comparable with, or greater than, the displacement of the ring itself, as shown in Fig. 4C.

Tracking the supersmarticle’s motion for a ring of mass m = 68 g (movie S3) revealed no correlation between final angular position between trials (e.g., Fig. 4D). We used σ2(t), the mean square displacement (MSD) of the ring, to characterize the motion: σ2(t) = 〈x2(t)〉 − 〈x(t)〉2, where σ2(t) ∝ tγ and γ specifies the type of diffusion the system undergoes. The supersmarticle exhibited different types of diffusion—normal (0 < γ ≤ 1), superdiffusive (1 < γ < 2), and even approximately ballistic (γ ≥ 2)—depending on the time scale observed (36). The short time scale regime was consistent with γ = 1 (Fig. 4E), indicating normal diffusive motion. The long time scale regimes were best fit with γ ≈ 1.45 representing directionally invariant superdiffusive motion.

We found that, if a smarticle near the boundary maintained a fixed straight shape or became “inactive” (Fig. 5A), the supersmarticle displayed directed drift on short time scales (movie S4). Because the angular position of the inactive smarticle around the ring was not fixed, drift in a constant direction was not observed on longer time scales in the laboratory frame (see fig. S2). When trajectories were examined in the frame of the inactive smarticle (Fig. 5B), the bias in drift toward the inactive smarticle became clear. In Fig. 5C, the cumulative displacements are shown in the continuously rotating frame attached to the center link of the inactive smarticle such that S(t)=i=0tΔsiR̂i and S(t)=i=0tΔsiR̂i. Here, Δsi denotes the vector connecting the center of the ring at consecutive instants in time, and R̂i,R̂i are the unit vectors specifying the local frame (Fig. 5B). As with the fully active supersmarticle, the dynamics of the supersmarticle containing an inactive smarticle were superdiffusive and, at short time scales, approximately ballistic, as indicated by γ ≈ 2.

Fig. 5 Biasing supersmarticle transport.

(A) Supersmarticle schematic, with the inactive smarticle in red. (B) Supersmarticle trajectory frame transformation from laboratory to inactive smarticle frame. (C) Supersmarticle trajectories rotated into the laboratory frame, where axes are now the perpendicular and parallel components to the frame of the inactive particle.

Statistical model

To understand the supersmarticle diffusion and its dependence on internal mechanical interactions, we developed a model based on kinetic theory. Formally, the average displacement of the ring would be given by ΔR=P(Ω)ΔR(Ω)dΩ, where Ω represents the microstate (i.e., position, orientation, and heading) of the supersmarticle constituents immediately before a collision, and ΔR is the ring displacement due to an individual collision. The resulting mean ring displacement, ΔR, is then computed by integrating displacements due to individual collisions over the microstate probability distribution. However, because we do not have access to the detailed relation of Ω to the complicated smarticle-smarticle and smarticle-ring collisions, this calculation is intractable, demanding the development of a simplified ensemble model.

We imagine active smarticles rattling inside the ring and colliding against the ring and the inactive smarticle. The role of the active smarticle in the supersmarticle is simplified as simple contacts around the ring. The contacts are abstracted as nudges (Fig. 6A). Each nudge has a uniform probability to act in any direction. As a result of symmetries in the system geometry, we may partition the space of Ω into six distinct types of collisions. Six collisions arise from two independent factors: whether or not the inactive smarticle is in contact with the ring and which of three regions of the ring the active smarticle contacts (denoted by roman numerals in Fig. 6B). Each of these individual collision types generates a unique response on the collective system.

Fig. 6 Statistical model of supersmarticle transport.

(A) Schematic of the theoretical collision model. (B) Three regions with distinct collision types for the theory as described in the text. (C) Theoretical (red) and experimental (black and blue) data for velocity versus mass ratio M, showing mean and standard deviation. The blue data point is offset in M for visibility and represents an experiment where the inactive particle was endogenously chosen by light (see text) for 40 trials. (D) Distributions of drift speed probabilities for M regimes.

This simplification leads to a model with two random variables. The first is an angle Θ that represents the direction of an individual nudge and takes a value between 0 and 2π. The second is a binary variable β that represents whether or not the inactive robot is in contact with the ring. These variables together represent the various microstates Ω of the system. Depending on the value of the random variables, an individual nudge can move either the ring or both the ring and the inactive particle. Then, focusing on the movement of the ring, we must determine RA and RB, which are the distances that a nudge will displace the ring when moving only the ring, and when moving both ring and inactive smarticle, respectively. The possibilities for ring movement are summarized in table S1.

With RA and RB, we can describe the simplified model of the supersmarticle ensemble in expectation, which we decompose into parallel and perpendicular components (see Fig. 5B). Denoting the proportion of time that the inactive smarticle is in contact with the ring as λ, the frequency of nudges as f, the amount of time the supersmarticle has been moving as T, the inactive smarticle’s angular diameter as Ψ (see Fig. 6B), and treating each nudge as an independent event, the expected component of the velocity of the ring along R̂ of the inactive smarticle isv=(f/π)[λ(RARB)(1sin(Ψ))(1λ)RAsin(Ψ)](1)

The perpendicular component of the ring velocity is simplified substantially because we know by symmetry that ΔR=i=16ΔRi=ΔR̂. This is to say that the mean displacement of the ring averaged over all distinct collisions is in the parallel direction; hencev=0(2)

However, the variance along this direction is non-zero. The corresponding variances to the expected parallel and perpendicular velocities, Var[v] and Var[v], are detailed fully in Materials and Methods.

Last, to completely specify this model, we must calculate RA and RB. To this end, we determine the relationship between the mass of the ring and the distance it moves from a nudge by modeling the active smarticles as pistons pushing on a sliding mass (see Materials and Methods). The predictions resulting from this model are plotted in Fig. 6C and fig. S3.

The theory correctly predicts the supersmarticle’s drift speed relative to that calculated experimentally as 〈v〉 = S(T)/T with T = 75τ. The theory predicts that the direction of 〈v〉 will reverse for a large enough ring mass. Directionality depends on the mass ratio M = msmarticle/mring between the inactive smarticle and the ring, with reversal at a critical value of M ≈ 0.8. To test this prediction, we conducted experiments for a series of different ring masses. The results are summarized in Fig. 6C. The theory closely predicts the mean velocity, including direction reversal (movie S5). Although the theory predicts 〈v〉 = 0 (fig. S3), we observed slight discrepancies, particularly at larger M. We attribute these discrepancies in variance to correlations between collisions, whereas the theory assumes that collisions are independent.

The model elucidates the physics governing the dependence of 〈v〉 on M, as a function of the ensemble’s internal mechanical interactions. Consider first the high-M limit. The three collision types involving the (light) ring but not the (heavy) inactive smarticle dominate the net motion (see Fig. 6B). Both of the forward collisions (region 2) are of this type, as is one rearward (region 3) collision, resulting in a relatively large positive ΔR̂. Conversely, in the low-M limit, five of six collision types give rise to nearly equal magnitude ring displacements, the exception being the forward collision (region 1) of the active smarticle with the inactive smarticle when the latter is not in contact with the ring, in which case the ring displacement is exactly zero. This deficit in the forward-directed ring displacement results in a (small) negative value for the net displacement.

The close agreement between theory and experiment for the drift speed velocity is perhaps unexpected: With only N = 4 active smarticles, it is not clear that a purely statistical kinetic theory approach should work. The theory overestimates the observed fluctuations in Fig. 6C, an indication of substantial correlations in the smarticle swarm collisions, which will be the focus of future work. Yet, despite this incongruity with the variance, the derived theoretical model is still capable of generating the directed motion of the ensemble observed experimentally in the frame of the inactive smarticle.

Directing a phototaxing supersmarticle

On the basis of the intuition gained from the kinetic model, we programmed smarticles to inactivate when light detected from its photosensor exceeded a threshold. When illuminated at low angles (i.e., in the plane of the smarticle light sensors), photoinactivated smarticles occlude light from neighbors further from the source (Fig. 7A, inset), creating a situation similar to that analyzed in the previous section. The inactivated smarticle occludes the light from its neighbors: The straightening and resulting occlusion of light serves as a decentralized and stigmergic directive. The inactive smarticle is affecting the motion of the ring by affecting the motion of the remaining smarticles. This decentralized strategy has been used in previous swarm robotic collectives to generate group movement and transport without requiring explicit communication between agents (37, 38).

Fig. 7 Endogenous supersmarticle phototaxis.

(A) Trajectory from an experiment of a self-directed (endogenously forced) photophilic supersmarticle tracking a static light source (movie S6). Inset: Schematic showing how a smarticle in the straight configuration can occlude light from smarticle behind it. (B) Map depicting when and which smarticles endogenously inactivate. (C) Lorenz plots detailing general equality of smarticle inactivity over a 25-min endogenous trial consisting of five separate excursions in different directions (see movie S6). Over the complete trial, we found G = 0.21, as shown in the bolded line. The unbolded lines are the Lorenz curves for the five separate excursions, where we found G = [0.28, 0.4, 0.42, 0.34, 0.49].

However, we found that rather than regulate the angular location of an individual inactive smarticle, the static light source induced a switching sequence of inactive smarticles, leading to supersmarticle phototaxis. Because collisions in the ring can cause an inactive smarticle’s position to shift, when an inactive smarticle was dislocated from its lighted position, it switched to the active state. Consequently, an active smarticle could then be nudged into a position to receive enough light to become inactive. Hence, the supersmarticle phototactic drift was via endogenous steering, that is, where smarticle immobilization was spontaneously selected for without external feedback (see Fig. 7A and movie S6).

The endogenously forced system drifted in a preferred direction in the laboratory frame with a similar 〈v〉 to that of the non–light-driven system (Fig. 6C), whose drift was only observable in the frame of the inactive smarticle. This is remarkable given the complex switching dynamics of the inactive smarticle: For example, depending on distance and orientation relative to the light, it was possible for multiple smarticles to be simultaneously inactive as depicted in Fig. 7B. Moreover, the rotational symmetry of the supersmarticle allows one to infer that if the supersmarticle can translate in one direction, it should be able to translate in another direction by selecting different inactive smarticles.

To further highlight the supersmarticle systems indifference toward which smarticle is inactive, we plotted the cumulative distribution of total inactivity time in the form of a Lorenz curve (Fig. 7C). The curve presents the share of inactive time covered by the smarticle spending the least time being inactive (39). The shape of the Lorenz curve reflects the inequality in the distribution of the inactivity times of the smarticles: The more concave the curve, the more unequal the distribution. To characterize the Lorenz curve, we introduce the Gini coefficient, G, defined as the ratio of the area between the Lorenz curve and a line representing equality to the total area under the line of equality (39, 40). A value of 0 represents equality, and a value of 1 is perfect inequality. In a single endogenous experiment lasting 25 min where the light changed directions five times (see movie S6), we found that G = 0.21 (see bolded line in Fig. 7C). The Lorenz curve using data from all trials shows that 57% of the inactivity time was accounted for by 43% of the smarticles.

By considering each of the five excursions independently, the Gini coefficient and Lorenz curve can change markedly (see unbolded lines in Fig. 7C). For singular excursions, certain smarticles may remain in the inactive position for extended periods of time with a static light source, thus giving the Lorenz curves high values of inequality. This is a result of aforementioned correlations that can happen in smarticle collisions. On shorter time scales, the correlations may incorrectly lead one to believe a smarticle hierarchy exists; however, on long time scales, it becomes apparent that the smarticles are indeed commutable.

Thus, although crude, the endogenously drifting supersmarticle result demonstrates that the collective can perform a task/behavior (40) such that locomotor control of the system is decentralized and offloaded completely to mechanical interactions (41) in response to highly structured environmental signals (i.e., smarticle inactivity patterns).

Discovering emergent control authority

Most control synthesis techniques from the past six decades rely on a deterministic understanding of actuation and its effect on system states (42). But to create an organized system out of disorganized components, it is necessary to understand what the collective can accomplish as a function of uncertain subsystem interactions (43). To enable the discovery of control strategies for collective locomotors, control must be synthesized with respect to ensemble properties rather than individual states.

We expect an ensemble’s control authority to be an emergent property rather than intrinsic. To address this, we introduce the notion of a candidate control signal to hypothesize actuation mechanisms based on broken symmetries in the system (44). Using control signals, we can take a system with symmetry—and associated conserved quantities—and apply control to break the symmetry, thereby asserting authority over otherwise conserved quantities. When actuation mechanisms are unknown, symmetry breaking can be used as a way to hypothesize candidate control signals contributing to a system’s emergent control authority.

Given a candidate control signal, we apply a nonparametric, unsupervised learning algorithm, dynamical system segmentation (DSS), to discover ensemble-level behaviors in relation to the signal. DSS extracts distinct system dynamics from the interactions of internal states and, when present, the effect of candidate control signals on states (45). Initially, the algorithm constructs a set of system models over sequential windows in time—each locally capturing the net effect of interactions between internal states and candidate control signals on the ensemble dynamics.

In constructing the set of models, we represent elements of the set using the Koopman operator, K—an infinite-dimensional linear operator describing measure-preserving nonlinear dynamical systems through the evolution of observables (46). This choice of model is important because the Koopman operator does not explicitly require state information to describe the evolution of the system. Instead, the operator depends on observables, which may be any time-varying sensor measurement or property of the system such as mass, volume, or temperature. Formally, observables g are real-valued functions drawn from an infinite-dimensional Hilbert space, ℍ, that take measurements as their argument. The evolution of an observable through the infinite-dimensional Koopman operator isKg(xk)=g(F(xk))=g(xk+1)(3)where K:HH acts directly on the observables in the function space. We approximate K in finite dimensions with a data-driven operator K : ℋ→ ℋ by choosing a basis for some subspace ℋ ⊂ ℍ and applying least-squares optimization to compute K. The finite-dimensional operator, K, is then an N × N matrix for a given choice of N-dimensional basis.

The algorithm then condenses the set of Koopman operators into a set of non-redundant exemplars by applying nonparametric clustering (47) directly onto the set of models—where each element is itself a matrix. The resulting compressed set contains all unique dynamical system behaviors observed in the dataset of sequential measurements. DSS achieves this without assuming how many behaviors the system exhibits—an important property when the cardinality is generally unknown a priori.

The output of DSS is a set of distinct, yet related, ensemble behaviors represented by a probabilistic graphical model G=(K,E). The graph’s node set is specified by the compressed set of system behaviors, and its edge set is determined empirically by the transitions observed in the training dataset. In this model, the ensemble behaviors, each of which is a deterministic description of the ensemble dynamics at a given configuration, are random variables whose joint probabilities are in G. We refer to the information encoded by G as the system’s behavioral patterns.

Although the literature of learning control is evolving rapidly, existing methods are not immediately well suited for a problem as ill-posed as discovering emergent control authority. For one, reliably designing payoffs to reward emergence may not be possible, making it difficult to directly apply most reinforcement learning approaches (4850). Moreover, techniques in inverse reinforcement learning, such as learning from demonstration (51) and imitation learning (52), typically suffer from a lack of generalizability, limiting the use of learned behaviors. DSS avoids these pitfalls by directly analyzing distinct system behaviors and constructing a predictive model from these subsystem interactions, leading to a generalizable model. In addition, DSS is extremely data efficient, which is critical given that emergence is typically a rare phenomenon.

Decentralized control of supersmarticles

On the basis of observations made in previous sections, we know the switching sequence of inactive smarticles (Fig. 7B) is causally related to system behavior via the breaking of symmetries in the internal collision distribution of the supersmarticle. We used this sequence as a candidate control signal and modeled its effect on supersmarticle dynamics with DSS. By taking data from a single endogenous phototaxis demonstration (such as Fig. 7A), we instantiated two separate models of the supersmarticle dynamics with DSS—one with candidate control information and one without—and studied their respective behavioral patterns. The basis functions used in DSS were selected on the basis of their ability to represent information about the relative locations of inactive smarticles within the ring and their effect on the motion of the supersmarticle (see the Supplementary Materials for details).

The resulting graphical models are shown in Fig. 8 (A to C) along with a graph constructed from the observed smarticle switching sequence. We refer to the switching sequence of inactive smarticles observed in the given experimental trial as the experimental inactivity patterns, shown in Fig. 8A as a graph. Each node in the graph of Fig. 8A represents a unique combination of smarticles that were experimentally observed to be simultaneously inactive, and the number in each node refers to the unique label of the respective inactive smarticles. For example, in Fig. 8A, the green node labeled 4 indicates that, at some point in the experimental trial, smarticle number 4 was inactive, and subsequently smarticle 1 also became inactive, changing the supersmarticle to the dark blue node labeled 1 and 4. We used the complexity (i.e., graph complexity) of the experimental inactivity patterns graph to represent an estimated baseline complexity of supersmarticle behaviors. If the candidate control signal was causally related to the ensemble dynamics, the system should have responded to actuation leading to behaviors identifiable by DSS.

Fig. 8 Comparison between observed and numerically generated behavioral patterns.

(A) Model shown is a graphical representation of a single inactive smarticle switching sequence observed in a phototactic experiment (such as Fig. 7A). We extracted graphical models using the DSS algorithm from the same experiment with and without candidate control information, shown in (C) and (B), respectively. (D) Simulated supersmarticle trajectories predict that the ensemble is capable of movement anywhere in the plane while receiving exogenous feedback from an external controller. (E) Experimental trajectory of a photophilic supersmarticle in which the system was exogenously steered through a maze by an experimenter (movie S7), validating the simulation’s predictions. The trajectory evolves in time from blue to red, and the black circles represent the initial and final ring configurations.

Without candidate control information, DSS is unable to identify a set of behaviors explaining the observed drift, as seen in the nominal behavioral patterns in Fig. 8B. However, when candidate control information is incorporated, DSS extracts a set of emergent behavioral patterns of equal cardinality to the inactive smarticle switching sequence (Fig. 8C), where numbered nodes correspond to distinct behaviors identified by DSS. We note that there is no one-to-one correspondence between identified behaviors and elements of the candidate control signal. This is to be expected because the ensemble dynamics are not exclusively determined by the candidate control signal and are also driven by uncertainty, which is captured by the probabilistic transitions in the generated graphs.

Through endogenous steering, we showed that the supersmarticle is capable of directed transport toward a fixed objective. To improve on the controllability of the ensemble shown in previous sections, we looked for alternative locomotion strategies in simulation by synthesizing control to directly manipulate the algorithmically extracted ensemble behaviors. Previous work in control of interconnected stochastic systems has shown that integral control can often be a simple and robust strategy (53). However, due to the discontinuous nature of supersmarticle inactivations, model-based control is necessary to directly optimize system actions.

We designed a simple decentralized model-predictive controller that greedily searches for inactive smarticle switching sequences to alter the ensemble’s behavior—as determined by its DSS model—to achieve a collective objective. The collective objective was expressed as a quadratic cost on the position of the supersmarticle centroid with respect to a desired goal location in the world frame. The supersmarticle centroid was calculated via distributed consensus in a fully connected topology (see the Supplementary Materials for additional details) (54). We conducted four sets of Monte Carlo simulations (40 trials each), over distinct goal locations—left, right, up, and down—with randomized initial conditions for a duration of T = 75τ per trial (where τ = 1.6s). The objectives were always equidistant and located directly vertically or horizontally from the initial conditions.

The resulting trajectories shown in Fig. 8D were exogenously steered in the world frame by the independent decision-making of individual smarticles. The simulation results confirm the symmetry-based theoretical predictions: The ensemble should be capable of locomotion anywhere in the plane via exogenously selected smarticle inactivations—even when we train the model using only a single trajectory moving in a single direction. The supersmarticle provides a test case for whether DSS can detect emergent behavior and whether DSS (or related algorithms) should be used in more general settings where symmetry-based inference about control authority is not possible. By allowing an external source of feedback to inactivate smarticles, the decentralized controller manipulated ensemble behaviors to achieve more complicated goals than the model trained on, thereby predicting entirely emergent behavior. As a result of the generalizability of our machine learning model, we are able to make predictions and control the supersmarticle in entirely new settings, thereby harnessing the system’s emergent control authority to accomplish brand new tasks. We note that extending the smarticle hardware to accommodate for the proposed control algorithm can be done in a practical and computationally efficient way and will be explored in future work.

On the basis of the simulation results, we experimentally validate the exogenous controllability predictions by guiding the supersmarticle through a simple maze using external feedback from an experimenter with a light source (Fig. 8E). Here, the experimenter is capable of directing the supersmarticle by freely shining a light source onto the ensemble, thereby using more complex inactivity sequences to achieve locomotion anywhere in the plane, just as the proposed decentralized control scheme did. Although the supersmarticle was provided with external guidance, it was able to achieve directed transport without state information or specifying individual objectives for its constituents. All movement was directly emergent from morphological computations in response to environmental signals (55). Hence, by framing the discovery of emergent control authority as a learning problem, we were able to hypothesize and model unconventional actuation leading to expressive controllable motion.


Inspired by a future in which a class of task-capable robots could be formed from myriad redundant and task-incapable components, we have created a primitive “robot made of robots” that can perform rudimentary phototaxis, despite none of its components—smarticles—having locomotor capabilities. A generic statistical model accurately captures the fundamental drift dynamics, rationalizing how the supersmarticle can sense an aspect of its environment—light—and use this to endogenously steer itself via asymmetric inactivation of individuals. Further, through the introduction of novel machine learning techniques, we constructed a data-driven model of the ensemble, which enabled discovery and proof-of-control alternatives for generating exogenous steering when agents are capable of computation.

We emphasize that, unlike other mobile robots, the supersmarticle displays phototaxis without a central processor or dedicated motor components. The key ingredient, and what differentiates our collective from other robot swarms and locomoting collectives, is that our system is made of components that have very low control authority—they cannot locomote individually—and are highly unpredictable—they create emergent behavior from the highly complex interactions of their internal degrees of freedom. Although such a system might seem idiosyncratic, we note that it bears similarities to cascades of conformational changes in the nanomachines that regulate many cellular processes: proteins (56). Given the ubiquity of such processes in these tiny machines, we posit that our model smarticle system could provide inspiration for the generation of substantially more complex task-capable ensembles like those pictured in Fig. 1, including perhaps three-dimensional (3D) collective locomotors and manipulators. Enabling robots to flexibly reconfigure to collectively perform tasks in the presence of environmental noise and individual component malfunction or degradation (23) could enhance robustness in robot swarms across scales, from intravenous delivery (5759) to search and rescue (60). Further, insights from collective robophysical systems (15, 61, 62) could elucidate principles by which biological collectives [like slime molds (63)] perform tasks in complex natural environments.


Smarticle robots

Each smarticle’s outer shell and arms, or outer links, are 3D printed. The arms are controlled by HD-1440A servomotors to a precision of <1 and with an accuracy of ±6. All processing and servomotor control is handled by an Arduino Pro Mini 328 (3.3 V/8 MHz model), which allows smarticles to be programmed to deform to specific configurations and gaits, where we define gaits as periodic trajectories in the configuration space (see Fig. 2B). When assembled, each single smarticle has a mass m = 34.8 ± 0.5 g. The system is powered by a 3.7-V, 150-mAh, 30-C LiPo battery (Venom; Rathdrum, ID) enabling hours of testing. Smarticle positions and orientation were tracked using an infrared video recording hardware/software suite (OptiTrack; Corvallis, OR). All experiments were conducted on a 60 cm–by–60 cm aluminum plate leveled flat to <0.1.

Smarticle experiments

The Gini coefficient (G) is a statistical measure derived from the shape of the Lorenz curve. A value of G = 0 represents a situation of perfect equality, or in the case of the supersmarticle, all smarticles spent an equal amount of time being inactive. Conversely, a value of G = 1 is a maximally unequal trial, or one where only a single smarticle was inactive over the course of the experiment.

Statistical model

Below, we detail the full form of the variance of vVar[v]=(f/4π3T)[π2sin(2Ψ)(RA2RB2λ)+RA2((4λ3+2λ2+2)Ψ+π3((4λ2+2)/(π2)+λ2+(2Ψ)/π))4RARBλ(2λ+1)(λΨ+Ψ+π)RB2λ(6(λ1)λΨ6πλ+2π2Ψ+π3)2(λΨ+Ψπ)(RARBλ)(4λ(RARB)sin(Ψ)+cos(2Ψ)(RARBλ))]

Furthermore, the full form of variance of v is shown belowVar[v]=(f/4πT)[RB2λ(π+2Ψ)RA2(π(λ2)+2Ψ)+(RA2RB2λ)sin(2Ψ)]

To calculate the values of RA and RB, we must start with masses m1 and m2 such that the relative distance between them, x1x2, is specified by the actuation of the smarticles. The first mass, m1, represents the arm of a smarticle, and m2 represents the body. The mass of the boundary they push on is mb. Both m2 and mb have friction between them and the surface they are sitting on. This is shown in Fig. 6A. On the basis of this model, we arrive at the following equations of motionF2fs=m2x¨2fsFb=(mb+m1)x¨1(4)where F2 and Fb are the friction force on m2 and mb, respectively, and fs is the force between m1 and m2. By specifying x1x2 = A0 sin (ωt + γ), F2 = (m2 + 2m1)gμ, and Fb = mbgμ, these equations can be integrated to find how far mb moves. Then, by plugging in for mb—the mass of just the ring—and the mass of the ring and the inactive smarticle, we can find RA and RB, respectively, as well as 〈v∥, ⊥〉 and Var[v∥, ⊥].

Dynamical system segmentation

The DSS algorithm is composed of three primary subroutines: (i) the calculation of Koopman operators over sequential windows of time via least-squares optimization, (ii) nonparametric clustering over the space of Koopman operators to determine unique system behaviors, and (iii) training a supervised learning model [e.g., support vector machine (SVM)] to learn relationships between system behaviors and construct the complete probabilistic graphical model. In the following sections, we expand on these subroutines, and a full outline of the algorithm can be found in the Supplementary Materials.

Koopman operators

The DSS algorithm first requires calculating finite-dimensional Koopman operators over sequential windows of the dataset. Although there are many ways to frame Koopman operator synthesis, we implement it as a least-squares optimization (64). Given a choice of nonlinear basis function ψ(x) and a data sample X = {x1, …, xM}, we can formulate the Koopman operator synthesis problem as solvingminK12k=1M1ψ(xk+1)Kψ(xk)2

This optimization has a closed form solution of the following formK= AGwhere † denotes the Moore-Penrose pseudoinverse, and the individual matrix components areG=1Mk=1M1ψ(xk)ψ(xk)TA=1Mk=1M1ψ(xk+1)ψ(xk)T

Nonparametric clustering

Given a set of Koopman operators, DSS looks for distinct dynamical behaviors by applying nonparametric clustering directly onto the set of operators. In particular, we apply hierarchical density-based spatial clustering of applications with noise (HDBSCAN) (47), which is a nonparametric clustering algorithm that specializes in problems subject to noisy and sparse measurements. By using HDBSCAN, we were able to discern distinct behaviors from the set of Koopman operators. From these clustered classes, we constructed class exemplars as a means of creating a set of distinct Koopman operators corresponding to observed system behaviors.

Supervised model

Once DSS has compiled a condensed set of exemplar behaviors, the algorithm must then determine the dependencies between each behavior and the states of the system. To this end, we trained an SVM. We did this by using the clustered class labels from HDBSCAN to label the state-space data. Then, using this newly labeled dataset, we trained a soft-margin SVM that assigned discerned behaviors to state observations. The SVM, in conjunction with the condensed set of exemplar behaviors, gave rise to the probabilistic graphical model, where the dynamics of the system are described by stochastically shifting Koopman operators.



Fig. S1. Unrotated center of mass trajectory of the smarticle cloud.

Fig. S2. Unrotated trajectories of the supersmarticle.

Fig. S3. Theoretical and experimental data for the perpendicular component of the supersmarticle drift speed.

Table S1. List of all six different types of collisions in the theoretical model.

Table S2. List of all parameters used in the theoretical model.

Algorithm S1. Dynamical system segmentation.

Movie S1. Individual smarticle performing square gait.

Movie S2. Smarticle cloud: Seven active smarticles.

Movie S3. Supersmarticle: M = 0.51, five active smarticles.

Movie S4. Supersmarticle: M = 0.51, one inactive, four active smarticles.

Movie S5. Supersmarticle: M = 3.6, one inactive, four active smarticles.

Movie S6. Supersmarticle: M = 3.6, endogenous phototaxing.

Movie S7. Supersmarticle: M = 3.6, exogenous phototaxing.


Acknowledgments: We thank M. Caveney for experimental assistance and D. Randall, A. W. Richa, A. Zangwill, J. M. Rieser, and P. Umbanhowar for helpful discussions. Funding: Support for W.S., S.L., R.W., and D.I.G. was provided by NSF PoLS-0957659, PHY-1205878, DMR-1551095, and ARO W911NF-13-1-0347. Funding support for Z.J. and K.W. was provided by PHY-1205878. Funding for T.A.B. and T.D.M. was provided by NSF CBET-1637764, and funding for A.P. was provided by an NSF GSRF. Author contributions: W.S. designed and performed all experiments with the help of R.W. and S.L.; Z.J. designed the theoretical model; and T.A.B. and A.P. designed the modeling algorithm and controller, and performed numerical simulations. W.S., Z.J., T.A.B., K.W.,T.D.M., and D.I.G. designed the research and wrote the paper; D.I.G. guided overall research program. Data and materials availability: All files needed for constructing smarticles can be found in, and all files needed for applying the DSS algorithm to the supersmarticle, along with representative data samples, can be found in

Stay Connected to Science Robotics

Navigate This Article