# Phononic Bandgap Programming in Kirigami

###### Abstract

This study investigates the programming of elastic wave propagation bandgaps
in a kirigami material by intentionally sequencing its constitutive mechanical
bits. Such sequencing exploits the multi-stable nature of the stretched
kirigami, allowing each mechanical bit to switch between two stable equilibria
with different external shapes (aka. “(0)” and “(1)” states). Therefore, by
designing the sequence of (0) and (1) bits, one can fundamentally change the
underlying periodicity and thus program the phononic bandgap frequencies. To
this end, this study develops an algorithm to identify the unique
periodicities generated by assembling “$n$-bit strings” consisting of $n$
mechanical bits. Based on a simplified geometry of these $n$-bit strings, this
study also formulates a theory to uncover the rich mapping between input
sequencing and output bandgaps. The theoretical prediction and experiment
results confirm that the (0) and (1) bit sequencing is a versatile approach to
programming the phonic bandgap frequencies. Moreover, one can additionally
fine-tune the bandgaps by adjusting the global stretch. Overall, the results
of this study elucidate new strategies for programming the dynamic responses
of architected material systems.

###### keywords:

Phononic Bandgap, Kirigami, Multi-Stability, Property Programming

Hesameddin Khosravi ∗, Suyi Li +

\+ Department of Mechanical Engineering, Virginia Tech
153 Durham Hall, 1145 Perry Street, Blacksburg, VA USA, 24060
* Department of Mechanical Engineering, Clemson University
224 Fluor Daniel Building, 216 South Palmetto Blvd., Clemson, SC USA, 29631

Correspondence: Prof. Suyi Li<EMAIL_ADDRESS>
##  1 Introduction

Advanced materials systems with on-demand property programming have enabled
new capabilities for many applications. Like in computer programming, where
one can design an input instruction to execute computation and obtain results,
a programmable material can provide a unique output response according to an
input set [1, 2]. It is worth emphasizing that these inputs should involve
multiple and typically local components. So programmable materials differ
fundamentally from the more traditional responsive materials that only react
to a global and ambient stimulus (such as heat-responsive shape memory alloys
[3, 4] or light-responsive liquid crystal polymers [5, 6, 7]). Moreover, users
should be able to reset and re-design the input set for different output
properties.

Because of the local and complex nature of the input set, almost all property
programming has been achieved in mechanical metamaterials consisting of an
assembly of unit cells with carefully designed geometries. In particular, if
each unit cell possesses two stable equilibria (aka. bi-stability), one can
discretize the input space for versatile programming. Intuitively, we can
refer to these bistable unit cells as “mechanical bits” and designate one
stable configuration as the state (0) and the other as the state (1) – similar
to the binary number system in digital computation. This way, we can program
the material property by inputting different combinations of (0) and (1) bits.
For example, in origami-based cellular structures with $n$ bi-stable
mechanical bits [8, 9], one can obtain the lowest elastic modulus if all bits
are in the state (0), the highest modulus when all bits are in the state (1),
and other $n-1$ levels of intermediate moduli if only portions of the
mechanical bits are in the state (1).

However, almost all material programming studies so far are limited to static
responses, so the permutation (or sequence) of inputs does not play a role. In
the examples above, the origami cellular material shows the same elastic
modulus as long as the total number of (1) bits does not change, and the
positions of these (1) bits have no influence. Therefore, we aim to explore
dynamic property programming by controlling the input sets’ _sequences_.
Specifically, we show how to program the phononic wave propagation bandgaps by
sequencing the mechanical bits in a kirigami-based metamaterial.

Kirigami material exhibiting phononic wave bandgap is an ideal platform to
investigate dynamic material property programming via input sequencing. The
phononic bandgaps – a unique phenomenon by which mechanical waves cannot
propagate through periodic media due to Braggs scattering [10] or local
resonance [11] – is one of the most fundamental dynamic properties in a
metamaterial system. The ability to manipulate wave propagation via bandgap
has enabled countless applications in noise mitigation [12], wave transmission
control [13], acoustic cloaking [14], non-reciprocity [15], and even wave-
based mechanical computation [16]. On the other hand, kirigami – the ancient
art of paper cutting – has become a versatile engineering framework for
designing and constructing multi-functional material systems [17]. Careful
cutting can impart stretchability to thin sheet materials like metal alloys,
graphene [18], and electronics [19], thus fostering new designs of flexible
sensors [20], robots [21, 22], and metamaterials [23, 24, 25]. In particular,
a stretched kirigami sheet can buckle out-of-plane and exhibit multi-stability
[26]. One can also use simple stretching to create and control periodicity to
generate tunable wave propagation bandgaps [27]. Therefore, the results of
this study can have broad applications to many engineering systems of vastly
different sizes and functionalities.

This study uses a simple kirigami sheet containing uniformly distributed
parallel cuts (Figure 1). When this kirigami is subject to significant in-
plane stretch, its slender facets could buckle and rotate out-of-plane,
generating a three-dimensional periodic topology. More importantly, each
elementary unit in the stretch-buckled kirigami exhibits bi-stability, so they
create a linear array of mechanical bits. This setup allows us to program
input – by sequencing the (0) and (1) bits – and fundamentally alter the
underlying periodic nature. Eventually, such periodic sequencing of mechanical
bits can change the elastic wave’s propagation direction from longitudinal to
out-of-plane, create a mechanical impedance mismatch at the junctions between
adjacent bits, and generate a wave bandgap at predictable frequencies.

In what follows, we elucidate the working principles of (0) and (1) bit input
sequencing, uncover the corresponding output bandgap programming, and validate
the results via extensive experiments.

Figure 1: Design, multi-stability, and input sequencing of the stretch-buckled
kirigami sheet. (a) The design parameters for a kirigami with parallel cuts.
Here, the definition of facets and kirigami bits are highlighted. (b) A
stretched kirigami sample containing five mechanical bits at either (0) or (1)
states. The $+$ and $-$ facet rotations in each bit are highlighted for
clarity. Note that the rotation direction must be the same between the two
adjacent facets at the kirigami bit’s boundary. Transverse elastic waves,
denoted by $U(x,t)$, are sent to propagate through the kirigami sheet. (c) The
unique input sequences generated by the 1, 2, and 3-bit strings, shown by a
paper-based kirigami sample. The inserted gray box summarizes the unique
periodicities introduced by 4 and 5-bit strings. Here, the vector
$\mathbf{e}_{1}$ defines the fundamental periodic unit underpinning input
sequence, and it is used for calculating the bandgap based on plane wave
expansion method.

##  2 Sequencing the Input

This study focuses on the transverse elastic wave propagation in a stretched
kirigami sheet with a “zig-zag” uniform distribution of parallel slit cuts
shown in Figure 1(a). A few simple design parameters — including the cut size
$l_{a}$, $l_{b}$, $l_{c}$, and the spacing between cuts $w$ — can fully define
the overall cutting pattern. When the kirigami sheet is stretched in-plane
along the direction perpendicular to the cuts (the $x-$direction), stress
starts to concentrate near the cut tips, eventually forcing the facets to
deform out-of-plane as the stretching reaches a critical level [28]. If the
kirigami sheet is stretched further after buckling, its long and slender
facets will rotate accordingly.

Importantly, these facet rotations can occur clockwise or counter-clockwise
due to the symmetric nature of kirigami (Figure 1b, for clarity, we denote the
counter-clockwise facets by “$+$,” and clockwise ones by “$-$” ). In practical
experimentations, the direction of facet rotation occurs randomly, dictated by
fabrication imperfections and loading conditions. However, if the global
stretch is constant, one could locally switch some facets’ rotation
directions, and the kirigami sheet can settle into the new configuration [29].
Therefore, the stretched kirigami sheet is multi-stable, and each stable
configuration corresponds to a unique sequence of $+$ and $-$ rotated facets.

To better describe the programmable periodicity created by such multi-
stability, we divide the stretch kirigami sheet into a linear array of
bistable mechanical bits (or “kirigami bit” hereafter), each consisting of two
slender facets (Figure 1a). A kirigami bit is at state (0) if its two facets
have the same rotation direction (aka. $++$ or $--$), or state (1) if the
rotation directions are opposite ($+-$ or $-+$). Although $++$ and $--$ facet
pairs are mirror-symmetric to each other with respect to the $x-y$ reference
plane, they have the same geometry for an incoming elastic wave. The same
principle applies to the $+-$ and $-+$ pairs.

Our experiment efforts also reveal a critical constraint regarding the
admissable facet rotation distributions: The adjacent facets’ rotation
direction must be the same at the boundary between two kirigami bits. For
example, to construct a pair of two (0) bits (labeled as “(00)”), one can set
the underpinning facet rotations to $++\;++$ (or $--\;--$). But the $++\;--$
(or $--\;++$) combination turn out to be unstable in our experiment. Likewise,
for an (11) pair, its underlying facet rotation must be $-+\;+-$ (or
$+-\;-+$). The same principle applies to (01) and (10) pairs (Figure 1b). If
we violate this constraint and force the two adjacent facets at the kirigami
bit’s boundary to rotate oppositely, the stretched kirigami sheet’s neutral
axis would shift locally away from the $x-y$ reference plane where the
external loads are located, creating an imbalance in the resultant internal
force and bending moment.

Then, we combine $n$ kirigami-bits into a string – which we refer to as
“$n$-bit string” hereafter – and infinitely repeat such a string into a
periodic topology to create the “input sequence.” This way, the arrangement of
(0) and (1) bits in the $n$-bit string offers a mechanism for bandgap
programming. The complexity of the periodicity in the input sequence is
directly related to the $n$-bit string’s length. For example, if the string
contains only one kirigami bit ($n=1$), there are only two possible
configurations: (0) or (1). So repeating a 1-bit string provides two periodic
input sequences $(\cdots 0\cdots)$ or $(\cdots 1\cdots)$ (Figure 1c). A 2-bit
string can show four possible configurations: (00), (01), (10), and (11).
However, repeating (00) or (11) gives the same periodicity as those from the
1-bit string. Moreover, infinitely repeating (01) and (10) creates the same
periodicity. Therefore, the 2-bit string offers one additional periodicity
than the 1-bit string, so there are three unique periodic input sequences in
total: $(\cdots 0\cdots)$, $(\cdots 1\cdots)$, and $(\cdots 01\cdots)$.
Similarly, a 3-bit string offers two new unique periodicities compared to the
2-bit string: $(\cdots 001\cdots)$ and $(\cdots 011\cdots)$, bringing the
total number of unique input sequences to 5.

As $n$ increases, the number of unique periodicities increases rapidly. A
generic $n$-bit string in the stretched kirigami contains a series of $n$
mechanical bits, and each can settle either in the stable state (0) or (1).
Therefore, the $n$-bit string is similar to an $n$-bit binary number and has
$2^{n}$ unique configurations (just like the $n$-bit binary number has $2^{n}$
unique values). However, as we just see in the 2 and 3-bit string examples,
when we repeat and assemble $n$-bit strings to create an infinite periodic
array as the input sequence, many show redundant periodicity. For example, a
4-bit string has 16($=2^{4}$) unique configurations, but it does not create 16
unique periodicities. There are two different reasons for the redundancy.

The first reason is that some configurations in the $n$-bit kirigami string
have the same periodicity as those generated by the smaller bit strings. For
example, one can configure the 4-bit string into $(0000)$. If we repeat and
assemble this configuration, we will end up with an infinite array of 0 bits:
$(\cdots 000000000000\cdots)$. However, its periodicity is already considered
by the 1-bit kirigami string, as shown in Figure 1(c). Similarly, $(0101)$ and
$(1010)$ configurations generate the same periodicity from the 2-bit string.

The second reason for the redundant periodicity is that some configurations in
the $n$-bit string create identical structures once they are repeated and
assembled into an infinite input sequence. Again in the 4-bit string example,
$(1110)$ and $(1011)$ are two different configurations. Repeating these two
configurations creates $(\cdots 111011101110\cdots)$ and $(\cdots
101110111011\cdots)$, respectively. However, by careful inspection, one can
easily see that these two infinite arrays have the same underlying periodicity
(Figure 2a-b).

Here, we borrow the masking and convolution concepts from image processing to
search for this second type of redundant periodicity. Suppose we select a
unique configuration in the $n$-bit kirigami string as the “basis,” we can
generate $n-1$ “similarities” by shifting an $n$-bit-long mask successively
along the corresponding infinite sequence (Figure 2c). Then, we can select
another configuration as the new basis and generate its $n-1$ similarities. If
identical items exist among these similarities, these two configurations are
considered to create the same periodicity (Figure 2c). Moreover, since each
similarity itself is also a unique configuration of an $n$-bit string, we can
deduce that for each configuration of the $n$-bit kirigami string, there exist
$n-1$ different configurations that could generate the identical periodic
sequence.

