Distributed prediction of unsafe reconfiguration scenarios of modular robotic Programmable Matter

We present a distributed framework for predicting whether a planned reconfiguration step of a modular robot will mechanically overload the structure, causing it to break or lose stability under its own weight. The algorithm is executed by the modular robot itself and based on a distributed iterative solution of mechanical equilibrium equations derived from a simplified model of the robot. The model treats inter-modular connections as beams and assumes no-sliding contact between the modules and the ground. We also provide a procedure for simplified instability detection. The algorithm is verified in the Programmable Matter simulator VisibleSim, and in real-life experiments on the modular robotic system Blinky Blocks.


I. INTRODUCTION
Materials able to autonomously assume any shape are the dream of engineers. Currently, the most advanced artificial systems possessing elements of this functionality are modular self-reconfigurable robots-machines composed of robotic units (modules) which can bond together, move over one another and communicate, as well as store and process information [1]. They may be compared to swarms of fire ants forming bio-mechanical structures from their own bodies [2]. Cooperation of millions of tiny, densely-packed modules is expected to produce the desired emergent shape change of the entire ensemble. A system of this kind would be a realization of the futuristic concept of Programmable Matter [3].
The operation of densely-packed self-reconfigurable robots is based on the movement of modules from one part of the robot to another. This poses not only the challenging hardware problem of designing miniaturized modules able to operate in large 3D ensembles, but also the software problem of controlling this motion. If the robot is autonomous, its operation must be collectively planned and controlled by the modules, taking into account geometric and mechanical constraints on reconfiguration at each stage of motion. Although offline approaches could be used to check these constraints efficiently for a number of selected configurations, the less efficient online approaches must be used in general to handle any shape or any possible interaction with obstacles.
Geometric constraints result from the fact that modules need sufficient space to make a planned move but also must constantly remain in physical contact with other modules [4]. Mechanical constraints, in turn, result from the requirements of integrity and stability of the entire robot-the structure cannot break at inter-modular junctions or lose balance during reconfiguration. The mechanical constraints can be neglected in some special cases, like reconfiguration under no external loading (weightlessness) or 2D reconfiguration on flat ground with perpendicular gravity as the only loading. Otherwise, they usually need to be considered.
Currently, almost all algorithms for planning and controlling self-reconfiguration of densely-packed systems take only geometric constraints into account, like in [5], [6], [7], [8], [9]. An up-to-date survey of self-reconfiguration algorithms for modular robots can be found in [10]. By contrast, works on the mechanical behavior of densely-packed modular structures are few, present mostly centralized procedures, and rarely discuss reconfiguration. Examples can be found in [11], focused on optimizing the compliance of modular tools, and in [12], aimed at predicting the mechanical behavior of structures produced by additive manufacturing. As a separate research direction, special types of modular structures were investigated in [13], [14] and [15] to check the possibility of using modular robots as collective actuators. Autonomous reconfiguration planning that takes into account both geometric and mechanical constraints remains a challenging open problem, some aspects of which are investigated in this work.
In the present paper, we develop the approach introduced in [16] much further into a more realistic framework. Arbitrary 3D structures are investigated, for which a linear-elastic FE model is again adopted, with the addition of unilateral contact conditions which represent interaction with the surroundings. Two failure modes are considered: overloading of intermodular connections and loss of balance, both checked in a distributed manner. Two methods of checking the loss of balance are proposed: (1) a simplified one, valid for structures standing on a flat surface, and (2) a model-based one, which is more general but requires solution of the mechanical balance equations with contact conditions. Verification is performed in a dedicated simulator VisibleSim [17], as well as experimentally on the real robotic modules Blinky Blocks [18]. The computational cost and several possible extensions of the applied weighted Jacobi iterative solver are also discussed.

A. Basic characteristics of a modular robot
For the purposes of algorithmic mechanical analysis, a modular robot will be represented only by its connection topology and inter-modular connection strengths. We will focus on structures built of cubic Blinky Blocks [18], which are arranged on the Cartesian grid. However, in principle, the proposed approach can handle other connection topologies.
From the information-theoretic point of view, a modular robot has a distributed and asynchronous computing architecture, with mesh connection topology corresponding to the communication network of modules [19]. In Blinky Blocks, information can only be directly transferred between adjacent modules, so the communication network has the same connection topology as the modular robot itself. Such neighborto-neighbor-only communication imposes strict limits on the information-passing abilities of the system and requires special algorithms for effective prediction of mechanical failures.

