# Distributed prediction of unsafe reconfiguration scenarios of modular
robotic Programmable Matter

Benoît Piranda  University of Bourgogne Franche-Comté, FEMTO-ST Institute,
CNRS, France  Paweł Chodkiewicz  Faculty of Automotive and Construction
Machinery Engineering
Warsaw University of Technology, Poland  Paweł Hołobut  Institute of
Fundamental Technological Research
Polish Academy of Sciences, Poland  Stéphane P.A. Bordas  Department of
Engineering
Faculty of Science, Technology and Medicine
University of Luxembourg  Julien Bourgeois  University of Bourgogne Franche-
Comté, FEMTO-ST Institute, CNRS, France  Jakub Lengiewicz  Institute of
Fundamental Technological Research
Polish Academy of Sciences, Poland  Department of Engineering
Faculty of Science, Technology and Medicine
University of Luxembourg

###### Abstract

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.

_Keywords:_ Self-reconfiguration, Modular robots, Distributed algorithms,
Mechanical constraints, Programmable Matter.

##  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 inter-modular
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.

##  II Distributed prediction of the mechanical state of a robot

###  II-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 neighbor-to-neighbor-only communication imposes strict
limits on the information-passing abilities of the system and requires special
algorithms for effective prediction of mechanical failures.

###  II-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.

(a)  (b)

Figure 1: Two types of failure after an additional module (gray dashed line)
is attached: (a) breakage of a connection; (b) loss of overall stability. The
red lines designate failing connections.

###  II-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:

  * •

The modular robot is represented by a Finite Element (FE) model; see Sec.
II-D.

  * •

Two types of connections are assumed: the inter-modular connection, modeled as
a linear-elastic beam, and the connection between a module and an external
support (e.g., the ground), modeled as a linear-elastic beam with unilateral
no-sliding contact conditions at the support end.

  * •

The mechanical state of a planned configuration with new modules attached (see
Sec. II-E) is obtained by solving a non-linear problem (non-linearities result
from the contact conditions; see Sec. II-F).

  * •

The problem is solved in a distributed fashion using the weighted Jacobi
iterative scheme (see Sec. II-G).

  * •

After the iterations converge, two failure criteria are checked: for the loss
of balance (see Sec. III-C) and overloading of connections (see Sec. III-D).

###  II-D Standard 3D frame model