Figure 2: Identifying the redundant periodicity in the $n$-bit kirigami
string, using $n=4$ as an example. (a, b) Here, we choose three unique
configurations in the 4-bit string as the bases and repeat these bases to
create the periodic input sequences. (c) Using a 4-bit-long mask, we make
3(=4-1) similarities for each basis by shifting the mask. For example, the
similarities of the (0011) basis are (0110), (1100), and (1001). Then, we can
identify redundant periodicities by searching for identical items among the
similarities. In this example, the two bases (1110) and (1101) generate the
same periodicity due to the many identical similarities, highlighted by the
color-coded boxes.

Therefore, denote $P_{n}$ as the number of _new_ periodicities introduced by
the $n$-bit kirigami string. Based on the discussions above, $P_{1}=2$,
$P_{2}=1$, and $P_{3}=2$. For a $n$-bit string, we have:

$P_{n}=2^{n}-\sum_{i}{iP_{i}}-P_{n}(n-1),$ (1)

where $2^{n}$ is the number of all possible configurations in the $n$-bit
string. $i$ are the smaller divisors of $n$ (i.e., $i\leq n/2$), so
$\sum_{i}{iP_{i}}$ corresponds to the first type of periodicity redundancy
from the smaller bit strings. Here, the coefficient $i$ before ${P_{i}}$
includes all of the similarities in the smaller bit strings. On the other
hand, $P_{n}(n-1)$ corresponds to the second type of redundancy mentioned
above. Simplifying Equation (1) yields

$P_{n}=\frac{2^{n}-\sum_{i}iP_{i}}{n},\quad\quad i=\text{divisor of }n\text{,
and }i\leq\frac{n}{2}.$ (2)

##  3 Mapping the Input Sequence to Output Bandgap

Once the input sequence is programmed, the stretched kirigami takes a long-
slender shape with a finite thickness in the out-of-plane $z-$direction.
Therefore, we assume negligible shear deformation and rotational inertia in
kirigami’s cross-sections so that the classical Euler-Burnulli bean equation
governs the elastic wave propagation:

$\frac{\partial^{2}}{\partial
x^{2}}\left[EI(x)\frac{\partial^{2}{U(x,t)}}{\partial x^{2}}\right]+\rho
A(x)\frac{\partial^{2}{U(x,t)}}{\partial t^{2}}=0,$ (3)

where $\rho$ and $E$ are the mass density and Young’s modulus of the
constitutive sheet material, respectively. $U(x,t)$ is the out-of-plane
transverse wave displacement field in the $z-$direction. $A(x)$ and $I(x)$ are
the spatial distribution of cross-section area and area moment of inertia,
respectively. Note that the in-plane stretching force does not play any
significant role here because right after the stretch-induced buckling occurs,
the kirigami’s reaction force does not increase significantly unless the
stretch becomes high enough to generate plastic deformation [26].

Calculating $A(x)$ and $I(x)$ is quite challenging because the kirigami facets
have complex-curved shapes orientated at an angle with respect to the neutral
$x-y$ plane. Therefore, we adopt the simplification method proposed in the
previous literature that uses virtual origami folds and flat surfaces to
approximate the stretch-buckled kirigami’s geometry [27] (Figure 3). This
method has been proven successful in providing reasonably accurate predictions
of bandgap frequency without incurring a high computational cost.

Figure 3: Geometric simplification of stretched kirigami (a) A close-up view
of a stretched kirigami, showing its complex curved facets. (b) The
corresponding simplified geometry using the virtual folds and flat surface
assumptions. Notice the definition of $\varphi$ angle. (c) A detailed view of
the simplified (0) and (1) bits. Here, the straight cross-sections at
different positions in the kirigami bit are highlighted. (d) Examples of
$A(x)$ and $I(x)$ distributions ($l_{a}=0.044$m, $l_{b}=0.022$m,
$l_{c}=0.086$m, $w=0.016$m, $\varphi_{\text{(0)}}=0.66$, and
$\varphi_{\text{(1)}}=1.23$).

According to the geometric simplification, the (0) bit consists of 8 flat
surfaces and 8 fold lines, and all these fold lines are perpendicular to the
slit cuts. On the other hand, a (1) bit has 2 additional fold lines aligned
with the cuts (Figure 3). In this way, the post-buckling deformation of the
kirigami becomes a rigid-folding motion with one kinematic degree of freedom.
To this end, we choose the dihedral angle ($\varphi$) between the rotated
facets in the kirigami bit and the $y-z$ reference plane as the independent
variable (Figure 3b). When the kirigami sheet is un-deformed (or flat in the
$x-y$ reference plane), $\varphi$ takes the maximum value
($\varphi_{max}=\pi/2$); when the buckled kirigami sheet is fully stretched,
$\varphi$ take the minimum value:

$\varphi_{min}=\tan^{-1}\left(\frac{2w}{l_{c}}\right).$ (4)

The length of the kirigami mechanical bits is:

$a=\frac{2w}{\sin\varphi}.$ (5)

Using the geometric simplification, one can represent the in-plane stretching
stiffness of the post-bucked kirigami bits by assuming that the virtual folds
behave like hinges with embedded torsional springs and that the surface
between them is rigid [30]. Intuitively, (1) bit is stiffer in stretching than
(0) bit because it has two additional virtual folds (thus two additional
torsional springs). Therefore, (1) bit is shorter than a (0) bit along the in-
plane $x-$direction because they are mechanically connected in series. We
stretched a kirigami sample containing four (0) bits and four (1) bits to
different global lengths and measured the kirigami bits’ lengths locally, and
found a simple linear relationship:

$a_{(1)}=\frac{a_{(0)}+33\text{mm}}{2.53}\pm 0.43\text{mm}.$ (6)

Moreover, the cross-sections of simplified kirigami bits are straight lines
(Figure 3c), so the calculation of $A(x)$ and $I(x)$ becomes manageable.
Figure 3(d) summarizes the results corresponding to the (0) and (1) bits.
Section 1 of the supplement materials details the underpinning formulations.

To analyze the wave propagation behavior in the stretch-buckled kirigami
sheet, we employ the plane wave expansion (PWE). First, we define
$\mathbf{R}=a_{1}\mathbf{e}_{1}$ as a translational vector in the _direct
space_ to represent the periodicity of the programmed kirigami sheet. In other
words, we can construct the kirigami sheet by a set of infinite translational
operations on the $n$-bit string based on $\mathbf{R}$. In this case, $a_{1}$
is an integer number, and the lattice vector $\mathbf{e}_{1}$ connects the two
ends of the $n$-bit string (Figure 1c). We then define the reciprocal lattice
vector $\mathbf{G}$ in that $e^{-i\;\mathbf{G}\cdot\mathbf{R}}=1$. Here,
$m_{1}$ is another integer, $\mathbf{G}=m_{1}\mathbf{b}_{1}$, and
$\mathbf{b}_{1}$ is the lattice vector in the _reciprocal space_. Equivalent
to the lattice vector $\mathbf{e}_{1}$ in the direct space, infinite
translation operations of the reciprocal lattice vectors $\mathbf{b}_{1}$
defines the _reciprocal space_ , which has similar symmetry as the direct
space. Like in the direct lattice, an equivalent periodic unit can then be
defined in the reciprocal space based on $\mathbf{b}_{1}$, giving us the First
Irreducible Brillouin Zone (FBZ).

Then, a separation of variables gives $U(x,t)=\tilde{U}(x)e^{-i\omega t},$
where $\omega$ is the harmonic oscillation frequency. Based on the Blotch
theorem, we can formulate the spatial term $\tilde{U}(x)$ for the whole
kirigami sheet into a product of plane-wave and periodic functions in the
First Irreducible Brillouin Zone so that

$\tilde{U}(x)=e^{i(k_{1}b_{1}x)}\sum_{n_{1}}{\hat{U}({n_{1}b_{1}})}e^{i(n_{1}b_{1}x)},$
(7)

where $b_{1}$ is the magnitude of reciprocal lattice vector
($b_{1}=|\mathbf{b}_{1}|$), $n_{1}$ is an integer, and ${k}_{1}$ is a wave
number within the FBZ in that ${k}_{1}\in[0,\frac{1}{2}]$. We further expand
the distribution of cross-section area $A(x)$ and area moment of inertia
$I(x)$ of an $n$-bit string into Fourier series so that

$\displaystyle I(x)$ $\displaystyle=\sum_{m_{1}}\hat{I}(m_{1}b_{1})e^{i(2\pi
m_{1}x)},$ (8) $\displaystyle A(x)$
$\displaystyle=\sum_{m_{1}}\hat{A}(m_{1}b_{1})e^{i(2\pi m_{1}x)}.$ (9)

By substituting the formulations above into the governing equations of motion
(3), we obtain the eigenvalue problem for calculating the dispersion curves
and bandgap frequencies.

$\sum_{n_{1}}E(n_{1}+k_{1})^{2}(k_{1}+m_{1}+n_{1})^{2}\hat{I}(m_{1}b_{1})-\omega^{2}\rho\hat{A}(m_{1}b_{1})\;=0.$
(10)

Section 2 of the supplement material summarizes the dispersion curves from 1,
2, and 3-bit strings, and interested readers can refer to the authors’
previous publication for a more detailed formulation underpinning the PWE
method [27].

##  4 Programming and Fine-Tuning the Output Bandgap

To validate the mapping between input programming and output bandgap, we
compare the theoretical prediction to experiment measurements based on a
Nylon-based kirigami sample with 13 mechanical bits ($E=5.09$GPa,
$\rho=1200$Kg/m3, material thickness $t=0.00175$m, with cut sizes
$l_{a}=0.044$m, $l_{b}=0.022$m, $l_{c}=0.086$m, and cut spacing $w=0.016$m).
Figure 8 in Section 3 of the supplement material details the experiment setup.
We attach one end of the kirigami sheet to a shaker (Labworks DB-140 with
Pa-141 amplifier) to provide the out-of-plane ($z-$direction) excitation. The
other end is attached to a rigid end fixture via thin tape to mimic the free
boundary conditions. To measure the wave transmissibility through the stretch-
buckled kirigami sheet, we use a signal generator (Tektronix AFG3022c) to
generate harmonic excitations with sweeping frequencies and use two laser
vibrometers (Polytec OFC-5000) to measure the input and transmitted waves at
the two ends of kirigami. Then, we calculate the wave transmissibility,
defined as $TR=\tilde{U}_{\text{out}}/\tilde{U}_{\text{in}}$, where
$\tilde{U}_{\text{in}}$ and $\tilde{U}_{\text{out}}$ are the maximum
displacement of incoming and transmitted wave, respectively. This
transmissibility has been used widely to evaluate bandgap frequency [27]. It
is worth noting that as we change the sequence of input (0) and (1) bits, we
also adjust the global stretch of the kirigami so that the length of (0) and
(1) bits remain the same among different tests: $a_{(0)}=0.052$m,
$a_{(1)}=0.034$m.

Figure 4: Phononic bandgap programming in kirigami by input sequencing. (a)
Test results from 1, 2, and 3-bit strings validate the bandgap programming
conectp, with good agreement between theoretically predicted bandgaps and
experimentally measured wave transmissibility. (b) Summary of theoretically
predicted bandgaps of 1, 2, 3, and 4-strings, arranged according to the
percentage of their mechanical bits in the state (1). That is, the input
sequence at the top has the most (0) bits, and the one at the bottom has the
most (1) bits.

The analytical theory’s predictions and experimental results all confirm the
occurrence of Bragg’s bandgap in the stretch buckled kirigami sheet (Figure
4). More importantly, these bandgap frequencies can be effectively programmed
by simply sequencing the (0) and (1) bits. Based on the theoretical
prediction, if we use the 1-bit string to sequence the input, the first
bandgap occurs at 517 to 846Hz when the 1-bit string is (0) (aka. all kirigami
bits are in state 0). However, the bandgap jumps to a much higher frequency at
3576 to 3720Hz when the 1-string is (1). By using the 2-bit string (01) to
sequence the input, we can achieve a new bandgap at 2661 to 3163Hz. Similarly,
the 3-bit string offers two additional unique periodicities compared to the
2-bit string (aka. $P_{3}=2$), and we obtained two new bandgaps accordingly
(Figure 4a). Figure 4(b) summarizes the theoretically predicted bandgap
frequencies up to the 4-bit string. The bandgap generally moves to higher
frequencies and becomes narrower when more kirigami bits are in the state (1).
However, if the total number of (1) bits stays the same, the bandgap would
jump to different frequencies as we change the sequence of these (1) bits,
elucidating the rich programmability of the bandgaps in kirigami. For example,
the 2-bit (01) and 4-bit (0011) both have 50% of their kirigami bits in the
state (1), and they show different bandgap frequencies due to their unique
sequencing. Our experiment results based on 1, 2, and 3-bit strings agree well
with the theoretical prediction, validating the bandgap programmability.