B. Problem to be solved
Reconfiguration of a modular robot can be broken down into simple steps. A single step consists of (a) releasing some inter-modular connections, (b) shifting modules which have been freed into nearby positions and, finally, (c) creating new connections at the new locations. As an alternative way of restructuring one may consider attaching new modules to a structure and discarding some old ones. Each reconfiguration step can potentially cause failure of the structure, which is in general irreversible. To avoid it, a prior mechanical analysis can be performed to ascertain that the planned new configuration is mechanically safe.
We restrict further presentation only to the case of attaching additional modules to an existing structure. Nevertheless, the whole idea can also be applied to a complete reconfiguration step. In such a case, the mechanical analysis would need to be performed for each of the above stages (a), (b) and (c) of the planned reconfiguration step. The step would be considered permissible if the structure was predicted to be safe after each of the three stages. Since this kind of reconfiguration cannot be easily performed on Blinky Blocks, we do not further develop this topic in the present paper.
Algorithmic prediction of two failure modes, shown in Fig. 1, will be analyzed, expanding on the ideas proposed in [16]. The first one is when addition of new blocks increases stresses at some inter-modular junctions beyond the holding capacity of connectors, as a result of which the structure breaks, Fig. 1a. The second one is when the structure loses balance when the new modules are added, Fig. 1b.

C. Overview of model-based failure prediction
We propose a distributed procedure which can predict both types of failure simultaneously. It approximately solves a special mechanical problem for the robot with new modules attached and can be customized to handle different module designs. The procedure has several distinctive features, described in detail in the succeeding sections:

D. Standard 3D frame model
In the adopted FE model, each module p is represented by a node with 6 degrees of freedom u p (3 displacements, u x , u y , u z , and 3 rotations, τ x , τ y , τ z ) and each pair of connected modules is represented by a beam joining their nodes. A module in contact with the ground is represented by a special beam between the module's node and a "ground" node g, with u g = 0, as it is explained in Sec. II-F. In Fig. 1, the beams are presented as lines joining the centers of adjacent modules. In a coordinate system CS whose z axis points upwards, the stiffness matrices for a beam joining module p and module q lying below it read K 11 pq = K 11 and K 12 pq = K 12 , where: while E, L, A, I x , I y and J ν are the elastic modulus, length, cross-sectional area, area moments of inertia in the x and y directions, and scaled torsion constant in the z direction of the beam, respectively (see Tab. I). If the two neighbors p and q are aligned with the z axis of another coordinate system CS', then K 11 pq and K 12 pq take the form: R pq is the 3 × 3 rotation matrix from CS' to CS, T denotes a transpose, and block matrix notation is used. The ensemble is in static equilibrium if, for every module p, the sums of reaction forces and torques between p and all its neighbors q are equal to the gravitational force and null external torque acting on p, respectively, given by the vector F ext p = [0, 0, 9.81 · M, 0, 0, 0] T . The equilibrium equations are: with unknown vectors u p and u q . After assembling the global vectors u, F ext and stiffness matrix K, incorporating all degrees of freedom of the structure, Eq. (5) takes the form Ku = F ext . It still needs to be extended to account for the modules planned to be added (Sec. II-E) and contact conditions (Sec. II-F). The respective quantities for the extended system, called the perturbed state, will be denoted by symbols with bars.

E. Adding virtual modules
To predict the state of the system one reconfiguration step ahead, the algorithm must take into account virtual modulesthe modules which are planned to be attached. This is done in a simple way: a system analogous to Eq. (5) is built assuming that virtual modules are present. During computation, virtual modules are emulated by their existing neighbors, which store and process variables and messages related to virtual modules.