In the adopted FE model, each module $p$ is represented by a node with 6
degrees of freedom ${\bf u}_{p}$ (3 displacements, $u_{x},u_{y},u_{z}$, and 3
rotations, $\tau_{x},\tau_{y},\tau_{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 ${\bf u}_{g}={\bf 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
${\bf K}_{pq}^{11}={\bf K}^{11}$ and ${\bf K}_{pq}^{12}={\bf K}^{12}$, where:

$\displaystyle{\bf
K}^{11}=\dfrac{E}{L^{3}}\begin{pmatrix}12I_{x}&0&0&0&-6I_{x}L&0\\\
0&12I_{y}&0&6I_{y}L&0&0\\\ 0&0&AL^{2}&0&0&0\\\ 0&6I_{y}L&0&4I_{y}L^{2}&0&0\\\
-6I_{x}L&0&0&0&4I_{x}L^{2}&0\\\ 0&0&0&0&0&J_{\nu}L^{2}\end{pmatrix},$ (1)
$\displaystyle{\bf
K}^{12}=\dfrac{E}{L^{3}}\begin{pmatrix}-12I_{x}&0&0&0&-6I_{x}L&0\\\
0&-12I_{y}&0&6I_{y}L&0&0\\\ 0&0&-AL^{2}&0&0&0\\\
0&-6I_{y}L&0&2I_{y}L^{2}&0&0\\\ 6I_{x}L&0&0&0&2I_{x}L^{2}&0\\\
0&0&0&0&0&-J_{\nu}L^{2}\end{pmatrix},$ (2)

while $E$, $L$, $A$, $I_{x}$, $I_{y}$ and $J_{\nu}$ 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 ${\bf K}_{pq}^{11}$ and
${\bf K}_{pq}^{12}$ take the form:

${\bf K}_{pq}^{11}={\bf\hat{R}}_{pq}{\bf
K}^{11}{\bf\hat{R}}_{pq}^{\text{T}}\quad,\quad{\bf
K}_{pq}^{12}={\bf\hat{R}}_{pq}{\bf K}^{12}{\bf\hat{R}}_{pq}^{\text{T}}\,,$ (3)

where

${\bf\hat{R}}_{pq}=\begin{pmatrix}{\bf R}_{pq}&{\bf 0}\\\ {\bf 0}&{\bf
R}_{pq}\end{pmatrix},$ (4)

${\bf R}_{pq}$ is the $3\times 3$ rotation matrix from CS’ to CS,
$\Box^{\text{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 ${\bf F}_{p}^{\text{ext}}=[0,0,\shortminus 9.81\cdot
M,0,0,0]^{\text{T}}$. The equilibrium equations are:

$\displaystyle\forall_{p}\left[\sum_{q}{\bf K}_{pq}^{11}{\bf u}_{p}+{\bf
K}_{pq}^{12}{\bf u}_{q}\right]={\bf F}_{p}^{\text{ext}}\,,$ (5)

with unknown vectors ${\bf u}_{p}$ and ${\bf u}_{q}$. After assembling the
global vectors ${\bf u}$, ${\bf F}^{\text{ext}}$ and stiffness matrix ${\bf
K}$, incorporating all degrees of freedom of the structure, Eq. (5) takes the
form ${\bf K}{\bf u}={\bf F}^{\text{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.

symbol | value |  description
---|---|---
$E$ | $100\text{\,}\mathrm{MPa}$ |  elastic modulus
$L,\quad A=L^{2}$ |  $40\text{\,}\mathrm{mm}$, $40\times{}40~{}$mm2 |  length & cross-sectional area
$I_{x}$, $I_{y}$ |  $L^{4}/12~{}$mm4 |  $x$ and $y$ moments of inertia
$J_{\nu}=\frac{J}{2(1+\nu)}$ |  $2.25\cdot L^{4}/41.6~{}$mm4 |  scaled $z$ torsion constant ($J$: $z$ torsion constant, $\nu$: Poisson’s ratio)
$M$ | $0.061\,06\text{\,}\mathrm{kg}$ |  mass of a block
$F^{\text{max}}_{\text{V}}$ | $11.98\text{\,}\mathrm{N}$ |  strength of a vertical connection
$F^{\text{max}}_{\text{L}}$ | $14.97\text{\,}\mathrm{N}$ |  strength of a lateral connection
TABLE I: Blinky Blocks’ geometric and material parameters

###  II-E Adding virtual modules

To predict the state of the system one reconfiguration step ahead, the
algorithm must take into account _virtual modules_ —the 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.

###  II-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 co-planar 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).

(a)  Figure 2: Regularized contact conditions: (a) a module in contact with
the ground, (b) contact-separation condition for the $f_{z}$ component, and
(c) tilting (stable/unstable) condition for the $m_{y}$ component. The dashed
lines indicate the respective “exact” (non-regularized) relationships for the
case of a rigid module in contact with rigid ground.

Without loss of generality, let us consider a module supported from below. The
module is loaded with the forces and torques ${\bf
F}=[f_{x},f_{y},f_{z},m_{x},m_{y},m_{z}]^{\text{T}}$, which correspond to the
kinematic variables ${\bf
u}=[u_{x},u_{y},u_{z},\tau_{x},\tau_{y},\tau_{z}]^{\text{T}}$. The contact
condition (see Fig. 2b) takes the form of the Signorini problem:

$f_{z}\leq 0\qquad\&\qquad u_{z}\geq 0\qquad\&\qquad f_{z}u_{z}=0.$ (6)

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 ($\tau_{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:

$\displaystyle\Phi_{x}\leq
0\quad\&\quad\tau_{x}\cdot\mathrm{sign}(\hat{m}_{x})\geq
0\quad\&\quad\Phi_{x}\cdot\tau_{x}=0\,,$
$\displaystyle\mathrm{where}\quad\hat{m}_{x}=m_{x}-f_{y}\cdot
L/2\,,\quad\Phi_{x}=|\hat{m}_{x}|+f_{z}\cdot L/2\,;$
$\displaystyle\Phi_{y}\leq
0\quad\&\quad\tau_{y}\cdot\mathrm{sign}(\hat{m}_{y})\geq
0\quad\&\quad\Phi_{y}\cdot\tau_{y}=0\,,$
$\displaystyle\mathrm{where}\quad\hat{m}_{y}=m_{y}+f_{x}\cdot
L/2\,,\quad\Phi_{y}=|\hat{m}_{y}|+f_{z}\cdot L/2\,;$ (8)

where $L$ is the module’s size. The tilting (bending) condition expresses the
fact that the torques $|\hat{m}_{x}|$ or $|\hat{m}_{y}|$ cannot exceed the
torque $\shortminus f_{z}\cdot L/2$ produced by the compressive force
$\shortminus 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–II-F).

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$):

${\bf
F}_{pq}^{\mathrm{tr}}=[f_{x}^{\mathrm{tr}},f_{y}^{\mathrm{tr}},f_{z}^{\mathrm{tr}},m_{x}^{\mathrm{tr}},m_{y}^{\mathrm{tr}},m_{z}^{\mathrm{tr}}]^{\text{T}}={\bf
K}^{11}_{pq}\bar{{\bf u}}_{p}\,,$ (9)

where the term ${\bf K}^{12}_{pq}\bar{{\bf u}}_{q}$ is absent because the
support is assumed to be immobile, $\bar{{\bf u}}_{q}={\bf 0}$. Then, a
corrected vector ${\bf F}_{pq}$ is determined, taking into account two
conditions:
_(i) the unilateral (normal) contact-separation condition combined with no-
sliding and no-twisting conditions_ :

$\begin{cases}\begin{array}[]{llr}{\bf F}_{pq}:=\gamma\cdot{\bf
F}_{pq}^{\mathrm{tr}}&\textrm{ for }f_{z}^{\mathrm{tr}}\geq 0&\textrm{
(separation)}\\\
(f_{x},f_{y},f_{z},m_{z})=(f_{x}^{\mathrm{tr}},f_{y}^{\mathrm{tr}},f_{z}^{\mathrm{tr}},m_{z}^{\mathrm{tr}})&\textrm{
otherwise}&\textrm{ (contact)}\end{array}\end{cases}$ (10)

_(ii) the tilting (bending) condition, used only when $f_{z}^{\mathrm{tr}}<0$,
i.e. when the module is in contact:_

$m_{x}:=\begin{cases}\begin{array}[]{lll}m_{x}^{\mathrm{tr}}&\hskip
1.99997pt\textrm{for }\Phi_{x}^{\mathrm{tr}}<0&\text{(stable)}\\\ \gamma\cdot
m_{x}^{\mathrm{tr}}+\bar{\gamma}\cdot\frac{L}{2}\left[f_{y}^{\mathrm{tr}}-\mathrm{sign}\left(\hat{m}_{x}^{\mathrm{tr}}\right)f_{z}^{\mathrm{tr}}\right]&\hskip
1.99997pt\textrm{otherwise}&\text{(tilting)}\end{array}\end{cases}$ (11)
$m_{y}:=\begin{cases}\begin{array}[]{lll}m_{y}^{\mathrm{tr}}&\hskip
1.99997pt\textrm{for }\Phi_{y}^{\mathrm{tr}}<0&\text{(stable)}\\\ \gamma\cdot
m_{y}^{\mathrm{tr}}-\bar{\gamma}\cdot\frac{L}{2}\left[f_{x}^{\mathrm{tr}}+\mathrm{sign}\left(\hat{m}_{y}^{\mathrm{tr}}\right)f_{z}^{\mathrm{tr}}\right]&\hskip
1.99997pt\textrm{otherwise}&\text{(tilting)}\end{array}\end{cases}$ (12)

where $\hat{m}_{x}^{\mathrm{tr}}$, $\hat{m}_{y}^{\mathrm{tr}}$,
$\Phi_{x}^{\mathrm{tr}}$ and $\Phi_{y}^{\mathrm{tr}}$ are computed as in Eqs.
(II-F) and (II-F), but using the components of ${\bf F}_{pq}^{\mathrm{tr}}$;
$\gamma=10^{-4}$ and $\bar{\gamma}=1-\gamma$.

The corrected contact tangent matrix $\bar{{\bf K}}{}^{11}_{pq}$ is obtained
as a derivative of ${\bf F}_{pq}$ with respect to $\bar{{\bf u}}_{p}$. Again,
the matrix $\bar{{\bf K}}{}^{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 $\gamma$, 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.

###  II-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 predictor-corrector scheme (Sec. II-F), reads:

$\displaystyle\bar{{\bf K}}\bar{{\bf u}}=\bar{{\bf F}}^{\text{ext}}.$ (13)

Eq. (13) is solved iteratively using the weighted Jacobi scheme [16]. A single
iteration $i\rightarrow i+1$ for module $p$ reads:

$\displaystyle\bar{{\bf u}}{}_{p}^{i+1}=\beta\bar{{\bf
D}}_{p}{}^{\\!\\!-1}\left(\bar{{\bf F}}_{p}^{\text{ext}}-\bar{{\bf
R}}_{p}\bar{{\bf u}}{}_{p}^{i}-\sum_{q}\bar{{\bf K}}_{pq}^{12}\bar{{\bf
u}}{}_{q}^{i}\right)+(1-\beta)\bar{{\bf u}}{}_{p}^{i}\,,\qquad$ (14)

where $\bar{{\bf D}}_{p}=\text{diag}\left(\sum_{q}\bar{{\bf
K}}_{pq}^{11}\right)$ and $\bar{{\bf R}}_{p}=\left(\sum_{q}\bar{{\bf
K}}_{pq}^{11}\right)-\bar{{\bf D}}_{p}$ are the diagonal and the remainder
parts of the respective stiffness sub-matrices, while $\beta=2/3$. Initially,
we set $\bar{{\bf u}}{}^{0}={\bf 0}$ (alternatively, the solution for the non-
perturbed state, if available, could be used instead of ${\bf 0}$, which would
reduce the necessary number of iterations). Note that in the iteration $i+1$
only the values of $\bar{{\bf u}}$ 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 $\bar{{\bf u}}_{p}^{0}={\bf
0}$ to all its neighbors and initializes its iteration counter,
$\textrm{iter}_{p}=0$. In a given iteration $\textrm{iter}_{p}=i$, $B_{p}$ can
receive from any of its neighbors $B_{q}$ a message containing $\bar{{\bf
u}}_{q}^{i}$ (displacements of $B_{q}$ calculated in iteration $i$), which is
then stored in $B_{p}$’s buffer. When $B_{p}$ has received $\bar{{\bf
u}}_{q}^{i}$ from all its neighbors, it computes $\bar{{\bf u}}_{p}^{i+1}$
(see Eq. 14), increments its counter,
$\textrm{iter}_{p}\leftarrow\textrm{iter}_{p}+1$, and sends $\bar{{\bf
u}}_{p}^{i+1}$ 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].

###  II-H Convergence properties and possible improvements

The Weighted Jacobi scheme converges if the spectral radius of the iteration
matrix ${\bf C}_{\beta}={\bf I}-\beta\bar{{\bf D}}^{-1}\bar{{\bf K}}$ is less
than 1:

$\displaystyle\rho({\bf
C}_{\beta})=\text{max}(|\lambda_{1}|,\ldots,|\lambda_{N}|)<1,$ (15)

where $\lambda_{i}$ are the eigenvalues of ${\bf C}_{\beta}$. 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.

| $\mathbb{C}_{\text{C}}$ |  $\mathbb{C}_{\text{T}}$, $\mathbb{C}_{\text{M}}$, $\mathbb{C}_{\text{SST}}$, $\mathbb{C}_{\text{MST}}$ | $\mathbb{C}_{\text{WJ}}$
---|---|---|---
Execution time | $O(d)$ | $O(\tilde{d})$ | $O(n^{2})$
No. of CPU operat./module | $O(1)$ | $O(1)$ | $O(n^{2})$
Total no. of messages sent | $O(n)$ | $O(n)$ | $O(n^{3})$
Memory usage per module | $O(1)$ | $O(1)$ | $O(1)$
TABLE II: Complexity assessments for subroutines of the algorithm.
$\mathbb{C}_{\text{C}}$ refers to centroid selection,
$\mathbb{C}_{\text{T}}$—to tree construction, $\mathbb{C}_{\text{M}}$—to
finding the center of mass, $\mathbb{C}_{\text{SST}}$—to the simplified
stability check, $\mathbb{C}_{\text{MST}}$—to the model-based stability check,
and $\mathbb{C}_{\text{WJ}}$—to the weighted Jacobi procedure.

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 stability—without
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, and $\tilde{d}$ is the depth of the
spanning tree (usually, $\tilde{d}$ and $d$ are of the same order; see Sec.
III-A for more details).

Figure 3: Flowchart of the algorithm.  (a)  Figure 4: (a) A modular robot
(unstable) with a selected centroid (thick blue line), the center of mass
(star) and a spanning tree (arrows). (b) Detection of whether a given point
(star) is inside the convex hull of the support points. The condition is
fulfilled if and only if all the straight angles (shaded) sum up to the full
angle (they do not in this example).

###  III-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 quasi-
optimal 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.

###  III-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.

  1. (a)

The center of mass of the robot is computed. Starting at the leaves of the
spanning tree, each node sums up the masses of all its subbranches and its
own, $m_{i}$, and likewise the weighted centers of mass of all its subbranches
and its own, $m_{i}{\bf X}_{i}$. The two sums are propagated to the parent
node and the process continues. At the centroid, the center of mass of the
robot is retrieved as $[X,Y,Z]=(\Sigma m_{i}{\bf X}_{i})/(\Sigma m_{i})$.

  2. (b)

$X$ and $Y$ are broadcast over the tree.

  3. (c)

For any supporting module $i$, its _safe angle range_ $[\alpha_{i},\beta_{i}]$
is determined as the sum of safe angle ranges of its corners. The
$180^{\circ}$ safe angle range of a corner ${\bf p}_{j}=[X_{j},Y_{j},0]$ is
swept by the planar vector $[X_{j}-X,Y_{j}-Y]$ when turning $90^{\circ}$ left
or right. It covers those directions in which the structure cannot tilt; see
Fig. 4b. The safe angle range of any corner ${\bf
p}_{j}=[X_{j},Y_{j},Z_{j}\neq 0]$ is assumed to be empty.

  4. (d)

The safe angle ranges are summed up over the tree, just like masses were in
step (a). The summation always gives a single interval, because all considered
ranges are either empty or not less than the straight angle.

  5. (e)

The structure is stable if and only if the aggregated angle range at the
centroid equals the full angle.

###  III-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 elasto-static
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:

  1. (a)

Solve the mechanical problem (see Sec. II) and initiate the stability check
(the spanning tree is used again).

  2. (b)

For every module, determine the set of its active corner points by checking
the contact conditions (Sec. II-F):

     * •

No contact $\to$ return the empty set.

     * •

Tilting in two directions $\to$ return a single corner—the common point of the
two edges of rotation.

     * •

Tilting in one direction $\to$ return two corners—the end points of the edge
of rotation.

     * •

Contact & no tilting $\to$ return a special ‘stable’ state.

  3. (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’.

     * •

If the points in $S$ are noncollinear then set $S=\text{`stable'}$.

     * •

If multiple points in $S$ are collinear then leave only two.

     * •

Pass $S$ up the spanning tree.

  4. (d)

The stability check ends at the centroid, with the result being either
‘stable’ or not.

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.

###  III-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:

$\displaystyle[f_{x},f_{y},f_{z},m_{x},m_{y},m_{z}]^{\text{T}}=\frac{1}{2}{\bf\hat{R}}_{pq}^{\text{T}}({\bf
F}_{pq}-{\bf F}_{qp})=$
$\displaystyle\qquad=\frac{1}{2}{\bf\hat{R}}_{pq}^{\text{T}}(\bar{{\bf
K}}^{11}_{pq}\bar{{\bf u}}_{p}+\bar{{\bf K}}^{12}_{pq}\bar{{\bf
u}}_{q}-\bar{{\bf K}}^{11}_{qp}\bar{{\bf u}}_{q}-\bar{{\bf
K}}^{12}_{qp}\bar{{\bf u}}_{p})\,,\qquad$ (16)

where ${\bf\hat{R}}_{pq}^{\text{T}}$ 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^{\textrm{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^{\textrm{max}}$ can take two values; see
Table I. The safety condition for both tension and bending is:

$F^{\textrm{max}}>2\cdot\max(|m_{x}|,|m_{y}|)/L+f_{z}\,.$ (17)

The check is performed for all connections and the results are aggregated by
the centroid over the spanning tree.

##  IV Implementation, simulations and experiments

###  IV-A Implementation details

The procedures have been implemented and tested in the integrated environment
developed at FEMTO-ST111Programmable Matter project at FEMTO-ST:
https://projects.femto-st.fr/programmable-matter/.. 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\text{\,}\mathrm{B}\mathrm{y}\mathrm{t}\mathrm{e}\mathrm{s}$ was used and
the maximum message size was set to
$17\text{\,}\mathrm{B}\mathrm{y}\mathrm{t}\mathrm{e}\mathrm{s}$ (messages
containing $6\times 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 Lego-like system
reinforced with a central magnet, and the latter by 4 magnets placed in the
corners of each face; see also Fig. 5.

(a)  Figure 5: Blinky Blocks: functional ones (left and right), and two pieces
of a dismantled one with a top view of the motherboard (in the middle).

###  IV-B Simulations and experiments

(a)  (b)

Figure 6: Experiments in VisibleSim and on Blinky Blocks. Computation times
and iteration numbers are shown in insets.

Six different modular configurations and failure scenarios were investigated
(see Fig. 6), both in VisibleSim and on Blinky Blocks 222Videos of selected
experiments: https://youtu.be/d3aE8GjbYd8. 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 model-
based 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$\to$b) is designed
to be mechanically safe, while the second one (b$\to$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.

(a)

(a)

(b)

(b)

Figure 7: Number of weighted Jacobi iterations necessary to attain a given
accuracy for fixed-arms of different sizes. Size 10 refers to the example in
Fig. 6-#4(a); other sizes have shorter/longer arms. (a) Convergence
characteristics. (b) A quadratic polynomial fit for two accuracy levels.

_CPU time and convergence properties._ Computing $\bar{{\bf u}}^{i}_{p}$ and
exchanging messages with neighbours takes a Blinky Block a nearly constant
time $T\approx$$110.5\text{\,}\mathrm{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 $\|\bar{{\bf
u}}^{i}-\bar{{\bf u}}^{*}\|/\|\bar{{\bf u}}^{*}\|$ as the number of iterations
$i$ grows, where $\bar{{\bf u}}^{*}$ 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.

_Acknowledgement._ This work was partially supported by the EU Horizon 2020
Marie Sklodowska Curie Individual Fellowship _MOrPhEM_ (grant No. 800150), by
the NCN Project “Micromechanics of Programmable Matter” (grant No.
2011/03/D/ST8/04089), by the ANR (ANR-16-CE33-0022-02), the French
Investissements d’Avenir program, ISITE-BFC project (ANR-15-IDEX-03), EIPHI
Graduate School (contract ANR-17-EURE-0002), Mobilitech project and EU Horizon
2020 research and innovation programme (grant No. 811099 TWINNING Project
DRIVEN for the Univ. of Luxembourg).

## References

  * [1] M. Yim, W. Shen, B. Salemi, D. Rus, M. Moll, H. Lipson, E. Klavins, and G. S. Chirikjian, “Modular self-reconfigurable robot systems,” IEEE Robotics Automation Magazine, vol. 14, no. 1, pp. 43–52, 2007.
  * [2] N. J. Mlot, C. A. Tovey, and D. L. Hu, “Fire ants self-assemble into waterproof rafts to survive floods,” Proceedings of the National Academy of Sciences, vol. 108, no. 19, pp. 7669–7673, 2011.
  * [3] S. C. Goldstein, J. D. Campbell, and T. C. Mowry, “Programmable matter,” IEEE Computer, vol. 38, no. 6, pp. 99–101, 2005.
  * [4] A. Nguyen, L. J. Guibas, and M. Yim, “Controlled module density helps reconfiguration planning,” in Proceedings of the 4th International Workshop on Algorithmic Foundations of Robotics, pp. 23–36, 2000.
  * [5] R. Fitch and Z. Butler, “Million module march: Scalable locomotion for large self-reconfiguring robots,” The International Journal of Robotics Research, vol. 27, no. 3-4, pp. 331–343, 2008.
  * [6] K. Støy, “Using cellular automata and gradients to control self-reconfiguration,” Robotics and Autonomous Systems, vol. 54, no. 2, pp. 135–141, 2006.
  * [7] B. Piranda and J. Bourgeois, “A distributed algorithm for reconfiguration of lattice-based modular self-reconfigurable robots,” in Proceedings of the 24th Euromicro International Conference on Parallel, Distributed, and Network-Based Processing (PDP), pp. 1–9, IEEE, 2016.
  * [8] A. Naz, B. Piranda, J. Bourgeois, and S. C. Goldstein, “A distributed self-reconfiguration algorithm for cylindrical lattice-based modular robots,” in 15th IEEE International Symposium on Network Computing and Applications, pp. 254–263, 2016.
  * [9] J. Lengiewicz and P. Hołobut, “Efficient collective shape shifting and locomotion of massively-modular robotic structures,” Autonomous Robots, vol. 43, no. 1, pp. 97–122, 2019.
  * [10] P. Thalamy, B. Piranda, and J. Bourgeois, “A survey of autonomous self-reconfiguration methods for robot-based programmable matter,” Robotics and Autonomous Systems, vol. 120, p. 103242, 2019.
  * [11] P. J. White, S. Revzen, C. E. Thorne, and M. Yim, “A general stiffness model for programmable matter and modular robotic structures,” Robotica, vol. 29, pp. 103–121, 2011.
  * [12] J. Hiller and H. Lipson, “Dynamic simulation of soft multimaterial 3d-printed objects,” Soft robotics, vol. 1, no. 1, pp. 88–101, 2014.
  * [13] J. Campbell and P. Pillai, “Collective actuation,” International Journal of Robotics Research, vol. 27, no. 3–4, pp. 299–314, 2008.
  * [14] P. Hołobut, M. Kursa, and J. Lengiewicz, “Efficient modular-robotic structures to increase the force-to-weight ratio of scalable collective actuators,” in Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 3302–3307, 2015.
  * [15] J. Lengiewicz, M. Kursa, and P. Hołobut, “Modular-robotic structures for scalable collective actuation,” Robotica, vol. 35, pp. 787–808, 2017.
  * [16] P. Hołobut and J. Lengiewicz, “Distributed computation of forces in modular-robotic ensembles as part of reconfiguration planning,” in Proceedings of the IEEE International Conference on Robotics and Automation, pp. 2103–2109, 2017.
  * [17] B. Piranda, “Visiblesim: Your simulator for programmable matter,” in Algorithmic Foundations of Programmable Matter, Dagstuhl Seminar 16271, p. 12, 2016.
  * [18] B. T. Kirby, M. Ashley-Rollman, and S. C. Goldstein, “Blinky blocks: A physical ensemble programming platform,” in CHI ’11 Extended Abstracts on Human Factors in Computing Systems, CHI EA ’11, (New York, NY, USA), pp. 1111–1116, ACM, 2011.
  * [19] A. Tanenbaum, Structured Computer Organization.  Pearson Prentice Hall, 5th ed., 2006.
  * [20] M. Raynal, Distributed algorithms for message-passing systems, vol. 500.  Springer, 2013.
  * [21] G. Strang, Linear algebra and its applications.  Thomson, Brooks/Cole, 2006.
  * [22] Y. Saad, Iterative methods for sparse linear systems.  Society for Industrial and Applied Mathematics, 2nd ed., 2003.
  * [23] J. Lengiewicz, S. P. A. Bordas, and P. Hołobut, “Autonomous model-based assessment of mechanical failures of reconfigurable modular robots with a conjugate gradient solver,” in Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 11696–11702, 2020.
  * [24] W. L. Briggs, V. E. Henson, and S. F. McCormick, A Multigrid Tutorial.  SIAM, 2000.
  * [25] T. A. Wiesner, A. Popp, M. W. Gee, and W. A. Wall, “Algebraic multigrid methods for dual mortar finite element formulations in contact mechanics,” International Journal for Numerical Methods in Engineering, vol. 114, no. 4, pp. 399–430, 2018.
  * [26] P. Kerfriden, P. Gosselet, S. Adhikari, and S. P. A. Bordas, “Bridging proper orthogonal decomposition methods and augmented newton–krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems,” Computer Methods in Applied Mechanics and Engineering, vol. 200, pp. 850–866, 2011.
  * [27] L. Beex, P. Kerfriden, T. Rabczuk, and S. P. A. Bordas, “Quasicontinuum-based multiscale approaches for plate-like beam lattices experiencing in-plane and out-of-plane deformation,” Computer Methods in Applied Mechanics and Engineering, vol. 279, pp. 348–378, 2014.
  * [28] A. Naz, B. Piranda, S. C. Goldstein, and J. Bourgeois, “Approximate-centroid election in large-scale distributed embedded systems,” in 2016 IEEE 30th International Conference on Advanced Information Networking and Applications (AINA), pp. 548–556, 2016.
  * [29] B. Piranda and J. Bourgeois, “Designing a quasi-spherical module for a huge modular robot to create programmable matter,” Autonomous Robots, vol. 42, no. 8, pp. 1619–1633, 2018.


# A global method for mixed categorical optimization with catalogs

Charlie Vanaret<EMAIL_ADDRESS>Chair of Software and Algorithms
for Discrete Optimization
Technische Universität Berlin, 10623 Berlin, Germany

###### Abstract

In this article, we propose an algorithmic framework for globally solving
mixed problems with continuous variables and categorical variables whose
properties are available from a catalog. It supports catalogs of arbitrary
size and properties of arbitrary dimension, and does not require any modeling
effort from the user. Our tree search approach, similar to spatial branch and
bound methods, performs an exhaustive exploration of the range of the
properties of the categorical variables ; branching, constraint programming
and catalog lookup phases alternate to discard inconsistent values. A novel
catalog-based contractor guarantees consistency between the categorical
properties and the existing catalog items. This results in an intuitive
generic approach that is exact and easy to implement. We demonstrate the
validity of the approach on a numerical example in which a categorical
variable is described by a two-dimensional property space.

###### keywords:

Global optimization , Categorical variables , Branch and bound , Constraint
programming , Catalog lookup

††journal: European Journal of Operational Research

##  1 Motivation

Categorical variables – discrete variables that take values from a finite set
of non-numerical values and thus have no natural ordering – arise in numerous
structural optimization applications: optimal truss design [16], design of
optical systems [2], tyre selection [21], aircraft design [5], product family
design [19], optimization of thermal insulation systems [1], etc. While
problems with discrete variables are noticeably harder to solve than
continuous problems, categorical variables add another level of complexity,
since the core technique for solving discrete optimization problems
(relaxation of the integrity constraints) does not hold.

###  1.1 Problem definition

We start by making two assumptions:

  * A.1

all catalog items are represented by the same set of properties ;

  * A.2

catalog items are not mutually exclusive (that is, there are no forbidden
combinations of items) ;

and consider the following mixed categorical optimization problem:

$\displaystyle\underset{x,y}{\text{min}}$ $\displaystyle f(x,y)$ (1) s.t.
$\displaystyle\underline{c}\leq c(x,y)\leq\overline{c}$ $\displaystyle u(y)\in
C$ $\displaystyle x\in[\underline{x},\overline{x}],\quad
y\in[\underline{y},\overline{y}]$

where

  * 1.

$x\in[\underline{x},\overline{x}]$ is a vector of continuous variables ;

  * 2.

$u(y)$ is a categorical variable described by a vector of properties
$y\in[\underline{y},\overline{y}]$ ;

  * 3.

$f$ is the (possibly nonlinear) objective function ;

  * 4.

$\underline{c}\leq c(x,y)\leq\overline{c}$ are (possibly nonlinear)
constraints ;

  * 5.

$u(y)\in C$ is a catalog constraint: the categorical variable $u$ described by
its properties $y$ should be assigned an item from the finite catalog $C$.

For the sake of clarity and without loss of generality, we describe the
problem with a single categorical variable.

###  1.2 Literature on methodology

Evolutionary and genetic algorithms have been among the first techniques to
automatize catalog-based optimization in industrial applications. While some
authors adopted a generic approach [9, 10], others tailored the encoding (a
hierarchical structure), crossover (a trial-and-error selection of
breakpoints) and mutation operators to catalog design [11].

Recently, modern approaches based on exhaustive search tree methods and model
reformulation have proven efficient. [1] described mixed-integer nonlinear
(MINLP) models by replacing the assignments of categorical variables with
binary variables and used standard branch and bound methods via the AMPL
modeling language. [19] combined a local search and a global search for
categorical optimization: a pattern search technique that iteratively improves
the current point by visiting its categorical neighborhood is embedded within
a discrete global descent framework. [5] presented a branch-and-bound approach
that instantiates a categorical variable at each level of the search tree. It
evaluates lower bounds by relaxing the constraints that depend on undecided
categorical variables and progressively inserts them back into the problem
when the corresponding categorical variables are assigned a value. In [16],
the authors interpreted the categorical properties as points in a high-
dimensional metric space. They performed dimensionality reduction by using a
manifold learning technique that preserves geodesic distances. They then
defined new operators for genetic algorithms based on shortest paths in the
resulting low-dimensional graph. [21] described a splitting strategy for
simulation-based optimization that exploits convex relaxations of the discrete
search space to evaluate approximate lower bounds. The bounds are then used in
a branching heuristic within a multilevel coordinate search.

###  1.3 Contributions and outline

The key idea of our approach is to notice that, while $u$ is a categorical
variable, its properties $y$ are numerical variables. Instead of optimizing
the categorical variables $u$ directly, we suggest to treat the underlying
properties $y$ as optimization variables. We propose a novel tree search
approach for mixed categorical optimization, similar to spatial branch and
bound methods [23], that manipulates continuous relaxations of the properties
$y$ and alternates between branching and constraint programming phases on both
the continuous variables $x$ and the categorical properties $y$.

To ensure consistency between the optimization variables $y$ and the existing
catalog entries, we introduce a new filtering operator (or contractor [12]),
Clutch (Catalog LookUp Then Convex Hull), that performs a catalog lookup and
reduces the current range of $y$ by considering the union of the existing
catalog items within this range.

We highlight the main advantages of our strategy:

  1. 1.

it is a generic approach that can be applied to a wide variety of applications
;

  2. 2.

it is exact and produces the global optimum within a given numerical tolerance
;

  3. 3.

since we directly optimize the properties of the categorical variables, no
modeling effort (e.g. reformulation with binary variables [1]) is required
from the user ;

  4. 4.

it is intuitive and easy to implement ;

  5. 5.

the exploration of the property space is not limited to branching and bounding
phases: filtering algorithms (contractors) based on relational logic may be
invoked to contract the ranges of the property variables ;

  6. 6.

gradients can be exploited to compute enclosures (e.g. linear under- and
overestimators) of the objective function and constraints. Note however that,
on account of the discrete nature of categorical optimization, first-order
optimality conditions do not hold at the solution ;

  7. 7.

constraint programming techniques, if implemented using interval arithmetic,
are robust to roundoff errors.

Section 2 describes our algorithmic framework and the Clutch contractor. In
Section 3, we present an example with a continuous variable and a categorical
variable with a two-dimensional property space, and we demonstrate the
validity of our global method on two scenarios.

##  2 A framework for mixed categorical global optimization with catalogs

###  2.1 Description of the framework

We describe the overall framework in Algorithm 1 ; its tree search structure
is similar to spatial branch and bound solvers. It essentially transforms the
global optimization problem into a numerical constraint satisfaction problem
that keeps track of the best known upper bound $\tilde{f}$ of the global
minimum $f^{*}$. The constraint $f\leq\tilde{f}$ is enforced on each subspace
and is dynamically updated whenever $\tilde{f}$ is improved. If a solution is
sought within a tolerance $\varepsilon>0$ of the global minimum, the dynamic
constraint can be written $f\leq\tilde{f}-\varepsilon$.

Input: box $\bm{X}\times\bm{Y}$, objective function $f$, constraints, catalog,
tolerance $\varepsilon$

Output: upper bound $\tilde{f}$ of the global minimum $f^{*}$ such that
$\tilde{f}-f^{*}\leq\varepsilon$

Set upper bound $\tilde{f}\leftarrow+\infty$

Initialize queue $\mathcal{L}\leftarrow\\{(\bm{X}\times\bm{Y},-\infty)\\}$

while  _$\mathcal{L}$ not empty_ do

Exploration strategy: extract $(\bm{X}_{k}\times\bm{Y}_{k},lb_{k-1})$ from
$\mathcal{L}$

if  _$\tilde{f}-\varepsilon <lb_{k-1}$_ then discard the subspace

Lower bounding: compute a lower bound $lb_{k}$ of the problem

if  _$\tilde{f}-\varepsilon <lb_{k}$_ then discard the subspace

Filtering: contract $\bm{X}_{k}\times\bm{Y}_{k}$ with respect to the
constraints

Clutch: perform catalog lookup and contract $\bm{Y}_{k}$

if  _$\bm{X}_{k}\times\bm{Y}_{k}$ is empty_ then discard the subspace

Upper bounding: improve the best known upper bound $\tilde{f}$

Branching: split $\bm{X}_{k}\times\bm{Y}_{k}$ in two subspaces
$(\bm{X}_{k})_{i}\times(\bm{Y}_{k})_{i},i\in\\{1,2\\}$

for  _$i\in\\{1,2\\}$_ do

Clutch: perform catalog lookup and contract $(\bm{Y}_{k})_{i}$

Insert $((\bm{X}_{k})_{i}\times(\bm{Y}_{k})_{i},lb_{k})$ into $\mathcal{L}$

return $\tilde{f}$

Algorithm 1 Framework for mixed categorical optimization

Our framework invokes a sequence of generic strategies described in Section
2.2: exploration of the search space, lower bounding, filtering of the current
subspace, upper bounding and branching. Lower bounding and filtering are
performed before upper bounding as they rely on refutation: the current
subspace should be discarded as soon as it is proven suboptimal or infeasible.
The catalog-based contractor Clutch is described in Section 2.3.

###  2.2 Generic strategies

Our framework is described in a generic fashion. In this section, we describe
possible strategies for exploration, lower bounding, filtering, upper bounding
and branching.

####  2.2.1 Exploration strategy

The priority queue $\mathcal{L}$ is usually implemented as a binary heap, that
is a complete binary tree in which the priority of a node is greater than the
priority of its children. Inserting an element or extracting the element with
the highest priority is carried out in logarithmic time. The choice of
priority determines the exploration strategy of the search space (depth-first
search, best-first search, breadth-first search, MaxDist [26], etc.)

####  2.2.2 Lower bounding

The most commonly used technique to evaluate a lower bound of the original
problem on a given subspace is the Reformulation Linearization Technique (RLT)
[24]. It consists in generating a convex outer approximation (by linearization
or convexification) of the problem, then solving the convex relaxed problem.
Its interval counterpart X-Newton [4] implements reliable linearizations using
interval arithmetic [25] and a cheap duality-based postprocessing step [22]
for evaluating rigorous lower bounds even in the presence of roundoff errors.

A much cheaper approach, albeit of coarser quality, is to compute an enclosure
of the range of the sole objective function on a given subspace using interval
arithmetic. It amounts to computing a lower bound of the unconstrained problem
and may therefore largely overestimate the actual range.

####  2.2.3 Filtering

Filtering algorithms (contractors) reduce the domains of the variables with
respect to individual constraints (local consistency) or all the constraints
simultaneously (global consistency) by discarding inconsistent values. They
usually stem from the mixed integer, constraint programming and numerical
analysis communities, and may exploit the syntax tree, range or monotonicity
of an individual constraint, or perform a linearization/convexification of the
whole system of constraints.

22$[2,2]$$[0,20]$$[0,40]$$[-10,10]$$[0,100]$$[0,16]$$[-100,16]$$z$$x$$y$$\times$$\wedge$$-$$=$2
(a) Bottom-up evaluation phase

2$x$$z$$y$2$[2,2]$Step 3Step 1Step 2Step 4Step 5Step
6$\wedge$$-$$=$$\times$$\mathbf{[0,8]}$$\mathbf{[0,16]}$$[0,40]$2$\mathbf{[-4,4]}$$[-10,10]$$\mathbf{[0,16]}$$[-100,16]$$[0,20]$$\mathbf{[0,16]}$$[0,16]$$[0,100]$$\mathbf{[0,16]}$
(b) Top-down propagation phase

Figure 1: FBBT/HC4Revise on the constraint $2x=z-y^{2}$

One of the most intuitive approaches, called evaluation-propagation algorithm
[20], FBBT [6] or HC4Revise [7], draws its inspiration from relational
interval arithmetic [13]. It carries out two traversals of the syntax tree of
a constraint: the bottom-up phase evaluates each node (an elementary
operation) using interval arithmetic, and the top-down phase propagates the
constraint downwards by using inverse (or projection) elementary functions.
Figure 1 illustrates the approach on the constraint $2x=z-y^{2}$: the
procedure contracts the initial range $[0,20]\times[-10,10]\times[0,16]$ to
$[0,8]\times[-4,4]\times[0,16]$ by removing values that cannot satisfy the
constraint.

In order to contract a box with respect to a system of constraints, a fixed-
point algorithm (an idempotent propagation loop) handles the constraints
individually and propagates range reductions through the system (Algorithm 2).

Input: box $\bm{X}$, constraints $\\{c_{1},\ldots,c_{m}\\}$

Output: contracted box $\bm{X}$

Initialize the set of awake constraints
$\mathcal{A}\leftarrow\\{c_{1},\ldots,c_{m}\\}$

while  _$\mathcal{A}$ is not empty_ do

Extract a constraint $c_{i}$ from $\mathcal{A}$

$\bm{X}^{\prime}\leftarrow$ contraction of $\bm{X}$ with respect to $c_{i}$

if  _$\bm{X}^{\prime}\neq\bm{X}$_ then

Add to $\mathcal{A}$ the constraints that involve contracted variables of
$\bm{X}^{\prime}$

$\bm{X}\leftarrow\bm{X}^{\prime}$

end if

end while

return $\bm{X}$

Algorithm 2 Fixed-point algorithm

Convexification techniques, called OBBT [6] or X-Newton [4], generate linear
under- and overestimators of all the constraints on a given subspace. The
lower (resp. upper) bound of the variable $x_{i}$ is reduced by minimizing
(resp. maximizing) $x_{i}$. Overall, the ranges of $n$ variables may be
reduced with $2n$ calls to a linear programming solver. For a comprehensive
survey, see [26].

####  2.2.4 Upper bounding

The objective value of any feasible point of the problem (provided that the
feasible set is not empty) is an upper bound of the global minimum. The best
known upper bound $\tilde{f}$ of the global minimum can be improved by local
search algorithms or metaheuristics in an integrated [27] or cooperative
manner [15, 8, 14, 3, 26]. A cheaper strategy consists in systematically
updating $\tilde{f}$ with the objective value of the midpoint of the current
subspace, provided that it is feasible [17].

####  2.2.5 Branching

Branching consists in partitioning the current subspace into several (usually
two) subspaces. There exist a variety of strategies, among which three stand
out:

  * 1.

the variable with the largest domain is partitioned ;

  * 2.

the variables are partitioned one after the other in a round-robin fashion ;

  * 3.

the smear heuristic picks the variable that has the largest estimated impact
(a measure based on the interval gradient) on the objective function [18].

###  2.3 Clutch catalog-based contractor

Clutch (Catalog LookUp Then Convex Hull) is a novel catalog-based contractor
that enforces local consistency with respect to the individual catalog
constraint $u(y)\in C$. Recall that in our framework, the range of the
categorical properties $y$ is reduced by the spatial branch and bound
approach. Clutch should therefore be invoked whenever the range of $y$ is
reduced (upon branching and filtering) to guarantee that the range of $y$ is
consistent with existing catalog items. This ensures that subspaces
$\bm{X}\times\bm{Y}$ that satisfy the standard constraints $c(x,y)$ or the
dynamic constraint $f(x,y)\leq\tilde{f}$ but cannot contain a catalog item
($u(y)\notin C$) are discarded by refutation.

More formally, given a box $\bm{Y}$ and a catalog $C$, Clutch returns a box
$\bm{Y}_{\mathit{clutch}}$ that satisfies the correction property:

$\bm{Y}\,\cap C\,\subset\bm{Y}_{\mathit{clutch}}\subset\bm{Y}$ (2)

where $\bm{Y}\,\cap\,C$ are the catalog items that are in $\bm{Y}$. Since our
framework manipulates compact boxes whose faces are parallel to the axes, the
best possible relaxation of $\bm{Y}\,\cap\,C$ is the smallest box that
contains all items of $\bm{Y}\,\cap\,C$:

$\bm{Y}_{\mathit{clutch}}=\square(\bm{Y}\cap C)$ (3)

where $\square$ is the convex hull operation. The procedure is described in
Algorithm 3.

Input: box $\bm{Y}$, catalog

Output: filtered box $\bm{Y}_{\mathit{clutch}}$

Initialize $\bm{Y}_{\mathit{clutch}}\leftarrow\varnothing$

foreach  _item $y$ in catalog_ do

if  _$y$ is in $\bm{Y}$_ then

Compute the convex hull
$\bm{Y}_{\mathit{clutch}}\leftarrow\square(\bm{Y}_{\mathit{clutch}},\\{y\\})$

return $\bm{Y}_{\mathit{clutch}}$

Algorithm 3 Clutch contractor

Upon completion of Clutch, three situations may occur:

  1. 1.

$\bm{Y}\,\cap\,C$ is empty, that is there are no catalog items that belong to
$\bm{Y}$: $\bm{Y}$ is inconsistent with the catalog constraint and is
discarded ;

  2. 2.

$\bm{Y}\,\cap\,C$ contains a single catalog item and is therefore equal to
$\square(\bm{Y}\,\cap\,C)$: the categorical variable is assigned to this item
;

  3. 3.

$\bm{Y}\,\cap\,C$ contains at least two distinct items:
$\square(\bm{Y}\,\cap\,C)$ is therefore a relaxation of $\bm{Y}\,\cap\,C$ and
further branching or filtering of the subspace is required.

The reduction of the range of the categorical properties $y$ can be propagated
to the other constraints $c(x,y)$ via the use of a fixed-point algorithm (see
Algorithm 2). We illustrate in Section 3 how constraint propagation involving
Clutch (for a catalog constraint) and HC4Revise (for a standard constraint)
can iteratively reduce the ranges of the categorical and continuous variables
on two scenarios.

##  3 Constraint propagation with CLUTCH: an example

We consider the following example:

$\displaystyle\underset{x,y}{\text{min}}$ $\displaystyle y_{1}^{3}$ (4) s.t.
$\displaystyle c_{1}(x,y):2y_{1}=x-y_{2}^{2}$ $\displaystyle c_{2}(y):u(y)\in
C$ $\displaystyle x\in[0,16],y_{1}\in[0,20],y_{2}\in[-10,10],$

where

  * 1.

$x$ is a continuous variable with range $[0,16]$ ;

  * 2.

the categorical variable $u$ is described by two properties $y_{1}$ and
$y_{2}$ with ranges $[0,20]$ and $[-10,10]$ ;

  * 3.

the catalog $C$ is described by pairs of $(y_{1},y_{2})$ properties ;

  * 4.

we seek the global solution with a tolerance of $\varepsilon=10^{-6}$.

We now introduce two scenarios, each with a different catalog. In scenario 1
(Section 3.1), the feasible set of the problem contains a single catalog item
; in this case, the assignment of $u(y)$ to a catalog item does not require
the use of the objective function. In scenario 2 (Section 3.2), the feasible
set of the problem contains two catalog items ; in this case, upper bounding
and the use of the dynamic constraint $f\leq\tilde{f}-\varepsilon$ is required
to discard one the two items.

###  3.1 Scenario 1: feasible set with a single catalog item

The catalog (Table 1) contains five items. The catalog items in their initial
domain $[0,20]\times[-10,10]$ are represented in Figure 2(a).

Table 1: Scenario 1: catalog $C$ Item | Property $y_{1}$ | Property $y_{2}$
---|---|---
1 | 4 | -8
2 | 3 | 2
3 | 7 | -3
4 | 14 | 8
5 | 19 | -8

Table 2 describes the successive filtering and branching phases (second
column) that reduce the ranges of $x$ and $y$ starting from
$(x,y_{1},y_{2})\in[0,16]\times[0,20]\times[-10,10]$: the third column
contains for each phase the current ranges of $x$ and $y$, while the fourth
column presents the resulting contracted ranges. The first column indicates
the corresponding subfigures in Figure 2. The branching phase and the two
resulting subspaces are highlighted in red.

Table 2: Scenario 1: alternance of filtering and branching phases | Phase | Current box | Result
---|---|---|---
(b) |  HC4Revise($c_{1}$) | $[0,16]\times[0,20]\times[-10,10]$ | $[0,16]\times[0,8]\times[-4,4]$
(c) |  Clutch($c_{2}$) | $[0,16]\times[0,8]\times[-4,4]$ | $[0,16]\times[3,7]\times[-3,2]$
|  HC4Revise($c_{1}$) | $[0,16]\times[3,7]\times[-3,2]$ | $[6,16]\times[3,7]\times[-3,2]$
(d) |  Branching on $y_{1}$ | $[6,16]\times[3,7]\times[-3,2]$ | $[6,16]\times{\color[rgb]{1,0,0}[3,5]}\times[-3,2]~{}\cup$
|  |  | $[6,16]\times{\color[rgb]{1,0,0}[5,7]}\times[-3,2]$
|  Clutch($c_{2}$) | $[6,16]\times{\color[rgb]{1,0,0}[5,7]}\times[-3,2]$ | $[6,16]\times[7,7]\times[-3,-3]$
|  HC4Revise($c_{1}$) | $[6,16]\times[7,7]\times[-3,-3]$ | $\varnothing$
(e) |  Clutch($c_{2}$) | $[6,16]\times{\color[rgb]{1,0,0}[3,5]}\times[-3,2]$ | $[6,16]\times[3,3]\times[2,2]$
(f) |  HC4Revise($c_{1}$) | $[6,16]\times[3,3]\times[2,2]$ | $\mathbf{[10,10]\times[3,3]\times[2,2]}$

In Figure 2, the blue domain is the feasible set of $c_{1}$ on the current box
$X\times Y_{1}\times Y_{2}$ (that is the set
$\\{2y_{1}=x-y_{2}^{2}~{}|~{}(x,y_{1},y_{2})\in X\times Y_{1}\times Y_{2}\\}$)
projected on the $(y_{1},y_{2})$ plane, the current box is represented as a
black rectangle, boxes contracted by HC4Revise or Clutch as red rectangles,
boxes partitioned during branching as green rectangles and discarded catalog
items as gray points.

(a) Initial range of properties $y_{1}$ and $y_{2}$

(b) HC4Revise with respect to $c_{1}$

(c) Clutch with respect to $c_{2}$

(d) Branching on $y_{1}$

(e) Clutch with respect to $c_{2}$

(f) HC4Revise with respect to $c_{1}$

Figure 2: Scenario 1: alternance of filtering and branching phases

These results show that the first contraction step discards three of the five
catalog items ; convergence to the optimal (and only feasible) point
$(x^{*},y_{1}^{*},y_{2}^{*})=(10,3,2)$ is achieved using constraint
propagation between HC4Revise (for $c_{1}$) and Clutch (for $c_{2}$) and a
single branching step on $y_{1}$. The first subspace
$[6,16]\times[5,7]\times[-3,2]$ generated by branching is quickly discarded,
while the second $[6,16]\times[3,5]\times[-3,2]$ is subsequently reduced to
the optimal solution.

###  3.2 Scenario 2: feasible set with several catalog items

In this scenario, we consider the catalog of the first scenario with an
additional item described by properties $y=(1,-1)$ (Table 3). The catalog
items in their initial domain $[0,20]\times[-10,10]$ are represented in Figure
3(a).

Table 3: Scenario 2: catalog $C$ Item | Property $y_{1}$ | Property $y_{2}$
---|---|---
1 | 4 | -8
2 | 3 | 2
3 | 7 | -3
4 | 14 | 8
5 | 19 | -8
6 | 1 | -1

Table 4 describes the successive filtering and upper bounding phases (second
column) that reduce the ranges of $x$ and $y$ starting from
$(x,y_{1},y_{2})\in[0,16]\times[0,20]\times[-10,10]$: the third column
contains for each phase the current ranges of $x$ and $y$, while the fourth
column presents the resulting ranges. The first column indicates the
corresponding subfigures in Figure 3. The shortcut $3^{-}$ stands for
$\sqrt[3]{27-\varepsilon}$, a quantity close to but strictly smaller than 3.

(a) Initial range of properties $y_{1}$ and $y_{2}$

(b) HC4Revise with respect to $c_{1}$

(c) Clutch with respect to $c_{2}$

(d) HC4Revise with respect to $f\leq\tilde{f}-\varepsilon$

(e) Clutch with respect to $c_{2}$

(f) HC4Revise with respect to $c_{1}$

Figure 3: Scenario 2: alternance of filtering and upper bounding phases Table 4: Scenario 2: alternance of filtering and upper bounding phases | Phase | Current box | Result
---|---|---|---
(b) |  HC4Revise($c_{1}$) | $[0,16]\times[0,20]\times[-10,10]$ | $[0,16]\times[0,8]\times[-4,4]$
(c) |  Clutch($c_{2}$) | $[0,16]\times[0,8]\times[-4,4]$ | $[0,16]\times[1,7]\times[-3,2]$
|  HC4Revise($c_{1}$) | $[0,16]\times[1,7]\times[-3,2]$ | $[2,16]\times[1,7]\times[-3,2]$
(d) |  HC4Revise($f\leq\tilde{f}-\varepsilon$) | $[2,16]\times[1,7]\times[-3,2]$ | $[2,16]\times[1,3^{-}]\times[-3,2]$
(e) |  Clutch($c_{2}$) | $[2,16]\times[1,3^{-}]\times[-3,2]$ | $[2,16]\times[1,1]\times[-1,-1]$
(f) |  HC4Revise($c_{1}$) | $[2,16]\times[1,1]\times[-1,-1]$ | $\mathbf{[3,3]\times[1,1]\times[-1,-1]}$

The same color code as in Figure 2 is used in Figure 3. The first contraction
step discards three of the catalog items. Constraint propagation between
HC4Revise and Clutch reduces the domain to the convex hull of catalog items 2,
3 and 6. An upper bounding phase (not shown in the figure) then selects the
feasible item 2 (with properties $(y_{1},y_{2})=(3,2)$) and the corresponding
value $x=10$ and updates the best known upper bound: $\tilde{f}\leftarrow 27$.
The subsequent contraction phase with respect to $f\leq\tilde{f}-\varepsilon$
eliminates a large part of the domain (including item 2), which is then
reduced to the optimal solution (item 6) after a Clutch phase. In this case,
no branching phase is required.

##  4 Conclusion

We introduced a generic framework for mixed categorical global optimization
with catalogs that produces the globally optimal solution. Based on spatial
branch and bound, it optimizes the underlying properties of the catalog items
and does not require any modeling effort (such as reformulation with binary
variables) from the user. The ranges of the categorical properties, similarly
to those of the continuous variables, are reduced or discarded during
constraint programming and branching. A novel catalog-based contractor called
Clutch, invoked upon filtering and branching, guarantees consistency between
the categorical properties that are being optimized and the existing catalog
items.

We presented an example with two different catalogs that demonstrates how an
alternance of constraint propagation between Clutch and classical contractors,
upper bounding and branching iteratively reduces the ranges of the continuous
and property variables until convergence to the global minimum.

The genericity and flexibility of our approach make it an excellent strategy
for a wide range of applications in which categorical variables arise. Since
it relies on the same mechanisms as spatial branch and bound codes, it can be
easily implemented in off-the-shelf solvers with minimal effort.

## References

  * Abhishek et al., [2010] Abhishek, K., Leyffer, S., and Linderoth, J. T. (2010).  Modeling without categorical variables: a mixed-integer nonlinear program for the optimization of thermal insulation systems.  Optimization and Engineering, 11(2):185–212.
  * Agurok, [2019] Agurok, I. (2019).  Multi-extremum optimization in lens design: navigation through merit function valleys maze.  arXiv preprint arXiv:1907.08676.
  * Alliot et al., [2012] Alliot, J.-M., Durand, N., Gianazza, D., and Gotteland, J.-B. (2012).  Finding and proving the optimum: Cooperative stochastic and deterministic search.  In ECAI, pages 55–60.
  * Araya et al., [2012] Araya, I., Trombettoni, G., and Neveu, B. (2012).  A contractor based on convex interval Taylor.  In International Conference on Integration of Artificial Intelligence (AI) and Operations Research (OR) Techniques in Constraint Programming, pages 1–16. Springer.
  * Barjhoux et al., [2017] Barjhoux, P.-J., Diouane, Y., Grihon, S., Bettebghor, D., and Morlier, J. (2017).  Mixed variable structural optimization: toward an efficient hybrid algorithm.  In World Congress of Structural and Multidisciplinary Optimisation, pages 1880–1896. Springer.
  * Belotti et al., [2010] Belotti, P., Cafieri, S., Lee, J., and Liberti, L. (2010).  Feasibility-based bounds tightening via fixed points.  In International conference on combinatorial optimization and applications, pages 65–76. Springer.
  * Benhamou et al., [1999] Benhamou, F., Goualard, F., Granvilliers, L., and Puget, J.-F. (1999).  Revising hull and box consistency.  In International Conference on Logic Programming, pages 230–244. MIT press.
  * Blum et al., [2011] Blum, C., Puchinger, J., Raidl, G. R., and Roli, A. (2011).  Hybrid metaheuristics in combinatorial optimization: A survey.  Applied soft computing, 11(6):4135–4151.
  * Brown and Hwang, [1993] Brown, D. R. and Hwang, K.-Y. (1993).  Solving fixed configuration problems with genetic search.  Research in Engineering Design, 5(2):80–87.
  * Carlson, [1996] Carlson, S. E. (1996).  Genetic algorithm attributes for component selection.  Research in Engineering Design, 8(1):33–51.
  * Carlson-Skalak et al., [1998] Carlson-Skalak, S., White, M. D., and Teng, Y. (1998).  Using an evolutionary algorithm for catalog design.  Research in Engineering Design, 10(2):63–83.
  * Chabert and Jaulin, [2009] Chabert, G. and Jaulin, L. (2009).  Contractor programming.  Artificial Intelligence, 173(11):1079–1100.
  * Cleary, [1987] Cleary, J. G. (1987).  Logical arithmetic.  Future computing systems, 2(2):125–149.
  * Cotta et al., [1995] Cotta, C., Aldana, J., Nebro, A. J., and Troya, J. M. (1995).  Hybridizing genetic algorithms with branch and bound techniques for the resolution of the TSP.  In Artificial neural nets and genetic algorithms, pages 277–280. Springer.
  * Gallardo et al., [2007] Gallardo, J. E., Cotta, C., and Fernández, A. J. (2007).  On the hybridization of memetic algorithms with branch-and-bound techniques.  IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), 37(1):77–83.
  * Gao et al., [2018] Gao, H., Breitkopf, P., Coelho, R. F., and Xiao, M. (2018).  Categorical structural optimization using discrete manifold learning approach and custom-built evolutionary operators.  Structural and Multidisciplinary Optimization, 58(1):215–228.
  * Ichida and Fujii, [1979] Ichida, K. and Fujii, Y. (1979).  An interval arithmetic method for global optimization.  Computing, 23(1):85–97.
  * Kearfott and Novoa III, [1990] Kearfott, R. B. and Novoa III, M. (1990).  Algorithm 681: INTBIS, a portable interval Newton/bisection package.  ACM Transactions on Mathematical Software (TOMS), 16(2):152–157.
  * Lindroth and Patriksson, [2011] Lindroth, P. and Patriksson, M. (2011).  Pure categorical optimization: a global descent approach.  Department of Mathematical Sciences, Division of Mathematics, Chalmers ….
  * Messine, [1997] Messine, F. (1997).  Méthodes d’optimisation globale basées sur l’analyse d’intervalle pour la résolution de problèmes avec contraintes.  PhD thesis, Toulouse, INPT.
  * Nedělková et al., [2019] Nedělková, Z., Cromvik, C., Lindroth, P., Patriksson, M., and Strömberg, A.-B. (2019).  A splitting algorithm for simulation-based optimization problems with categorical variables.  Engineering Optimization, 51(5):815–831.
  * Neumaier and Shcherbina, [2004] Neumaier, A. and Shcherbina, O. (2004).  Safe bounds in linear and mixed-integer linear programming.  Mathematical Programming, 99(2):283–296.
  * Ryoo and Sahinidis, [1995] Ryoo, H. S. and Sahinidis, N. V. (1995).  Global optimization of nonconvex NLPs and MINLPs with applications in process design.  Computers & Chemical Engineering, 19(5):551–566.
  * Sherali and Adams, [2013] Sherali, H. D. and Adams, W. P. (2013).  A reformulation-linearization technique for solving discrete and continuous nonconvex problems, volume 31.  Springer Science & Business Media.
  * Sunaga, [1958] Sunaga, T. (1958).  Theory of interval algebra and its application to numerical analysis.  RAAG memoirs, 2(29-46):209.
  * Vanaret, [2015] Vanaret, C. (2015).  Hybridization of interval methods and evolutionary algorithms for solving difficult optimization problems.  PhD thesis, INP Toulouse.
  * Zhang and Liu, [2007] Zhang, X. and Liu, S. (2007).  A new interval-genetic algorithm.  In Third International Conference on Natural Computation (ICNC 2007), volume 4, pages 193–197. IEEE.