While the experiment results and theoretical predictions show good agreement
regarding the bandgap frequencies for all input sequences, some minor
differences still occur. These discrepancies come from a combination of
different factors. For example, the simplified geometry used in our
theoretical prediction approximates the actual shapes of the stretch-buckled
kirigami bits, so it would inevitably introduce errors, especially for the (1)
bits. Moreover, the analytical model based on the plane wave expansion method
assumes an infinite number of mechanical bits in the kirigami sheet, while the
experiment prototype only has thirteen bits. The omission of the tension force
in the governing equation is another probable cause for discrepancies.
Regardless, these results firmly validate our assumption that by sequencing
the (0) and (1) mechanical bits, one can program the elastic wave propagation
bandgaps with considerable freedom. Moreover, our analytical model elucidates
the physical principles underpinning the mapping between input sequencing and
output bandgap frequency.

The input sequencing of (0) and (1) kirigami bits offer an avenue to program
the phononic bandgaps to different frequency ranges. In addition, we can fine-
tune the bandgap by slightly adjusting the global stretch. In other words, if
we consider the (0) and (1) bit sequencing as a “digital” and discrete method
to program bandgap frequency, adjusting the kirigami stretch is an “analog”
and continuous approach to fine-tune bandgaps. In a follow-up experiment, we
adjust the global stretch for the kirigami sheet – based on a 3-string (011)
input sequence – and continuously move the bandgap frequency (Figure 5).
Therefore, to achieve a targeted phononic bandgap frequency, one can first use
the (0) ad (1) bits input sequencing to program the bandgap near the target
and then adjust the global stretching to reach the target precisely.

Figure 5: Fine-tuning the bandgap by adjusting the global external stretch. By
lowering the global stretch from $a_{(1)}=0.034$m to $a_{(1)}=0.032$m, or
increasing the stretch so that $a_{(1)}=0.038$m, one can continuously adjust
the bandgap frequencies.

##  5 Summary

This study theoretically and experimentally examines the phonic bandgap
programming in stretched kirigami material with intentionally sequenced
mechanical inputs. The kirigami consists of a linear array of mechanical bits,
each exhibiting two distant stable equilibria due to the internal stress
generated by stretching. We refer to them as (0) and (1) states, inspired by
the binary number systems underpinning computer programming. Therefore, by
designing the sequence of these (0) and (1) bits, one can fundamentally change
the underlying periodicity in kirigami and thus shift the phononic bandgap
frequencies. We develop an algorithm to identify the unique periodicities
generated by repeating and assembling a string of $n$ kirigami bits. An
analytical theory is also formulated to uncover the rich mapping between input
sequencing and output bandgaps. This theory leverages a simplified geometry of
stretched kirigami using the virtual fold and flat surface assumption and
incorporates the plane wave expansion method to calculate the bandgap
structures. Our theoretical prediction and experiment results confirm that (0)
and (1) bit sequencing is a versatile and effective approach to programming
the phonic bandgap frequencies. Moreover, one can additionally fine-tune the
bandgaps by adjusting the global stretch. Overall, the results of this study
elucidate new approaches to programming the dynamic responses of architected
material systems.

Acknowledgements

The authors acknowledge the support from National Science Foundation
(CMMI-2240211 CAREER, 2240326).

## References

  * [1] Tie Jun Cui et al.  “Coding metamaterials, digital metamaterials and programmable metamaterials”  In _Light: science & applications_ 3.10 Nature Publishing Group, 2014, pp. e218–e218
  * [2] Shuo Liu and Tie Jun Cui  “Concepts, working principles, and applications of coding and programmable metamaterials”  In _Advanced optical materials_ 5.22 Wiley Online Library, 2017, pp. 1700624
  * [3] Junjiang Chen, Lauren E Hamilton, Patrick T Mather and James H Henderson  “Cell-Responsive Shape Memory Polymers”  In _ACS Biomaterials Science & Engineering_ ACS Publications, 2022
  * [4] Vagner Candido de Sousa, David Tan, Carlos De Marqui Jr and Alper Erturk  “Tunable metamaterial beam with shape memory alloy resonators: theory and experiment”  In _Applied Physics Letters_ 113.14 AIP Publishing LLC, 2018, pp. 143502
  * [5] Gülistan Koçer et al.  “Light-responsive hierarchically structured liquid crystal polymer networks for harnessing cell adhesion and migration”  In _Advanced materials_ 29.27 Wiley Online Library, 2017, pp. 1606407
  * [6] Lucia M Korpas, Rui Yin, Hiromi Yasuda and Jordan R Raney  “Temperature-responsive multistable metamaterials”  In _ACS Applied Materials & Interfaces_ 13.26 ACS Publications, 2021, pp. 31163–31170
  * [7] Jun Wu et al.  “Liquid crystal elastomer metamaterials with giant biaxial thermal shrinkage for enhancing skin regeneration”  In _Advanced Materials_ 33.45 Wiley Online Library, 2021, pp. 2106175
  * [8] Sattam Sengupta and Suyi Li  “Harnessing the anisotropic multistability of stacked-origami mechanical metamaterials for effective modulus programming”  In _Journal of Intelligent Material Systems and Structures_ 29.14 SAGE Publications Sage UK: London, England, 2018, pp. 2933–2945
  * [9] Larissa S Novelino et al.  “Untethered control of functional origami microrobots with distributed actuation”  In _Proceedings of the National Academy of Sciences_ 117.39 National Acad Sciences, 2020, pp. 24096–24101
  * [10] Andrea Bergamini et al.  “Phononic crystal with adaptive connectivity”  In _Advanced Materials_ 26.9 Wiley Online Library, 2014, pp. 1343–1347
  * [11] Gang Wang et al.  “Two-dimensional locally resonant phononic crystals with binary structures”  In _Physical review letters_ 93.15 APS, 2004, pp. 154302
  * [12] Roy H Olsson and I El-Kady  “Microfabricated phononic crystal devices and applications”  In _Measurement science and technology_ 20.1 IOP Publishing, 2008, pp. 012002
  * [13] Guangxin Liao et al.  “Acoustic metamaterials: A review of theories, structures, fabrication approaches, and applications”  In _Advanced Materials Technologies_ 6.5 Wiley Online Library, 2021, pp. 2000787
  * [14] Huanyang Chen and Che Ting Chan  “Acoustic cloaking and transformation acoustics”  In _Journal of Physics D: Applied Physics_ 43.11 IOP Publishing, 2010, pp. 113001
  * [15] Hussein Nassar et al.  “Nonreciprocity in acoustic and elastic materials”  In _Nature Reviews Materials_ 5.9 Nature Publishing Group, 2020, pp. 667–685
  * [16] Farzad Zangeneh-Nejad, Dimitrios L Sounas, Andrea Alù and Romain Fleury  “Analogue computing with metamaterials”  In _Nature Reviews Materials_ 6.3 Nature Publishing Group, 2021, pp. 207–225
  * [17] Jiayue Tao, Hesameddin Khosravi, Vishrut Deshpande and Suyi Li  “Engineering by Cuts: How Kirigami Principle Enables Unique Mechanical Properties and Functionalities”  In _Advanced Science_ Wiley-VCH, 2022, pp. to appear
  * [18] Melina K Blees et al.  “Graphene kirigami”  In _Nature_ 524.7564 Nature Publishing Group, 2015, pp. 204–207
  * [19] Braden M Li et al.  “Kirigami-inspired textile electronics: KITE”  In _Advanced Materials Technologies_ 4.11 Wiley Online Library, 2019, pp. 1900511
  * [20] Anne Katherine Brooks, Sudesna Chakravarty, Maryam Ali and Vamsi K Yadavalli  “Kirigami-Inspired Biodesign for Applications in Healthcare”  In _Advanced Materials_ 34.18 Wiley Online Library, 2022, pp. 2109550
  * [21] Ahmad Rafsanjani et al.  “Kirigami skins make a simple soft actuator crawl”  In _Science Robotics_ 3.15 American Association for the Advancement of Science, 2018, pp. eaar7555
  * [22] Hesameddin Khosravi, Steven M Iannucci and Suyi Li  “Pneumatic Soft Actuators With Kirigami Skins”  In _Frontiers in Robotics and AI_ 8 Frontiers Media SA, 2021
  * [23] Yigil Cho et al.  “Engineering the shape and structure of materials by fractal cut”  In _Proceedings of the National Academy of Sciences_ 111.49 National Acad Sciences, 2014, pp. 17390–17395
  * [24] Katia Bertoldi, Vincenzo Vitelli, Johan Christensen and Martin Van Hecke  “Flexible mechanical metamaterials”  In _Nature Reviews Materials_ 2.11 Nature Publishing Group, 2017, pp. 1–11
  * [25] Yichao Tang et al.  “Programmable kiri-kirigami metamaterials”  In _Advanced Materials_ 29.10 Wiley Online Library, 2017, pp. 1604262
  * [26] Ahmad Rafsanjani and Katia Bertoldi  “Buckling-induced kirigami”  In _Physical review letters_ 118.8 APS, 2017, pp. 084301
  * [27] Hesameddin Khosravi and Suyi Li  “Tunable Wave-Propagation Band gap via Stretching Kirigami Sheets”  In _Physical Review Applied_ 17.6 APS, 2022, pp. 064054
  * [28] Midori Isobe and Ko Okumura  “Initial rigid response and softening transition of highly stretchable kirigami sheet materials”  In _Scientific reports_ 6.1 Nature Publishing Group, 2016, pp. 1–6
  * [29] Yi Yang, Marcelo A Dias and Douglas P Holmes  “Multistable kirigami for tunable architected materials”  In _Physical Review Materials_ 2.11 APS, 2018, pp. 110601
  * [30] Suyi Li et al.  “Architected origami materials: how folding creates sophisticated mechanical properties”  In _Advanced materials_ 31.5 Wiley Online Library, 2019, pp. 1805282

## Supplement Material

### 1\. Simplified Geometry of Kirigami Bits

To analytically investigate the band gap distribution and its correlation to
input sequencing, we simplify the kirigami sheet geometry by assuming that
each mechanical bit consists of flat surfaces connected by virtual fold lines
(Figure 3). In this way, the post-buckling deformation of the kirigami becomes
a rigid-folding motion with one kinematic degree of freedom.

Figure 6 details the simplified geometry of a kirigami bit. We choose the
dihedral angle ($\varphi$) between surface B and the $y-z$ reference plane as
the independent variable. When the kirigami sheet is un-deformed (or flat in
the $x-y$ reference plane), $\varphi$ takes the maximum value
($\varphi_{max}=\pi/2$); when the buckled kirigami sheet is fully stretched,
$\varphi$ take the minimum value:

$\varphi_{min}=\tan^{-1}\left(\frac{2w}{l_{c}}\right).$ (11)

We also denote the dihedral angle between surface C defined in Figure 6 and
the $x-z$ reference place as $\beta$, so that

$\beta=\cos^{-1}\left(\frac{2w}{l_{c}\tan\varphi}\right).$ (12)

Clearly, $\beta=\pi/2$ when the kirigami bit is un-deformed (flat), and
$\beta=0$ when the Kirigami sheet is fully stretched. The overall length of
the unit cell is

$a=\frac{2w}{\sin\varphi}.$ (13)

Another important geometric variable is the _projected_ length of surface C
along the $y$ axis as illustrated in Figure 6(b, e):

$d=w\frac{\tan\beta}{\tan\varphi}=\sqrt{\frac{l_{c}^{2}}{4}-\frac{w^{2}}{\tan^{2}\varphi}}=\frac{l_{c}}{2}\sin\beta.$
(14)

Based on the geometric variables defined above, we can derive the cross-
section area $A(x)$ and the corresponding area moment of inertia $I(x)$ of the
kirigami bit. Here, we discuss two different cases for (1) bit and (0) bit.