F. Contact with external supports
We use a simplified contact model of a cubic module with an external support, in which the support must be flat and coplanar with one of the module's facets. Also, we only analyze initially existing contact interfaces and assume that no new ones appear under load. We distinguish two conditions: (i) a unilateral contact-separation condition (coaxial mode) combined with no-sliding/no-twisting (shear and torsional modes), (ii) a tilting condition (bending mode). We say that a module is in contact if the axial force in the beam representing the contact is compressive. Otherwise, the module is in separation. When in contact, the module can only tilt if the bending torque in the beam exceeds a limit torque which is proportional to the compressive force (see the explanation below).
Without loss of generality, let us consider a module supported from below. The module is loaded with the forces and Fig. 2b) takes the form of the Signorini problem: Additionally, we require that when the contact is active ( f z < 0), there is no tangential slip (u x = 0 & u y = 0) or twisting (τ z = 0), and the module can only tilt if at least one of the bending torques exceeds a limit torque. The tilting conditions (see Fig. 2c) can be conveniently written in the following form: where L is the module's size. The tilting (bending) condition expresses the fact that the torques |m x | or |m y | cannot exceed the torque f z · L/2 produced by the compressive force f z about any of the facet's edges; see also Fig. 2a.
A supported module is modeled as a beam with the same elastic constants as two connected modules. This gives the force-displacement relationship for the case of stable contact (when the module adheres to the support with an entire facet). This relationship is then used by a special predictor-corrector scheme, described below, which provides regularization to the contact conditions given by Eqs. (6)(7)(8).
In the regularized contact scheme for module p, a trial state is computed first, assuming a linear-beam-type connection with each support q (here, the contact direction is z): where the term K 12 pqūq is absent because the support is assumed to be immobile,ū q = 0. Then, a corrected vector F pq is determined, taking into account two conditions: (i) the unilateral (normal) contact-separation condition combined with no-sliding and no-twisting conditions: ) (ii) the tilting (bending) condition, used only when f tr z < 0, i.e. when the module is in contact: x ,m tr y , Φ tr x and Φ tr y are computed as in Eqs. (7) and (8), but using the components of F tr pq ; γ = 10 −4 andγ = 1 − γ. The corrected contact tangent matrixK 11 pq is obtained as a derivative of F pq with respect toū p . Again, the matrixK 12 pq is disregarded because supports are assumed to be immobile.
Remark: Ill-posedness of the problem is avoided by introducing a very weak spring, characterized by γ, that prevents free rigid-body motion of the structure. The drawback of this approach is poor conditioning of the resulting system of equations when the robot is unstable, which can deteriorate the convergence rate of the iterative solver; see Sec. II-H. However, as it is shown in Sec. III-C, the knowledge of which supports are active suffices for assessing stability and those are usually identified long before convergence is achieved.

G. Weighted Jacobi solution scheme
The global system of equations for the perturbed system, in which virtual modules are included in the model (Sec. II-E) and contact conditions are accounted for by the predictorcorrector scheme (Sec. II-F), reads: Eq. (13) is solved iteratively using the weighted Jacobi scheme [16]. A single iteration i → i + 1 for module p reads: whereD p = diag ∑ qK 11 pq andR p = ∑ qK 11 pq −D p are the diagonal and the remainder parts of the respective stiffness sub-matrices, while β = 2/3. Initially, we setū 0 = 0 (alternatively, the solution for the non-perturbed state, if available, could be used instead of 0, which would reduce the necessary number of iterations). Note that in the iteration i + 1 only the values ofū from the iteration i are used, and only those from p and its direct neighbors q. Thus only local communication is involved and the memory complexity is constant.
In the present implementation, the weighted Jacobi procedure is initiated by the centroid module which sends an Init message down the spanning tree, broadcasting the number of iterations to be done (see Sec. III-B for the centroid module and the spanning tree). Having received the Init message, each Blinky Block B p sends its initial vectorū 0 p = 0 to all   its neighbors and initializes its iteration counter, iter p = 0. In a given iteration iter p = i, B p can receive from any of its neighbors B q a message containingū i q (displacements of B q calculated in iteration i), which is then stored in B p 's buffer. When B p has receivedū i q from all its neighbors, it computes u i+1 p (see Eq. 14), increments its counter, iter p ← iter p + 1, and sendsū i+1 p to all its neighbors. The process continues until the prescribed number of iterations is reached.
Remark. The weighted Jacobi procedure behaves like the Alpha local synchronizer [20].