Case 1. (1) bit: The first case corresponding to the slightly stretched
Kirigami sheet after buckling in that $\pi/4<\varphi<\pi/2$ (assuming
$l_{c}>2W$). We can divide half of the mechanical bit into three sections
because the cross-section areas take distinct shapes in these three sections
(Figure 6c).

Figure 6: The simplified geometry of (1) and (0) kirigami bits. (a-c):
Isometric, front, and side views of the (1) bit, respectively.
(b${}_{\text{i}}$-b${}_{\text{iii}}$): Close-up front views of the cross-
section areas corresponding to different $x$ values. Notice the solid, dashed,
and dotted black lines are the cross-section area of the (1) bit in the three
different sections. (d-f) The different views of a (0) bit, where the closed-
up front views of its cross-section areas are detailed in
(e${}_{\text{i}}$-e${}_{\text{iii}}$).

The first section corresponds to $0<x<a/2-W\sin\varphi$, and the cross-
sectional areas involves surface B and surface C (solid lines in Figure 6bi).
The cross-section area is the summation of these two parts so that

$A_{i}(x)=\frac{t}{\sin\varphi}\left(l_{b}+2L_{i}^{C}\right)=\frac{t}{\sin\varphi}\left(l_{b}+\frac{x}{w\sin\varphi}\sqrt{4w^{2}+l_{c}^{2}\tan^{4}\varphi}\right),$
(15)

where $t$ is the thickness of the Kirigami sheet material. The bending moment
of inertia with respect to the $y$-axis also includes two parts:

$I_{i}(x)=I^{B}+l_{b}\frac{t}{\sin\varphi}\left(\frac{x}{\tan\varphi}\right)^{2}+2I_{i}^{C}+2\frac{x\tan\beta}{\cos\varphi}\frac{t}{\sin\varphi}\left(\frac{x\tan\varphi-x\cot\varphi}{2}\right)^{2},$
(16)

where the first two terms come from the cross-section in surface B (using
parallel axis theorem), and the second two terms come from the cross-section
in surface C. Here, $I^{B}$ is the bending moment of inertia of the surface B
cross-section with respect to its own center of mass in that

$I^{B}=\frac{l_{b}}{12}\left(\frac{t}{\sin\varphi}\right)^{3}.$ (17)

Similarly, $I_{i}^{C}$ is the bending moment of inertia with respect to its
own mass center in that

$I_{i}^{C}(x)=\frac{t}{12\sin\varphi}\frac{x\tan\beta}{\cos\varphi}\left[\frac{t^{2}}{\sin^{2}\varphi}+x^{2}(\tan\varphi+\cot\varphi)^{2}\right].$
(18)

The second section corresponds to $a/2-w\sin<x<w\sin\varphi$, and the Kirigami
cross-section area here involves surface A, B and C.

$A_{ii}(x)=\frac{t}{\sin\varphi}\left(l_{a}+l_{b}+2L_{ii}^{C}\right)=\frac{t}{\sin\varphi}\left(l_{a}+l_{b}+\sqrt{\frac{4w^{2}}{\tan^{4}\varphi}+l_{c}^{2}}\right).$
(19)

The bending moment of inertia with respect to the $y$-axis includes three
parts:

$\begin{split}I_{ii}(x)&=I^{A}+l_{a}\frac{t}{\sin\varphi}\left(\frac{a/2-x}{\tan\varphi}\right)^{2}+I^{B}+l_{b}\frac{t}{\sin\varphi}\left(\frac{x}{\tan\varphi}\right)^{2}\ldots\\\
&+2I_{ii}^{C}+2d\frac{t}{\sin\varphi}\left(\frac{a/4-x}{\tan\varphi}\right)^{2}.\end{split}$
(20)

Here, $I^{A}$ is the bending moment of inertia of the cross-section area in
surface A with respect to its own mass center so that

$I^{A}=\frac{l_{a}}{12}\left(\frac{t}{\sin\varphi}\right)^{3},$ (21)

and

$I_{ii}^{C}=\frac{t}{12\sin\varphi}\frac{W\tan\beta}{\tan\varphi}\left(\frac{t^{2}}{\sin^{2}\varphi}+\frac{a^{2}}{4\tan^{2}\varphi}\right).$
(22)

Finally, the third section of the half unit cell corresponds to
$\sin\varphi<x<a/2$. In this section, the cross-section area involve surface A
and C in that

$A_{iii}(x)=\frac{t}{\sin\varphi}\left(l_{a}+2L_{iii}^{C}\right)=\frac{t}{\sin\varphi}\left(l_{a}+\frac{a/2-x}{W\sin\varphi}\sqrt{4W^{2}+l_{c}^{2}\tan^{4}\varphi}\right)$
(23)

The corresponding bending moment of inertia with respect to the $y-$axis
includes two components in that,

$\begin{split}I_{iii}(x)&=I^{A}+l_{a}\frac{t}{\sin\varphi}\left(\frac{a/2-x}{\tan\varphi}\right)^{2}\ldots\\\
&+2I_{iii}^{C}+2\frac{(a/2-x)\tan\beta}{\cos\varphi}\frac{t}{\sin\varphi}\left[\frac{(a/2-x)(\tan\varphi-\cot\varphi)}{2}\right]^{2}.\end{split}$
(24)

where,

$I_{iii}^{C}(x)=\frac{t}{12\sin\varphi}\frac{(a/2-x)\tan\beta}{\cos\varphi}\left[\frac{t^{2}}{\sin^{2}\varphi}+(a/2-x)^{2}(\tan\varphi+\cot\varphi)^{2}\right].$
(25)

Case 2. (0) bit: This case corresponds to
$\varphi\in[\varphi_{min}\ldots\pi/4]$ (Figure 6d-f). Again, we can divide
half of the kirigami bit into three different sections based on the shape of
its cross-section. Here, the cross-sectional area $A(x)$ and bending moment of
inertia $I(x)$ in the first section $(0<x<w\sin\varphi)$ and third section
$(a/2-w\sin\varphi<x<a/2)$ are the same as those in the previous case. The
only difference is in section ii, where the cross-section only involves
surface C so that,

$A^{*}_{ii}(x)=2\frac{t}{\sin\varphi}L_{ii}^{C*}=\frac{t}{\sin\varphi}\sqrt{4w^{2}+l_{c}^{2}\tan^{4}\varphi},$
(26)

and

$I^{*}_{ii}(x)=2I_{ii}^{C*}+2\frac{t}{\sin\varphi}\left[(a/4-x)\tan\varphi\right]^{2},$
(27)

where

$I_{ii}^{C*}=\frac{t}{12\sin\varphi}w\tan\varphi\tan\beta\left(\frac{t^{2}}{\sin^{2}\varphi}+\frac{w^{2}}{\cos^{2}\varphi}\right)$
(28)

### 2\. Dispersion Relationship Based on 1, 2, and 3-Bit Strings

Figure 7: The dispersion relationship of kirigami sheets based on 1, 2, and
3-bit strings. The underpinning kirigami design and material properties are
summarized at the beginning of Section 4.

### 3\. Experiment Setup

Figure 8 details the complete experiment setup for measuring the transverse
elastic wave transmissibility over the stretch-buckled kirigami. It is worth
noting that the unit cells in Figure1 c appear different from each other, but
this is an optical illusion. All unit cells in the kirigami prototype, except
for the first and last one at the boundary, have the same shape and dihedral
angle. However, when the camera is placed in front of the kirigami prototype,
every unit cell is at a slightly different angle with respect to the camera
lens, so they appear different in the picture.

Moreover, we did not observe significant global plastic deformation in our
stretched kirigami sample because the zig-zag distributed parallel cuts are
intended to offer high stretchability and elastic response. Due to stress
concentration, there are small and local plastic deformations at the tip of
the slit cuts. These plastic deformations can be avoided by optimizing the cut
tip shape, and they should not significantly affect the elastic wave bandgaps
at low frequencies.

Figure 8: A schematic diagram of the test apparatus.


# A Robustness Analysis to Structured Channel Tampering
over Secure-by-design Consensus Networks

Marco Fabris and Daniel Zelazo M. Fabris is with the Department of Information
Engineering, University of Padua, Padua, Italy. D. Zelazo is with the Faculty
of Aerospace Engineering, Technion, Haifa, Israel. This research was made
possible thanks to the support of Ms. Beverly Bavil, Mr. David Dibner and
another anonymous donor. Corresponding author: D. Zelazo,
<EMAIL_ADDRESS>
###### Abstract

This work addresses multi-agent consensus networks where adverse attackers
affect the convergence performances of the protocol by manipulating the edge
weights. We generalize [1] and provide guarantees on the agents’ agreement in
the presence of attacks on multiple links in the network. A stability analysis
is conducted to show the robustness to channel tampering in the scenario where
part of the codeword, corresponding to the value of the edge weights, is
corrupted. Exploiting the built-in objective coding, we show how to compensate
the conservatism that may emerge because of multiple threats in exchange for
higher encryption capabilities. Numerical examples related to semi-autonomous
networks are provided.

###### Index Terms:

Agents-based systems, Network analysis and control, Secure consensus protocols

##  I Introduction

The consensus problem, consisting in the design of networked control
algorithms under which all individuals of a given multi-agent system (MAS)
attain an agreement on a certain quantity of interest [2], is commonly tackled
when it is required to achieve a global prefixed task. However, due to the
openness of communication protocols and the complexity of networks, the
agreement of MASs may be vulnerable to malicious cyber-attacks [3]. In
particular, if the agent sensors are threatened by an attacker, the measured
data may be unreliable or faulty. Indeed, the attack signals can even disrupt
the control performance of the group of agents through the communication
topology. Therefore, resilient solutions are required to ensure that MASs
fulfill consensus under security hazards [4]. Consequently, the secure control
of MASs is now a crucial issue to be investigated [5, 6, 7].

Several recent studies illustrate the importance of giving guarantees against
cyber-threats while these are attempting to disrupt a MAS that is trying to
reach consensus. In [8, 9], deception attackers injecting false data are
assumed to attack the agents or communication channels. To counteract this
kind of disruption, classic observers, impulsive control methods and event-
triggered adaptive cognitive control have been leveraged. Denial of service
(DoS) attacks then represent another challenging class of cyber-threats:
robust control techniques have been developed in [10, 11, 12] to ensure
sufficient levels of agreement over MASs under DoS providing guarantees based
on the maximum “quality of service” or Lyapunov theory.

As already introduced in [1], we embrace a different perspective. Instead of
developing tools to secure existing networks (see [13]), we provide inherently
secure embedded measures through the adoption of a network manager to
guarantee robust consensus convergence. For privacy and safety concerns, such
a manager is not allowed to access local states. Rather, it only intervenes in
an initial phase to ensure desired convergence performance via edge weight
assignment (see e.g. [14]) in a secure way. Nonetheless, differently from our
previous work, this paper is meant to generalize the secure-by-design
consensus (SBDC) dynamics towards multi-agent consensus networks where adverse
attackers affect the convergence performance through a structural hit to the
communication between manager and agents, thus corrupting edge weights in more
than a single link of the system. We summarize our main contributions below.

  * •

Two new guarantees based on the small gain theorem for the robust stability of
the agents’ agreement are given in both continuous and discrete time domains
when multiple network edges undergo a structured weight deviation from their
nominal values.

  * •

We introduce the notion of a resilience gap used to characterize the
conservatism of the robustness analysis. We show that for spanning trees, the
resilience gap is always zero even in multi-attack scenarios. We also discuss
how this gap can be reduced by modulating the encryption capabilities used in
the network.

The organization of the paper follows. Sec. II introduces the preliminary
notions for multi-agent consensus and the SBDC networks. Sec. III provides a
robustness analysis for the latter when subject to channel tampering modeled
as multiple edge weight perturbations. A numerical case study on semi-
autonomous networks is reported in Sec. IV to assess such theoretical results
and, lastly, concluding remarks and future works are sketched in Sec. V.

#### Notation

The set of real numbers, the $l$-dimensional (column) vector whose elements
are all ones and the $l$-dimensional (column) null vector are denoted by
$\mathbb{R}$, ${\mathds{1}}_{l}\in\mathbb{R}^{l}$ and
$\mathbf{0}_{l}\in\mathbb{R}^{l}$, respectively, while
$I_{l}\in\mathbb{R}^{l\times l}$ refers to as the identity matrix. Let
$\Omega\in\mathbb{R}^{l\times l}$ be a square matrix. Relation $\Omega\succeq
0$ means that $\Omega$ is symmetric positive semidefinite. Notation
$[\Omega]_{ij}$ addresses the entry of matrix $\Omega$ at row $i$ and column
$j$, while $\Omega^{\top}$, $\mathrm{tr}(\Omega)$ and $\left\|\Omega\right\|$
denote respectively its transpose, its trace and its spectral norm. Operator
$\mathrm{col}_{l}[\Omega]$ represents the $l$-th column of $\Omega$ and its
$i$-th eigenvalue is denoted by $\lambda_{i}^{\Omega}$. The vector space
spanned by a vector $\omega\in\mathbb{R}^{l}$ is denoted with
$\left<{\omega}\right>$. The infinity norm of $\omega$ is identified by
$\left\|\omega\right\|_{\infty}$. Also,
$\boldsymbol{\omega}=\mathrm{vec}_{i=1}^{l}(\omega_{i})$ defines the
vectorization operator that stacks vectors $\omega_{i}$, $i=1,\dots,l$ as
$\boldsymbol{\omega}=\begin{bmatrix}\omega_{1}^{\top}&\dots&\omega_{l}^{\top}\end{bmatrix}^{\top}$;
whereas, $\mathrm{diag}_{i=1}^{l}(\varsigma_{i})$ is a block diagonal matrix
made up with $\varsigma_{i}\in\mathbb{R}$, $i=1,\dots,l$, on the diagonal and
$\mathrm{diag}(\boldsymbol{\omega})=\mathrm{diag}_{i=1}^{l}(\omega_{i})$.
Lastly, the Kronecker product is denoted with $\otimes$.

##  II Preliminary results

Consensus models for MASs and preliminary notions are here given along with a
brief overview of robustness results in consensus networks having multiple
uncertain edges. The SBDC protocol is then briefly recalled.

###  II-A Overview on uncertain consensus networks

An $n$-agent network can be modeled through a weighted graph
$\mathcal{G}=\left(\mathcal{V},\mathcal{E},\mathcal{W}\right)$ so that each
element in the vertex set $\mathcal{V}=\left\\{1,\dots,n\right\\}$ is
associated to an agent in the group, while the edge set
$\mathcal{E}=\\{e_{k}\\}_{k=1}^{m}\subseteq\mathcal{V}\times\mathcal{V}$
characterizes the agents’ information exchange. In addition,
$\mathcal{W}=\\{w_{k}\\}_{k=1}^{m}$, with $m=|\mathcal{E}|$, denotes the set
of weights attributed to each link. Throughout this work, bidirectional
interactions among agents are assumed; hence, $\mathcal{G}$ is an undirected
graph. A cycle is defined as a non-empty and distinct sequence of edges
joining a sequence of vertices, in which only the first and last vertices are
equal. If a graph does not contain cycles it is said acyclic and if it is also
connected it is called a tree. The neighborhood of node $i$ is defined as the
set
$\mathcal{N}_{i}=\left\\{j\in\mathcal{V}\setminus\\{i\\}\;|\;(i,j)\in\mathcal{E}\right\\}$,
while the degree of node $i$ is defined through the cardinality
$d_{i}=|\mathcal{N}_{i}|$ of neighborhood $\mathcal{N}_{i}$. Moreover, the
incidence matrix is denoted as $E\in\mathbb{R}^{n\times m}$, in which each
column $k\in\\{1,\ldots,m\\}$ is defined via the $k$-th (ordered) edge
$(i,j)\in\mathcal{E}$, where $i<j$ is chosen w.l.o.g., and for edge $k$
corresponding to $(i,j)$ one has $[E]_{lk}=-1$, if $l=i$; $[E]_{lk}=1$, if
$l=j$; $[E]_{lk}=0$, otherwise. For all $k=1,\ldots,m$, the weight
$w_{k}:=w_{ij}=w_{ji}\in\mathbb{R}$ is related to $k$-th edge $(i,j)$, and
$W=\mathrm{diag}_{k=1}^{m}(w_{k})$ is the diagonal matrix of edge weights.
Additionally, the Laplacian matrix, incorporating the topological information
about $\mathcal{G}$, is denoted as $L(\mathcal{G})=EWE^{\top}$ (see [15]).
Henceforward, we also suppose that graph $\mathcal{G}$ is connected and
$L(\mathcal{G})\succeq 0$, thus having eigenvalues $\lambda_{i}^{L}$, for
$i=1,\ldots,n$, such that
$0=\lambda_{1}^{L}<\lambda_{2}^{L}\leq\cdots\leq\lambda_{n}^{L}$. A sufficient
condition to satisfy the latter requisite, which is adopted throughout this
paper, is to take $w_{ij}>0$ for all $(i,j)$. With an appropriate labeling of
the edges, we can always assume that the incidence matrix
$E=\begin{bmatrix}E_{\mathcal{T}}&E_{\mathcal{C}}\end{bmatrix}$ can be
partitioned into the incidence matrix $E_{\mathcal{T}}$, relative to a
spanning tree $\mathcal{T}\subseteq\mathcal{G}$ with $\tau=n-1$ edges, and the
incidence matrix $E_{\mathcal{C}}$, associated to
$\mathcal{C}=\mathcal{G}\setminus\mathcal{T}$. Consequently, we define the
cut-set matrix of $\mathcal{G}$ (see [16]) as
$R_{(\mathcal{T},\mathcal{C})}=\begin{bmatrix}I_{\tau}&T_{(\mathcal{T},\mathcal{C})}\end{bmatrix}$,
with
$T_{(\mathcal{T},\mathcal{C})}=(E_{\mathcal{T}}^{\top}E_{\mathcal{T}})^{-1}E_{\mathcal{T}}^{\top}E_{\mathcal{C}}$.

A summary of the weighted consensus problem in MASs follows. Let us consider a
group of $n$ homogeneous agents, modeled by a weighted and connected graph
$\mathcal{G}$, and assign a continuous-time state
$x_{i}=x_{i}(t)\in\mathbb{R}^{D}$ to the $i$-th agent, for $i=1,\dots,n$. The
full network state is given by $\mathbf{x}=\mathrm{vec}_{i=1}^{n}(x_{i})\in
X\subseteq\mathbb{R}^{N}$, with $N=nD$. The weighted consensus for a MAS can
be characterized as follows.

######  Definition II.1 (Weighted Consensus [15])

An $n$-agent network achieves _consensus_ if
$\underset{t\rightarrow+\infty}{\lim}~{}\mathbf{x}(t)\in\mathcal{A}$, where,
for some $\omega\in\mathbb{R}^{D}$,
$\mathcal{A}=\left<{{\mathds{1}}_{n}}\right>\otimes\omega$ is termed the
agreement set.

For a connected graph $\mathcal{G}$ with positive weights, it is well known
that the linear weighted consensus protocol, given by

$\dot{\mathbf{x}}=-\mathbf{L}(\mathcal{G})\mathbf{x},$ (1)

where $\mathbf{L}(\mathcal{G})=(L(\mathcal{G})\otimes I_{D})$, satisfies
$\mathbf{x}(t)\in\mathcal{A}$ as $t\rightarrow+\infty$.

In this direction, we also revisit a robustness result for the consensus
protocol with small-magnitude perturbations on the edge weights [16]. Within
this framework, we take into consideration the perturbed Laplacian matrix
$L(\mathcal{G}_{\Delta^{W}})=E(W+\Delta^{W})E^{\top}$ for a structured norm-
bounded perturbation

$\displaystyle\Delta^{W}\in\boldsymbol{\Delta}^{W}=\\{\Delta^{W}\,:\,\Delta^{W}=\mathrm{diag}_{k=1}^{m}(\delta^{w}_{k}),\|\Delta^{W}\|\leq\bar{\delta}^{W}\\}.$
(2)

Letting
$\mathcal{E}_{\Delta}:=\\{e_{1}^{\Delta},e_{2}^{\Delta},\ldots\\}\subseteq\mathcal{E}$
be the (nonempty) subset of uncertain edges, we can define the matrix
$P\in\\{0,1\\}^{|\mathcal{E}|\times|\mathcal{E}_{\Delta}|}$ that selects the
uncertain edges in $\mathcal{E}$, with $[P]_{ij}=1$, if $i$ and $j$ satisfy
$e_{i}=e_{j}^{\Delta}$; and $[P]_{ij}=0$, otherwise. This leads to

###### Lemma II.1

Let
$\Delta^{W}\in\mathbb{R}^{|\mathcal{E}_{\Delta}|\times|\mathcal{E}_{\Delta}|}$,
with $\Delta^{W}$ diagonal, and consider the nominal weighted consensus
protocol (1). Then the perturbed consensus protocol

$\dot{\mathbf{x}}=-(L(\mathcal{G}_{\Delta^{W}})\otimes I_{D})\mathbf{x}$ (3)

achieves consensus $\forall\Delta^{W}\in\boldsymbol{\Delta}^{W}$ (defined in
(2)), if

$\left\|\Delta^{W}\right\|<{\underbrace{\left\|P^{\top}R_{(\mathcal{T},\mathcal{C})}^{\top}(R_{(\mathcal{T},\mathcal{C})}WR_{(\mathcal{T},\mathcal{C})}^{\top})^{-1}R_{(\mathcal{T},\mathcal{C})}P\right\|}_{=:\mathcal{R}_{\mathcal{E}_{\Delta}}(\mathcal{G})}}^{-1}.$
(4)

###### Proof:

By [16, Theorem V.2], the edge-agreement version of (3) is asymptotically
stable. Consequently, as $\mathcal{G}_{\Delta^{W}}$ is connected, one has
$L(\mathcal{G}_{\Delta^{W}})\succeq 0$ with a simple eigenvalue at $0$.
Therefore, Def. II.1 is satisfied
$\forall\Delta^{W}\in\boldsymbol{\Delta}^{W}$. ∎

In a slight abuse of convention, we refer to Lemma II.1 as _robust stability_
result; see [16] for more discussion on this notion. When
$|\mathcal{E}_{\Delta}|=1$, it was shown in [16] that the bound in (4) is
tight and that $\mathcal{R}_{\left\\{(u,v)\right\\}}(\mathcal{G})$ can be
interpreted as the effective resistance between a pair of nodes $(u,v)$.
However, for multi-edge attacks this bound is inherently conservative; this
aspect will be elaborated upon in Section III-C.

###  II-B Secure-by-design consensus dynamics

In this work, we consider MASs that are subject to the same key principles
assumed in [1], that is with presence of tasks, objective coding, information
localization and a network manager111The network manager does not govern the
agents’ dynamics, namely it should not be intended to fulfill the role of a
global controller. Instead, it is chosen to precisely serve as an encryption
mechanism to set up and secure distributed algorithms running on the
underlying MAS.. These elements, which comprise the basic setup of the SBCD
dynamics, are recalled in the following lines.

Tasks are described by an encoded parameter $\theta$ that we term the codeword
and the space of all tasks is denoted as $\Theta$. Each agent in the network
then decodes this objective using its objective decoding function, defined as
$p_{i}:\Theta\to\Pi_{i}$, where $\Pi_{i}$ depends on the particular
application (e.g. $\Pi_{i}\subseteq\mathbb{R}^{n}$ within the consensus
setting). For $\theta\in\Theta$, $p_{i}(\theta)$ is called the _localized
objective_. Instead, if $\theta\notin\Theta$, $p_{i}(\theta)$ may not be
computable; nonetheless, any agent collecting such a codeword may send an
alert. More precisely, the objective coding is established via the non-
constant functions
$p_{i}(\theta):\Theta\rightarrow\Pi_{i}\subseteq\mathbb{R}^{n}$, such that
$[p_{i}(\theta)]_{j}:=p_{ij}(\theta)$ with

$p_{ij}(\theta)=\begin{cases}w_{ij},\quad\text{ if }(i,j)\in\mathcal{E};\\\
0,\quad~{}~{}~{}\\!\text{ otherwise}.\end{cases}$ (5)

The values $w_{ij}$ in (5) coincide with the nominal desired consensus weights
that are assigned, encoded and broadcast by the network manager, to achieve
desired convergence rates or other performance metrics in consensus networks.
Moreover, the information localization about the global state $\mathbf{x}$ is
expressed by means of $h_{i}(\mathbf{x}):X\rightarrow
Y_{i}\subseteq\mathbb{R}^{D\times n}$, such that
$\mathrm{col}_{j}[h_{i}(\mathbf{x})]:=h_{ij}(\mathbf{x}(t))\in\mathbb{R}^{D}$
with $h_{ij}(\mathbf{x})=x_{i}-x_{j}$, if $(i,j)\in\mathcal{E}$;
$h_{ij}(\mathbf{x})=\mathbf{0}_{D}$, otherwise. So, the $i$-th agent’s
dynamics for the SBDC is determined by