H. Convergence properties and possible improvements
The Weighted Jacobi scheme converges if the spectral radius of the iteration matrix C β = I − βD −1K is less than 1: where λ i are the eigenvalues of C β . Although in the cases analyzed here convergence is achieved, the number of iterations is very high, which is a well-known drawback of the method (the convergence rate tends to 1 when the system grows [21]). For the one-dimensional spring-in-series system of size n, we have analytically assessed the number of iterations necessary to attain an arbitrary relative error to be O(n 2 ). This is also confirmed numerically in Sec. IV-B for more complex structures. The assessment shows the low efficiency of the scheme and underlies the complexities provided in Table II. The framework presented in this work is not restricted to the weighted Jacobi scheme though. Its efficiency can potentially be significantly improved by adapting another method to solving the considered contact problem in a distributed way. We will briefly outline the three most promising directions.
Direction 1: The Krylov subspace methods [22] guarantee that the maximal number of iterations is at most equal to the number of degrees of freedom of the system, if the problem to be solved is linear. This can be further improved by appropriate preconditioning (see our preliminary study [23]). However, the need for global data aggregation and the non-linearities introduced by the contact problem require special treatment, which can deteriorate the time and memory efficiency.
Direction 2: Multigrid techniques [24] can more easily capture long-wave modes of the solution, which should improve the convergence rate. A special version must be devised, however, taking into account the contact conditions (like the one recently proposed [25]) and the specific computing architecture of the modular robot.

Direction 3:
The number of degrees of freedom of the system can be reduced by applying multi-scale methods or model order reduction techniques [26], [27]. However, it may be hard to find a suitable reduced space online efficiently.

III. MECHANICAL STABILITY AND OVERLOAD CHECK
Below we describe computational methods using which a robot can autonomously predict the two types of failure shown in Fig. 1, one reconfiguration step ahead. The methods utilize a spanning tree, which is discussed first in Sec. III-A. Sec. III-B describes a simple method of checking stabilitywithout iterations, but restricted to robots standing on flat ground. Sec. III-C discusses stability verification in the general case, using the iterative scheme of Sec. II. Finally, in Sec. III-D, conditions for inter-modular connection breakage are presented, utilizing the results of the iterative scheme. The flowchart of the procedure is shown in Fig. 3.
Assessments of complexities of the subroutines of the algorithm are presented in Table II, where n is the number of modules, d is the radius of the connection graph of the robot, andd is the depth of the spanning tree (usually,d and d are of the same order; see Sec. III-A for more details).

A. Spanning tree
A spanning tree allows efficient communication inside the robot. Its construction begins with a choice of the centroid module, serving as the root, which is selected near the robot's topological center; see e.g. Fig. 4a. This can be done automatically [28], but we chose the centroid manually in all examples.
The tree is extended to all modules, starting at the centroid which sends a Tree message to its neighbours. When a module receives the Tree message for the first time, it becomes a next-level node and sends the Tree message further. This usually leads to the construction of BFS-like trees of quasioptimal depth, without the need for synchronization. However, in rare cases the resulting tree may be far from optimal, with a Hamiltonian path over the structure as the worst-case scenario. The algorithms relying on the tree are then adversely affected.

B. Loss of balance in a simplified case
It is assumed that the robot is rigid and stands on flat ground under vertical gravity (Fig. 4a). The modules know their own masses and positions in a common Cartesian coordinate system, with z = 0 being the ground level. They also simulate the behavior of their virtual neighbors (Sec. II-E).
The stability check reduces to verifying that the center of mass of the robot lies over the convex hull of the points of support. If it does, the robot is stable, otherwise, it is not. The algorithm proceeds as follows.