$\displaystyle\dot{x}_{i}=-{\textstyle\sum}_{j\in\mathcal{N}_{i}}p_{ij}(\theta)h_{ij}(\mathbf{x}).$
(6)

Remarkably, (6) reproduces exactly the linear consensus protocol already
introduced in (1) and, defining
$\mathbf{p}(\theta)=\mathrm{vec}_{i=1}^{n}(p_{i}(\theta))\in\mathbb{R}^{n^{2}}$
and
$\mathbf{H}(\mathbf{x})=\mathrm{diag}_{i=1}^{n}(h_{i}(\mathbf{x}(t)))\in\mathbb{R}^{N\times
n^{2}}$, dynamics (6) can be rewritten as

$\dot{\mathbf{x}}=-\mathbf{H}(\mathbf{x})\mathbf{p}(\theta),$ (7)

leading to the following theoretical result.

######  Lemma II.2 ([1, Lemma III.1])

The SBDC protocol (7) reaches agreement for any given objective decoding
function $\mathbf{p}$ that satisfies (5).

Likewise, we consider and investigate the well-known discrete-time consensus
dynamics described by

$\mathbf{x}(t+1)=\mathbf{x}(t)-\epsilon\mathbf{L}(\mathcal{G})\mathbf{x}(t),$
(8)

where $\epsilon$ is a common parameter shared among all agents and designed to
belong to the interval $(0,2/\lambda_{n}^{L})$, as shown in [17]. According to
the characterization in (7), the SBCD dynamics adopted in discrete time in (8)
is also given by

$\mathbf{x}(t+1)=\mathbf{x}(t)-\epsilon\mathbf{H}(\mathbf{x}(t))\mathbf{p}(\theta),$
(9)

since
$\mathbf{H}(\mathbf{x})\mathbf{p}(\theta)=\mathbf{L}(\mathcal{G})\mathbf{x}$
holds thanks to Lem. II.2.

##  III Robustness of the SBDC protocol to structured channel tampering

The original contribution provided by this study aims at the design of secure
network systems to structured channel tampering while achieving consensus
task. To this aim, the system is embedded with security measures that allow to
render the network robust to small signal perturbations on some of the links.
Also, a description for the structured channel tampering is given along with
the relative robustness analysis for the SBDC protocol under multiple threats.

###  III-A Models, problem statement and key assumptions

The structured channel tampering problem under analysis is formulated as
follows. Similarly to the model adopted in [1], the prescribed codeword
$\theta$ is subject to a perturbation
$\delta^{\theta}\in\boldsymbol{\Delta}^{\theta}=\\{\delta^{\theta}\,:\,\|\delta^{\theta}\|_{\infty}\leq\bar{\delta}^{\theta}\\}$.
We let $\Theta$ be a Euclidean subspace, in particular
$\Theta=\Theta_{11}\times\Theta_{12}\times\cdots\times\Theta_{nn}\subseteq\mathbb{R}^{n^{2}}$,
and allow a codeword $\theta=\mathrm{vec}_{i=1}^{n}(\theta_{i})\in\Theta$ to
be divided into (at most) $n(n-1)/2$ relevant “subcodewords”
$\theta^{(k)}:=[\theta_{i}]_{j}=\theta_{ij}$, with $k=1,\ldots,m$, such that
$\theta_{ij}=\theta_{ji}$, if $i\neq j$, and $\theta_{ii}$ is free to vary,
for $i=1,\ldots,n$. Each $\theta_{ij}\in\Theta_{ij}\subseteq\mathbb{R}$ can be
seen as the $j$-th component of the $i$-th codeword fragment $\theta_{i}$,
with $i=1,\ldots,n$. Such subcodewords influence the value of $p_{ij}(\theta)$
directly if and only if $j\in\mathcal{N}_{i}$, i.e., it holds that
$p_{ij}(\theta)=p_{ij}(\theta_{ij})$ for all $(i,j)\in\mathcal{E}$.

Therefore, the consensus description to support this investigation is such
that the $i$-th nominal dynamics in (7) is modified as

$\dot{x}_{i}=-{\textstyle\sum}_{j\in\mathcal{N}_{i}}p_{ij}(\theta_{ij}+\delta^{\theta}_{ij})h_{ij}(\mathbf{x}),\quad
i=1,\ldots,n,$ (10)

with $\delta^{\theta}_{ij}=[\delta^{\theta}_{i}]_{j}$ and
$\delta^{\theta}_{i}$ satisfying
$\delta^{\theta}=\mathrm{vec}_{i=1}^{n}(\delta^{\theta}_{i})$. Analogously,
the $i$-th perturbed discrete time dynamics in (9) can be written as:

$x_{i}(t+1)=x_{i}(t)-\epsilon{\textstyle\sum}_{j\in\mathcal{N}_{i}}p_{ij}(\theta_{ij}+\delta^{\theta}_{ij})h_{ij}(\mathbf{x}(t)),$
(11)

where $\epsilon>0$ must be selected. In light of the previous discussion, the
following design problem can be now stated.

###### Problem III.1

Design objective functions $p_{ij}$ such that (10) (resp., (11) in discrete
time) reaches consensus, independently from the codeword
$\theta\in\Theta\subseteq\mathbb{R}^{n^{2}}$, while the underlying MAS is
subject to a structured attack
$\delta^{\theta}\in\boldsymbol{\Delta}^{\theta}$ on multiple edges (belonging
to $\mathcal{E}_{\Delta}$), i.e., with $\delta_{ij}^{\theta}=0$ for all
$(i,j)\in\mathcal{E}\setminus\mathcal{E}_{\Delta}$. In addition, determine the
largest $\rho^{\theta}_{\Delta}$ ensuring robust stability of (10) (resp.,
(11) in discrete time).

Within this framework, it is possible to leverage Lem. II.1 and yield the main
theoretical contribution of this paper, represented by the robustness
guarantees for system (10) when the target of a cyber-physical attack is a
multitude of edges. In this direction, we study how the robust stability of
(10) is affected by perturbations on all the weights
$p_{uv}(\theta_{uv})=w_{uv}$ attached to the connections
$(u,v)\in\mathcal{E}_{\Delta}$ that are caused by the deviations of each
subcodeword $\theta_{uv}$.

As clarified later in more detail, the same three additional assumptions on
the $p_{i}$’s proposed in [1] are sufficient to tackle Problem III.1, namely,
this robustness analysis is again restricted to a particular choice for the
objective coding, that is for concave and Lipschitz continuous differentiable
functions $p_{i}$. Thus, we let the $i$-th objective decoding function adopted
in model (10) have the following properties.

###### Assumption III.1

Each $p_{i}:\Theta\rightarrow\Pi_{i}$, with $i=1,\ldots,n$, has the subsequent
characterization:

  1. (i)

values $[p_{i}(\theta)]_{j}=p_{ij}(\theta_{ij})$, with
$\theta_{ij}=[\theta_{i}]_{j}$, satisfy (5) for all $(i,j)\in\mathcal{E}$ and
are not constant w.r.t. $\theta_{ij}$;

  2. (ii)

$p_{ij}$ is concave $\forall\theta\in\Theta$, i.e.,
$p_{ij}(\varsigma\eta_{1}+(1-\varsigma)\eta_{2})\geq\varsigma
p_{ij}(\eta_{1})+(1-\varsigma)p_{ij}(\eta_{2})$, $\varsigma\in[0,1]$,
$\forall\eta_{1},\eta_{2}\in\Theta_{ij}$;

  3. (iii)

$p_{ij}$ is Lipschitz continuous and differentiable w.r.t. $\theta$, implying
$\exists K_{ij}\geq 0:~{}|p_{ij}^{\prime}(\theta_{ij})|\leq K_{ij}$,
$\forall(i,j)\in\mathcal{E}$.

Also, to provide analytical guarantees to multiple attacks striking the
network, we will make use of the global quantity

$K_{\Delta}=\max_{(u,v)\in\mathcal{E}_{\Delta}}\left\\{K_{uv}\right\\}.$ (12)

###  III-B Guarantees for multiple attacks

The guarantees in [1, Theorem IV.1] can be extended as follows for a
continuous-time multiple-attack scenarios.

###### Theorem III.1

Assume that the characterization for objective decoding functions $p_{i}$ in
Asm. III.1 holds. For a structured injection attack
$\delta^{\theta}\in\boldsymbol{\Delta}^{\theta}$ affecting all edges in
$\mathcal{E}_{\Delta}$ define
$\mathcal{R}_{\mathcal{E}_{\Delta}}(\mathcal{G})$ and $K_{\Delta}$ as in (4)
and (12), respectively. Then the perturbed consensus protocol (10) is stable
and achieves agreement for all $\delta^{\theta}$ whenever

$\left\|\delta^{\theta}\right\|_{\infty}<\rho^{\theta}_{\Delta}=(K_{\Delta}\mathcal{R}_{\mathcal{E}_{\Delta}}(\mathcal{G}))^{-1},$
(13)

independently from the values taken by any codeword $\theta\in\Theta$.

###### Proof:

Similarly to the single edge attack case, Asm. III.1 brings to each ordered
logical step to conclude the thesis through Lem. II.1. Indeed, (i)-(iii) lead
to the determination of quantity $K_{ij}|\delta^{\theta}_{ij}|$, which can be
seen as the maximum magnitude of an additive perturbation
$\delta^{w}_{ij}:=p_{ij}(\theta_{ij}+\delta^{\theta}_{ij})-p_{ij}(\theta_{ij})$
affecting each $p_{ij}$, $\forall(i,j)\in\mathcal{E}$, independently from the
transmitted codeword $\theta$. Consequently, the fact that
$|\delta_{uv}^{w}|\leq K_{uv}|\delta_{uv}^{\theta}|$ holds for each edge
$(u,v)\in\mathcal{E}_{\Delta}$ implies that the following chain of
inequalities is verified:

$\left\|\delta^{w}\right\|_{\infty}\leq\max_{(u,v)\in\mathcal{E}_{\Delta}}\left\\{K_{uv}|\delta^{\theta}_{uv}|\right\\}\leq
K_{\Delta}\left\|\delta^{\theta}\right\|_{\infty}.$ (14)

Therefore, imposing inequality
$K_{\Delta}\left\|\delta^{\theta}\right\|_{\infty}<\mathcal{R}^{-1}_{\mathcal{E}_{\Delta}}(\mathcal{G})$,
in accordance with (4), leads to the thesis. ∎

We notice that the leftmost inequality in (14) translates into an essential
conservatism on the attack magnitude w.r.t. (4), similarly to the case
$|\mathcal{E}_{\Delta}|=1$. This can be modulated by choosing constants
$K_{ij}$ properly. Nonetheless, as soon as $|\mathcal{E}_{\Delta}|>1$ holds,
another essential conservatism that hinges on the attack’s scale
$|\mathcal{E}_{\Delta}|$ may arise w.r.t. its one-dimensional version. In the
sequel, this concept is referred to as the $\mathcal{E}_{\Delta}$-gap, or
resilience gap, since it depends to the fact that, in general, one has
$\mathcal{R}_{\mathcal{E}_{\Delta}}(\mathcal{G})\geq\max_{(u,v)\in\mathcal{E}_{\Delta}}\left\\{\mathcal{R}_{\left\\{(u,v)\right\\}}(\mathcal{G})\right\\}$
(see discussion in Sec. III-C for more details).

Concerning, instead, the discrete-time guarantees for system (11) provided in
[1, Theorem VI.1, Corollary VI.1], the following generalization can be made.

###### Corollary III.1

Assume that the characterization for objective decoding functions $p_{i}$ in
Asm. III.1 holds. Denote respectively with
$\bar{w}_{i}=\sum_{j\in\mathcal{N}_{i}}|w_{ij}|$ and
$\Psi_{\mathcal{G}}=\max_{i=1,\ldots,n}\left\\{\bar{w}_{i}\right\\}$ the
weighted degree of the $i$-th node and the maximum weighted degree of the
underlying graph $\mathcal{G}$. Let a structured injection attack
$\delta^{\theta}\in\boldsymbol{\Delta}^{\theta}$ affect all edges in
$\mathcal{E}_{\Delta}$ and define

$\displaystyle\psi_{i}(\delta^{\theta}_{uv})$
$\displaystyle=\bar{w}_{i}+K_{uv}|\delta^{\theta}_{uv}|,\quad\forall(u,v)\in\mathcal{E}_{\Delta},~{}i=u,v.$
(15)