C. Loss of balance in the model-based approach
In the general case with arbitrarily placed supports, there is no simple method to predict stability. Sometimes, if a model of the robot is simple, like under the rigid-body or elasto-static assumptions used here, there may even be no unique answer. Since more accurate modeling goes beyond the scope of the present paper, we will show how to utilize the proposed elastostatic model with contact and the iterative solution scheme to check stability in more general cases. The method, however, does not assure finding the correct solution in difficult cases (e.g., when the solution is non-unique by definition).
The method is based on the observation that, when the iterative solution scheme has converged, the local state of the contact conditions is fully determined. It is only necessary to check whether active supports prevent rigid-body rotation of the structure (at least three noncollinear support points must be active), which is achieved by the following procedure: (a) Solve the mechanical problem (see Sec. II) and initiate the stability check (the spanning tree is used again). (b) For every module, determine the set of its active corner points by checking the contact conditions (Sec. II-F): • No contact → return the empty set.
• Tilting in two directions → return a single corner-the common point of the two edges of rotation. • Tilting in one direction → return two corners-the end points of the edge of rotation. • Contact & no tilting → return a special 'stable' state. (c) Aggregate information starting at the leaves and moving up the spanning tree towards the centroid: • At each module, sum the sets of active points of its subbranches and its own, obtaining set S. By convention, adding any set to 'stable' gives 'stable'. Excluding the expensive phase of determining active supports (weighted Jacobi iterations), the complexity of the approach is the same as that of the simplified case; see Table II. The memory complexity per module is constant because each returned set of points has at most two elements.