Then the perturbed consensus protocol (11) is stable for all $\delta^{\theta}$
such that both (13) and

$\phi_{\mathcal{G}}(\delta^{\theta}):=\max\left\\{\Psi_{\mathcal{G}},\max_{i\in\left\\{u,v\right\\},(u,v)\in\mathcal{E}_{\Delta}}\psi_{i}(\delta^{\theta}_{uv})\right\\}<\epsilon^{-1}$
(16)

hold for any fixed $\epsilon$, independently from the values taken by any
codeword $\theta\in\Theta$.

###### Proof:

The proof is a generalization of that for Theorem VI.1 in [1], with the only
difference that $\phi_{\mathcal{G}}(\delta^{\theta})$ in (16) is harder to be
computed w.r.t. its single-edge-attack version
$\phi_{\mathcal{G}}(\delta^{\theta}_{uv})$ because the combinatorial search
space for the maximization of quantities $\psi_{i}(\delta^{\theta}_{uv})$
defined in (15) is, in general, larger in this kind of multiple-attack
scenario. ∎

###### Corollary III.2

Under all the assumptions adopted in Cor. III.1, defining $K_{\Delta}$ as in
(12), $\mathcal{R}_{\mathcal{E}_{\Delta}}(\mathcal{G})$ as in (4) and setting
$\epsilon<\Psi_{\mathcal{G}}^{-1}$, the perturbed consensus protocol (11) is
stable for all $\delta^{\theta}$ such that

$\left\|\delta^{\theta}\right\|_{\infty}<\rho^{\theta}_{\Delta}=K_{\Delta}^{-1}\min\\{\mathcal{R}_{\mathcal{E}_{\Delta}}^{-1}(\mathcal{G}),(\epsilon^{-1}-\Psi_{\mathcal{G}})\\}$
(17)

independently from the values taken by any codeword $\theta\in\Theta$. In
particular, if $\epsilon$ is selected as follows

$\epsilon\leq\epsilon_{\Delta}^{\star}:=(\Psi_{\mathcal{G}}+\mathcal{R}_{\mathcal{E}_{\Delta}}^{-1}(\mathcal{G}))^{-1}$
(18)

then $\rho^{\theta}_{\Delta}$ in (17) is maximized as $\epsilon$ varies and
condition (13) needs to be fulfilled solely to guarantee robust stability.

###### Proof:

Relation in (17) is the combined result of guarantee in (13) and that one
obtainable imposing
$\Psi_{\mathcal{G}}+K_{\Delta}\left\|\delta^{\theta}\right\|_{\infty}<\epsilon^{-1}$
to satisfy (16), since $\phi_{\mathcal{G}}(\delta^{\theta})$ can be upper
bounded as
$\phi_{\mathcal{G}}(\delta^{\theta})\leq\Psi_{\mathcal{G}}+K_{\Delta}\left\|\delta^{\theta}\right\|_{\infty}$.
On the other hand, relation (18) is derived enforcing
$\mathcal{R}_{\mathcal{E}_{\Delta}}^{-1}(\mathcal{G})\leq\epsilon^{-1}-\Psi_{\mathcal{G}}$
with the purpose to maximize $\rho_{\Delta}^{\theta}$ as $\epsilon$ varies. ∎

We conclude this subsection with the following remark.

###### Remark III.1

Observe that no additional assumptions on the decoding functions $p_{i}$
w.r.t. (i)-(iii) given in Subsec. III-A are required to solve Problem III.1,
i.e., to generalize the guarantees already yielded in [1] to a multiple-attack
scenario.

###  III-C Analysis of the resilience gap

In light of the previous theoretical results, we discuss here how decoding
functions can be seen as a useful tool to compensate against the resilience
gap (the $\mathcal{E}_{\Delta}$-gap). This analysis starts by recalling the
following preliminary proposition.

######  Proposition III.1 ([16, Proposition V.3])

For any weighted undirected graph $\mathcal{G}$ it holds that

$\mathcal{R}_{\mathcal{E}_{\Delta}}^{\star}(\mathcal{G})\leq\mathcal{R}_{\mathcal{E}_{\Delta}}(\mathcal{G})\leq\mathcal{R}_{\mathcal{E}_{\Delta}}^{tot}(\mathcal{G}),$
(19)

where $\mathcal{R}_{\mathcal{E}_{\Delta}}(\mathcal{G})$ is defined as in (4),

$\mathcal{R}_{\mathcal{E}_{\Delta}}^{\star}(\mathcal{G})=\underset{(u,v)\in\mathcal{E}_{\Delta}}{\max}\left\\{\mathcal{R}_{\left\\{(u,v)\right\\}}(\mathcal{G})\right\\}$
(20)

and
$\mathcal{R}_{\mathcal{E}_{\Delta}}^{tot}(\mathcal{G})=\mathrm{tr}(P^{\top}R_{(\mathcal{T},\mathcal{C})}^{\top}(R_{(\mathcal{T},\mathcal{C})}WR_{(\mathcal{T},\mathcal{C})}^{\top})^{-1}R_{(\mathcal{T},\mathcal{C})}P)$.

Prop. III.1 suggests us to define the $\mathcal{E}_{\Delta}$-gap as

$g(\mathcal{G},\mathcal{E}_{\Delta})=1-\mathcal{R}_{\mathcal{E}_{\Delta}}^{\star}(\mathcal{G})/\mathcal{R}_{\mathcal{E}_{\Delta}}(\mathcal{G}),$
(21)

with $\mathcal{R}_{\mathcal{E}_{\Delta}}^{\star}(\mathcal{G})$ and
$\mathcal{R}_{\mathcal{E}_{\Delta}}(\mathcal{G})$ defined by (20) and (4).
Indeed, the emerging conservatism related to the fact that multiple edges may
be under attack ($|\mathcal{E}_{\Delta}|>1$) grows as the value of
$g(\mathcal{G},\mathcal{E}_{\Delta})\in[0,1)$ increases. Also, exploiting
(21), inequality (13) can be rewritten as

$\left\|\delta^{\theta}\right\|_{\infty}<(1-g(\mathcal{G},\mathcal{E}_{\Delta}))/(K_{\Delta}\mathcal{R}^{\star}_{\mathcal{E}_{\Delta}}(\mathcal{G})),$
(22)

so that the quantity
$\rho_{\Delta^{\star}}^{\theta}:=(K_{\Delta}\mathcal{R}^{\star}_{\mathcal{E}_{\Delta}}(\mathcal{G}))^{-1}$
can be seen as the maximum value taken by $\rho_{\Delta}^{\theta}$ in (13) as
$\mathcal{G}$ and $\mathcal{E}_{\Delta}$ vary. In particular, if
$g(\mathcal{G},\mathcal{E}_{\Delta})=0$ holds (indicating the absence of this
kind of conservatism), there exists an edge
$(u^{\star},v^{\star})\in\mathcal{E}_{\Delta}$ for which, by (22), the network
$\mathcal{G}$ is robustly stable for all perturbations $\delta^{\theta}$ such
that

$\left\|\delta^{\theta}\right\|_{\infty}<(K_{\Delta}\mathcal{R}_{\left\\{(u^{\star},v^{\star})\right\\}}(\mathcal{G}))^{-1}=\rho_{\Delta^{\star}}^{\theta}.$
(23)

Remarkably, inequality (23) expresses the same guarantee provided by (13) as
if $\mathcal{E}_{\Delta}=\left\\{(u^{\star},v^{\star})\right\\}$ was assigned,
similarly to the one-dimensional case $|\mathcal{E}_{\Delta}|=1$ debated in
[1]. Since the case where only one edge is attacked clearly represents a
scenario in which the conservatism due to $|\mathcal{E}_{\Delta}|>1$ is
lacking, if the $\mathcal{E}_{\Delta}$-gap corresponding to a given setup
$(\mathcal{G},\mathcal{E}_{\Delta})$ vanishes then the robustness of
$\mathcal{G}$ subject to perturbations $\delta^{\theta}$ striking the subset
$\mathcal{E}_{\Delta}$ is maximized. The latter observation leads us to wonder
which are all the possible circumstances where
$g(\mathcal{G},\mathcal{E}_{\Delta})=0$ is satisfied? A partial yet extensive
answer is given by the next proposition.

###### Proposition III.2

Let $\mathcal{G}=(\mathcal{V},\mathcal{E},\mathcal{W})$,
$\mathcal{E}_{\Delta}\subseteq\mathcal{E}$ the subset of perturbed edges in
$\mathcal{E}$, and $g(\mathcal{G},\mathcal{E}_{\Delta})$ the
$\mathcal{E}_{\Delta}$-gap (21). Then $g(\mathcal{G},\mathcal{E}_{\Delta})=0$
if either one of the following conditions hold: i)
$~{}|\mathcal{E}_{\Delta}|=1$, or ii) $~{}2\leq|\mathcal{E}_{\Delta}|\leq
n-1=|\mathcal{E}|$.

###### Proof:

We prove each point of the statement separately.

i) If $|\mathcal{E}_{\Delta}|=1$ then there exists only one edge
$(u^{\star},v^{\star})\in\mathcal{E}_{\Delta}$ for which the mid inequality in
(19) must hold tightly, since
$\mathcal{R}^{\star}_{\mathcal{E}_{\Delta}}(\mathcal{G})=\mathcal{R}_{\mathcal{E}_{\Delta}}(\mathcal{G})$.
It follows that $g(\mathcal{G},\mathcal{E}_{\Delta})=0$.

ii) In this case, since $\mathcal{G}$ is connected by assumption then it has
to be a spanning tree; consequently, it can be also denoted as
$\mathcal{G}=\mathcal{T}(\mathcal{E},\mathcal{V},\mathcal{W})$. Hence, one has
$R_{(\mathcal{T},\mathcal{C})}=I_{n-1}$ and, therefore, formula (4) yields

$\mathcal{R}_{\mathcal{E}_{\Delta}}(\mathcal{T})=\underset{k:~{}e_{k}\in\mathcal{E}_{\Delta}}{\max}\left\\{|w_{k}^{-1}|\right\\},\quad\forall\mathcal{E}_{\Delta}\subseteq\mathcal{E}.$
(24)

Observe that (24) also holds for subsets $\mathcal{E}_{\Delta}$ such that
$|\mathcal{E}_{\Delta}|=1$ and the $\max$ over
$\left\\{k:~{}e_{k}\in\mathcal{E}_{\Delta}\right\\}$ in (24) is taken exactly
as in (20). We conclude that $g(\mathcal{G},\mathcal{E}_{\Delta})=0$ occurs
even for $|\mathcal{E}_{\Delta}|\geq 2$ whenever
$\mathcal{G}=\mathcal{T}(\mathcal{E},\mathcal{V},\mathcal{W})$, as the mid
inequality in (19) is tightly satisfied under these conditions. ∎

To conclude, we note that the $\mathcal{E}_{\Delta}$-gap in (22) can be
mitigated by an appropriate re-design procedure of $K_{\Delta}$.

###### Proposition III.3

Denote with $\rho_{\Delta^{\star}}^{\theta}(K_{\Delta})$ the value of
$\rho_{\Delta^{\star}}^{\theta}$ as $K_{\Delta}$ varies. Set
$K_{\Delta}^{\prime}:=(1-g(\mathcal{G},\mathcal{E}_{\Delta}))K_{\Delta}$. Then

$\rho_{\Delta^{\star}}^{\theta}(K_{\Delta}^{\prime})>\rho_{\Delta^{\star}}^{\theta}(K_{\Delta}),\quad\forall
K_{\Delta}>0.$

Prop. III.3 introduces an interesting tradeoff, where conservatism is reduced
at the expense of shrinking the image of the decoding function (see also Fact
IV.1 in [1]). This leads to an increasing demand of the encryption
capabilities.

##  IV Numerical examples

We here focus on a potential application of the proposed technique to semi-
autonomous networks (SANs) consisting of leader-follower autonomous agents
[18]. Let us assume that subset $\mathcal{V}_{l}\subseteq\mathcal{V}$,
$\mathcal{V}_{l}\neq\varnothing$, collects all the leader agents of a SAN
$\mathcal{S}=(\mathcal{G},\mathcal{U})$, where
$\mathcal{G}=(\mathcal{V},\mathcal{E},\mathcal{W})$ is an undirected and
connected graph representing the network interactions of $\mathcal{S}$ and
$\mathcal{U}=\left\\{\mathbf{u}_{1}(t),\ldots,\mathbf{u}_{|\mathcal{V}_{l}|}(t)\right\\}$,
$\mathbf{u}_{\ell}(t)\in\mathbb{R}^{D}$, denotes the set of external inputs
that directly influence each leader agents in $\mathcal{V}_{l}$. Typically,
setting $\mathbf{u}(t)=\mathrm{vec}_{i\in\mathcal{V}_{l}}(\mathbf{u}_{i}(t))$,
the continuous-time dynamics for SANs are yielded by

$\dot{\mathbf{x}}=(L_{B}(\mathcal{G})\otimes I_{D})\mathbf{x}+(B\otimes
I_{D})\mathbf{u},$ (25)

where
$L_{B}(\mathcal{G})=L(\mathcal{G})+\mathrm{diag}(B{\mathds{1}}_{|\mathcal{V}_{l}|})$,
with $B\in\mathbb{R}^{n\times|\mathcal{V}_{l}|}$ such that $[B]_{i\ell}>0$, if
agent $i$ belongs to the leader set $\mathcal{V}_{l}$; $[B]_{i\ell}=0$,
otherwise. Also, a discrete-time version of (25) can be evidently provided by

$\mathbf{x}(t+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1})=(I_{N}-\epsilon(L_{B}(\mathcal{G})\otimes
I_{D}))\mathbf{x}(t)+(B\otimes I_{D})\mathbf{u}(t),$ (26)

where the quantity $\epsilon$ preserves the same meaning of (11). The
stability of SANs endowed with the discussed security mechanisms can be
analyzed by observing that the positive semi-definiteness of $L(\mathcal{G})$
clearly implies the positive semi-definiteness of $L_{B}(\mathcal{G})$.
Consequently, it is sufficient to ensure Thm. III.1 and Cor. III.1 to hold in
order to respectively guarantee the robust stability of protocols (25) and
(26).

Figure 1: (a): underlying topology of the given SAN (leaders in orange and
followers in light blue); (b): Objective decoding functions considered for
this case study.

In the following lines, we examine the SAN
$\mathcal{S}=(\mathcal{G},\mathcal{U})$ whose interconnections are modeled
through the graph $\mathcal{G}=(\mathcal{V},\mathcal{E},\mathcal{W})$
illustrated in Fig. 1: the sole leader (node $1$) is marked in orange, the
followers (remaining nodes) are marked in light blue, and the edge weights
$\mathcal{W}=\left\\{w_{k}\right\\}_{k=1}^{5}$ given by $w_{1}=w_{12}$,
$w_{2}=w_{35}$, $w_{3}=w_{46}$, $w_{4}=w_{24}$, $w_{5}=w_{23}$ are chosen
according to the same picture. Moreover, we set $D=1$, $[B]_{i\ell}=1$ for
$\ell=1$ if $i=1$, and $\mathcal{U}=\left\\{u_{1}\right\\}$, with
$u_{1}=-0.5$. We test the proposed security methods both in continuous and
discrete time by uniformly adopting for all the edges the following objective
decoding function

$p_{\Delta}(\eta)=\begin{cases}K_{\Delta}\left(\frac{4}{13}\sqrt{\eta+1}+1\right),&\quad\text{if
}\eta\geq 3;\\\
K_{\Delta}\left(-\frac{2}{13}\eta^{2}+\eta\right),&\quad\text{if
}0\leq\eta<3;\\\ K_{\Delta}\eta,&\quad\text{if }\eta<0;\end{cases}$ (27)

where $K_{\Delta}>0$ is a suitable Lipschitz constant for $p_{\Delta}(\eta)$,
to be selected in order to ensure that
$\left\|\delta^{w}\right\|_{\infty}\leq\rho_{\Delta}^{\theta}=0.5$. More in
detail, we consider two kinds of malicious attacks hitting part of
$\mathcal{W}$. The first external perturbation is delivered against edge
$(1,2)$ (red strike in Fig. 1). Whereas, the second one involves not only edge
$(1,2)$ but also edges $(3,5)$ and $(4,6)$ (yellow strikes in Fig. 1). Hence,
under this specific adversarial setup, we denote the set of edges under attack
as $\mathcal{E}_{1}=\left\\{(1,2)\right\\}$ and
$\mathcal{E}_{2}=\left\\{(1,2),(3,5),(4,6)\right\\}$ for the first and second
kind of perturbation, respectively.
Observing that $\mathcal{G}$ in Fig. 1 is an undirected tree (i.e.
$\mathcal{G}=\mathcal{T}$) then
$\mathcal{R}_{\mathcal{E}_{\Delta}}(\mathcal{G})$ can be computed as in (24).
One can set
$\epsilon=(\mathcal{R}_{\mathcal{E}_{1}}^{-1}(\mathcal{G})+\Psi_{\mathcal{G}})^{-1}=(2w_{1}+w_{4}+w_{5})^{-1}=0.1$
in order to satisfy (18) for both $\mathcal{E}_{1}$, $\mathcal{E}_{2}$; thus,
the Lipschitz constant $K_{\Delta}$ of $p_{\Delta}$ can be finally designed
according to (24) and (13) so that it respectively takes values $K_{1}=6$ and
$K_{2}=2$ for perturbations on edges $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$.
In Fig. 1, the decoding function (27) is then depicted (solid curves) for such
values of $K_{\Delta}$ along with the corresponding (dashed) lines having
slope equal to its Lipschitz constant.

We now discuss the robust stability of the SAN under investigation. In
particular, we focus on the trajectories of (25) and (26) for $\mathcal{S}$
subject to the following perturbations:
$\delta_{1}^{\theta}=-0.5\rho_{\Delta}^{\theta}[1~{}\\!\\!0~{}\\!\\!0~{}\\!\\!0~{}\\!\\!0]^{\top}$,
$\delta_{2}^{\theta}=-0.5\rho_{\Delta}^{\theta}[1~{}\\!\\!1~{}\\!\\!1~{}\\!\\!0~{}\\!\\!0]^{\top}$,
where $\delta_{1}^{\theta}$ strikes the edge in $\mathcal{E}_{1}$ and
$\delta_{2}^{\theta}$ strikes the three edges in $\mathcal{E}_{2}$. It is
crucial to note that if $K_{\Delta}\\!\\!=\\!\\!K_{1}$ is selected as in the
single-case attack scenario proposed in [1] adopting $p_{\Delta}$ to counter
the perturbation $\delta_{2}^{\theta}$, then no stability guarantees can be
given (see Figs. 2(a)-2(b)). Indeed, despite this choice allows to mitigate
the malicious effects of $\delta_{1}^{\theta}$, the additional perturbations
on edges $(3,5)$ and $(4,6)$ are capable of disrupting the network state
convergence towards $u_{1}$ as $t\\!\rightarrow\\!+\infty$ because
$K_{\Delta}=K_{1}$ does not match condition (13). Instead, Figs. 2(c)-2(d)
show that such an issue is overtaken if $K_{\Delta}=K_{2}$ is set ensuring
proper state convergence in accordance with (13).

(a) Divergence (continuous time).

(b) Divergence (discrete time).

(c) Convergence (continuous time).

(d) Convergence (discrete time).

Figure 2: Behaviors of the SAN trajectories as $\mathcal{E}_{\Delta}$ and
$p_{\Delta}$ vary. According to $\epsilon$, each discrete time stamp is set to
$0.1~{}s$. (a)-(b): stability guarantees for attack $\delta_{1}^{\theta}$ on
$\mathcal{E}_{1}$ but not for $\delta_{2}^{\theta}$ on $\mathcal{E}_{2}$;
(c)-(d): stability guarantees for the attacks $\delta_{i}^{\theta}$ on
$\mathcal{E}_{i}$, $i=1,2$.

In light of this, we can appreciate that the extension of the SBDC protocol
towards multiple-attack scenarios requires, in general, more encryption
capabilities and attention to the design phase in order to ensure robust
stability compared to the single-attack case. Indeed, in this simulation, we
have seen that only by considering
$\mathcal{R}_{\mathcal{E}_{\Delta}}(\mathcal{G})$ computed as in (4) and
lowering the value for $K_{\Delta}$ the agreement can be reached.

##  V Final remarks and future directions

The paper broadens the approach devised in [1] for secure consensus networks
to a multiple attack scenario, both in the continuous and discrete time
domains. Even in this framework, small-gain-theorem-based stability guarantees
are yielded to this aim and, remarkably, no additional assumption is needed to
provide such a generalization. In addition, the conservatism arising from a
multiplicity of threats has been addressed and analyzed. Future works will
involve new developments towards other multi-agent protocols, such as
nonlinear consensus and distance-based formation control.

## References

  * [1] M. Fabris and D. Zelazo, “Secure consensus via objective coding: Robustness analysis to channel tampering,” _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ , pp. 1–13, 2022.
  * [2] R. Olfati-Saber, J. A. Fax, and R. M. Murray, “Consensus and cooperation in networked multi-agent systems,” _Proceedings of the IEEE_ , vol. 95, no. 1, pp. 215–233, Jan 2007.
  * [3] M. Meng, G. Xiao, and B. Li, “Adaptive consensus for heterogeneous multi-agent systems under sensor and actuator attacks,” _Automatica_ , vol. 122, p. 109242, 2020.
  * [4] J. Wang, X. Deng, J. Guo, and Z. Zeng, “Resilient consensus control for multi-agent systems: A comparative survey,” _Sensors_ , vol. 23, no. 6, 2023\.
  * [5] D. Zhang, G. Feng, Y. Shi, and D. Srinivasan, “Physical safety and cyber security analysis of multi-agent systems: A survey of recent advances,” _IEEE/CAA Journal of Automatica Sinica_ , vol. 8, no. 2, pp. 319–333, 2021\.
  * [6] Y. Yang, Y. Xiao, and T. Li, “Attacks on formation control for multiagent systems,” _IEEE Transactions on Cybernetics_ , pp. 1–13, 2021.
  * [7] S. Huo, H. Wu, and Y. Zhang, “Secure consensus control for multi-agent systems against attacks on actuators and sensors,” _International Journal of Robust and Nonlinear Control_ , vol. 32, no. 8, pp. 4861–4877, 2022.
  * [8] Z.-W. Liu, Y.-L. Shi, H. Yan, B.-X. Han, and Z.-H. Guan, “Secure consensus of multi-agent systems via impulsive control subject to deception attacks,” _IEEE Transactions on Circuits and Systems II: Express Briefs_ , pp. 1–1, 2022.
  * [9] A.-A. Gao, A.-H. Hu, and Z.-X. Jiang, “Secure consensus of multi-agent systems based on event-triggered impulsive control,” _Journal of Computer Applications_ , p. 0, 2022.
  * [10] S. M. Elkhider, S. El-Ferik, and A.-W. A. Saif, “Denial of service attack of qos-based control of multi-agent systems,” _Applied Sciences_ , vol. 12, no. 9, 2022.
  * [11] J. Wang, Y. Li, Z. Duan, and J. Zeng, “A fully distributed robust secure consensus protocol for linear multi-agent systems,” _IEEE Transactions on Circuits and Systems II: Express Briefs_ , vol. 69, no. 7, pp. 3264–3268, 2022\.
  * [12] M. Sathishkumar and Y.-C. Liu, “Resilient memory event-triggered consensus control for multi-agent systems with aperiodic dos attacks,” _International Journal of Control, Automation and Systems_ , pp. 1–14, 2022\.
  * [13] C. Gao, X. He, H. Dong, H. Liu, and G. Lyu, “A survey on fault-tolerant consensus control of multi-agent systems: trends, methodologies and prospects,” _International Journal of Systems Science_ , vol. 53, no. 13, pp. 2800–2813, 2022.
  * [14] L. Xiao and S. Boyd, “Fast linear iterations for distributed averaging,” _Systems & Control Letters_, vol. 53, no. 1, pp. 65–78, 2004.
  * [15] J. Lunze, _Networked control of multi-agent systems_. Edition MoRa, 2019.
  * [16] D. Zelazo and M. Bürger, “On the robustness of uncertain consensus networks,” _IEEE Transactions on Control of Network Systems_ , vol. 4, no. 2, pp. 170–178, 2017.
  * [17] M. Fabris, G. Michieletto, and A. Cenedese, “A general regularized distributed solution for system state estimation from relative measurements,” _IEEE Control Systems Letters_ , vol. 6, pp. 1580–1585, 2022.
  * [18] A. Chapman, _Semi-Autonomous Networks: Effective Control of Networked Systems through Protocols, Design, and Modeling_. Springer, 2015.