D. Overloading of inter-modular connections
Connection overloading is checked when the iterative scheme has sufficiently converged and after checking that the structure is stable. The forces and torques which act between module p and its neighbor q are predicted as follows: whereR T pq rotates the resulting vector into a coordinate system in which axial forces are aligned with the z axis.
To avoid connection breakage, the tensile force f z and torques m x and m y , computed in Eq. (16) in the middle of the connection, must not overpower the magnetic forces F max binding the modules. (Shear and torsion are omitted, because Blinky Blocks' connectors are assumed to be resistant to those modes of breakage.) The vertical and lateral connections of Blinky Blocks differ, so that F max can take two values; see Table I. The safety condition for both tension and bending is: The check is performed for all connections and the results are aggregated by the centroid over the spanning tree.

A. Implementation details
The procedures have been implemented and tested in the integrated environment developed at FEMTO-ST 1 . It consists of the virtual test bed VisibleSim [17] and the reconfigurable modular robot Blinky Blocks [18], so the same implementation could be executed on both platforms. The software was adjusted to be compatible with the real Blinky Blocks Version 1 hardware: reduced floating-point precision of 4 Bytes was used and the maximum message size was set to 17 Bytes (messages containing 6 × 4 Byte-long vectors were split in half).
The program's flowchart is shown in Fig. 3, and the consecutive steps of the algorithm are discussed in the previous sections. The choice between the simplified (Sec. III-B) and the full (Sec. III-C) stability check to be performed is preset by the user. In the case of the Blinky Blocks hardware, a preliminary step is additionally performed, in which the same main program is loaded into each Blinky Block, and a common coordinate system for a given configuration of Blinky Blocks is propagated, starting from a special block with preset coordinates.
The material parameters of the Blinky Block model described in Sec. II-D are provided in Table I. All have been assessed experimentally, except Young's modulus which was chosen arbitrarily (its exact value is not essential for assessing overload and stability). The dimensions and mass have been measured directly. Connection strengths have been obtained in a simplified manner, by finding the maximum number of Blinky Blocks that their magnets could hold hanging in a vertical alignment. The top/bottom and lateral connection strengths differ because the former is produced by a Legolike system reinforced with a central magnet, and the latter by 4 magnets placed in the corners of each face; see also Fig. 5.

B. Simulations and experiments
Six different modular configurations and failure scenarios were investigated (see Fig. 6), both in VisibleSim and on Blinky Blocks 2 . The sizes of structures ranged from 8 modules in test #1(a) to 29 modules in test #6(c). Experiments on substantially larger structures would be difficult due to the high computational complexity of the algorithm combined with the limited processing and communication speed of the Version 1 of Blinky Blocks. In the physical experiments, reconfiguration of Blinky Blocks was done by manually attaching new modules to an existing structure. In all the presented cases, the modelbased analysis was involved (addressing both overloading and loss of stability), which required execution of the weighted Jacobi iterative scheme. Additionally, in the loss-of-stability scenarios (Fig. 6 #1-#3), the results obtained with the simplified and the model-based stability checks were successfully cross-validated.  Because it is in general difficult to automatically assess the necessary number of weighted Jacobi iterations, this number was adjusted manually case by case. The criterion was to make the number of iterations possibly low while obtaining correct predictions at the same time. The problem of how the necessary number of iterations scales with the system size in a stable case is briefly discussed later. In an unstable case, this number is expected to be much higher. Following the conclusion in the final Remark of Sec. II-F, in unstable cases we stopped computations just after all contact states stabilized, but before numerical convergence was achieved.
Each of the six tests in Fig. 6 shows the results of execution of the program for two consecutive construction steps of a particular Blinky Blocks structure. In every figure, the first construction step (a→b) is designed to be mechanically safe, while the second one (b→c) to result in failure, which is then demonstrated in the third part of the figure (c). Additionally, in the top row, VisibleSim results are shown for the tests #1 and #4. The tests are split into loss-of-stability (left column) and overloading (right column) scenarios. From top to bottom, the scenarios are ordered by complexity, i.e., 2D cases in tests #1 & #4, 3D cases in tests #2 & #5, and 3D cases with more complex connection topologies in tests #3 & #6.
The results of calculations are displayed using colors: the color of a block corresponds to the highest tensile/bending stress level in any of its connections, as given by the right hand side of Eq. (17). Green to orange colors represent the safe stress range, while red indicates potentially overstressed connections. Blinky Blocks were programmed to blink in red when tensile stresses in some of their connections exceeded the critical level, while global imbalance of a structure was signaled by the centroid module blinking in purple. Blue Blinky Blocks are fixed-they are attached to the floor.
In all tests except #6, the predictions are confirmed by physical experiments. In test #6, breakage is correctly predicted but ill-localized. This can be observed in Fig. 6-#6(b) which indicates breakage of the pillars, while the actual breakage occurs as it is shown in Fig. 6-#6(c). There are two possible reasons for the observed discrepancy. The first one is that the assumed mechanical model of the modular robot is too simple. The second one is the omission of twisting torques from the adopted criterion of breakage. It was also very difficult to keep the structures #6(a) and #6(b) operational-an effect which was not expected. In both cases, weight-induced deformations caused separation of electrical connectors, despite the structures did not break. This necessitated using additional supports just to perform computations. We view test #6 as one of benchmark cases for future research on more accurate models and failure criteria.
CPU time and convergence properties. Computingū i p and exchanging messages with neighbours takes a Blinky Block a nearly constant time T ≈110.5 ms (9.05 runs per second). Because communication is local, T is also the global time of a single weighted Jacobi iteration, independent of the configuration. Since the time cost of the other steps of the algorithm is negligible (see Tab. II), the overall execution time can be assessed by multiplying the number of iterations by T .
The number of iterations needed to attain a given accuracy greatly depends on the system's configuration and generally grows with the number of modules n. Assessment of the number of iterations is generally not straightforward, even without considering unilateral contact conditions; see also Sec. II-H. In Fig. 7 we demonstrate the expected trends for a given family of configurations. Fig. 7a shows linear convergence of the relative error ū i −ū * / ū * as the number of iterations i grows, whereū * is the numerically exact solution. Fig. 7b presents the necessary numbers of iterations from Fig. 7a for two example relative errors, displaying quadratic growth with n and confirming the assessments in Tab. II.

V. CONCLUSIONS AND FUTURE RESEARCH
We presented a distributed algorithm for checking if a modular robot will retain its mechanical integrity and stability after new modules are attached to it at prescribed positions. The algorithm can be used to assess the mechanical safety of a reconfiguration step planned by a self-reconfigurable robot. The procedure is designed to run on the modular robot itself, and we have verified its predictions through tests in the dedicated simulator VisibleSim and on the real modular system Blinky Blocks. To our knowledge, this is the first time three-dimensional modular-robotic structures compute their mechanical state in a fully distributed manner.
The algorithm can be improved towards: adopting faster iterative schemes, as discussed in Sec. II-H and tried in [23]; extending the application range to soft modular robots; checking the construction/reconfiguration several steps ahead; as well as addressing other module geometries and broader module-support contact scenarios. Future experimental validation will use the currently produced new version of Blinky Blocks with faster CPUs and communication, and possibly quasi-spherical catoms [29] having up to 12 neighbors per module and electrostatic connectors.