# Beyond-Mean-Field with an Effective Hamiltonian Mapped from an Energy
Density Functional

J. Ljungberg  J. Boström  B. G. Carlsson  A. Idini and J. Rotureau
Mathematical Physics, LTH, Lund University, S-22100 Lund, Sweden
<EMAIL_ADDRESS>
###### Abstract

A method for beyond-mean-field calculations based on an energy density
functional is described. The main idea is to map the energy surface for the
nuclear quadrupole deformation, obtained from an energy density functional at
the mean-field level, into an effective Hamiltonian expressed as a many-body
operator. The advantage of this procedure is that one avoids the problems with
density dependence which can arise in beyond-mean-field methods. The effective
Hamiltonian is then used in a straightforward way in the generator-coordinate-
method with the inclusion of projections onto good particle numbers and
angular momentum. In the end, both spectra and wave functions are obtained. As
an example of the method, calculations for the nucleus 62Zn is performed with
three different parametrizations of the Skyrme functional. The results are
compared with experiment.

##  1 Introduction

A successful route to understanding and describe atomic nuclei is the mean-
field approach. In this view, all the nucleons build up a common potential in
which they move independently. This potential can, for example, be found in a
self-consistent way through the Hartree-Fock (HF) method or, with the
inclusion of pairing correlations, the Hartree-Fock-Bogoliubov (HFB) method.
In those methods, the state of the nucleus is one single Slater-determinant of
single-particle or quasiparticle nature.

The mean-field approach is able to describe bulk properties of nuclei, such as
the binding energy, with high accuracy. Also, it is common that the solutions
break one or more symmetries of the underlying interaction. The most prominent
is the breaking of spherical symmetry into a deformed nucleus. This symmetry
breaking provides a way to gain physical insight into nuclear phenomena. For
example, by interpreting the deformation as a static shape of the nucleus in
its intrinsic frame, one can understand the emergence of rotational bands
which show up in many nuclear spectra.

However, despite the success of the mean-field description, it can not be the
complete description of the system. Due to the breaking of symmetries, it is
impossible for the solutions to be true eigenstates of the nuclear
Hamiltonian. One would expect a true eigenstate, which respects the
symmetries, to be a super position of different shapes with different
orientations. Furthermore, a deformed state can not be labeled with good
quantum numbers for the angular momentum. And that is crucial in order to get
proper theoretical spectra and to calculate transitions between states;
without the use of extra assumptions for going from the intrinsic frame into
the laboratory frame.

The method presented here belongs to the class of beyond-mean-field methods.
The common feature for those methods is that the state of the nucleus can be
expressed as a linear combination of several Slater-determinants; as opposite
to the single one in mean-field methods.

In our model, we want to restore the broken symmetries of particle numbers and
angular momentum with the use of projections. Also, we want to include
correlations associated with large amplitude collective motions; that is:
vibrations and rotations. This is done through the generator-coordinate-method
(GCM) together with the cranking method. As the generator coordinates we use
the $\beta$ and $\gamma$ deformation parameters.

In the search of a theoretical method to describe and predict low energy
properties of atomic nuclei, a common starting point is to introduce a general
two- and three-body interaction. In order to perform practical computations,
the three-body interaction is often converted to a density dependent two-body
interaction; turning the interaction into an energy density functional (EDF)
[1]. A commonly used EDF is based on the Skyrme expansion [2].

Using an EDF in the framework of HFB will reproduce experimental binding
energies with high precisions. Therefore, it can be a good idea to use those
results as input to a beyond-mean-field method.

However, extending the mean-field models based on EDF:s with projections and
GCM has been shown to give rise to formal and technical issues. Firstly, some
effective interactions, which are used for the low-momentum regime of nuclei,
do also contain a non-physical high-momentum part. This part is not probed at
the mean field level; whereas it can be so in GCM when mixing different many-
body states. Those interactions can give rise to ultraviolet divergences [3].

Secondly, how to treat the density in the mixing process is not well defined;
although many approximate ways have been proposed [4, 5, 6, 7].

Therefore it can be important to have the interaction expressed as a many-body
Hamiltonian with no density dependence when aiming for beyond-mean-field
calculations.

##  2 Model

In our model, we use results from an EDF as input to a beyond-mean-field
calculation through the frameworks of GCM and projections. In the end, the
model produces spectra and wave functions. The wave functions can be used to
calculate observables; such as electromagnetic transitions. In the future, we
want to use the wave functions as input to reaction studies by constructing an
optical potential out of them; see Ref. [8] in this volume.

A detailed description of the model can be found in Ref. [9]. Here follows a
summary.

###  2.1 The Effective Hamiltonian

The starting point for this model is to postulate an effective Hamiltonian of
the form

$H_{eff}=\sum_{i}e_{i}a^{\dagger}_{i}a_{i}-\\\
\frac{1}{4}G\sum_{ijkl}P_{ij}P_{kl}a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}-\\\
\frac{1}{4}\chi\sum_{\mu}\sum_{ijkl}\left[\tilde{Q}^{2\mu}_{ik}\tilde{Q}^{2\mu*}_{lj}-\tilde{Q}^{2\mu}_{il}\tilde{Q}^{2\mu*}_{kj}\right]a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}.$
(1)

Here, the first term is taken from a spherical HF calculation based on an EDF.
This calculation constructs a single-particle basis for which
$a^{\dagger}_{i}$ is a creation operator for the spherical orbital $i$ with
energy $e_{i}$. The second term is a pairing term. We use the seniority
pairing interaction with the strength $G$; which is determined from the
uniform spectra method [10]. The third term is a modified quadrupole
interaction adopted from Ref. [11]. At the mean-field level, this interaction
gives rise to a deformed Wood-Saxon potential. The strength of the quadrupole
interaction, $\chi$, is a free parameter. It will be determined through a fit
to constrained EDF calculations.

The form of $H_{eff}$ has been chosen in such a way that it includes the most
important features of atomic nuclei. Namely: spherical single-particle
energies, pairing correlations and quadrupole deformations.

###  2.2 Mapping an EDF into $H_{eff}$

The parameter $\chi$ is selected such that $H_{eff}$ reproduces results
obtained from an EDF. That is: Mapping an EDF into $H_{eff}$. Here we choose
to calculate the binding energy as a function of the quadrupole deformation,
$E(\beta_{2})$, from constrained HF calculations; both with an EDF and with
$H_{eff}$. Then $\chi$ is selected such that the two curves for $E(\beta_{2})$
coincides as much as possible.

###  2.3 GCM and Projections

With $\chi$ fixed, $H_{eff}$ can be used in the framework of GCM in a
straightforward way. In GCM, the eigenenergies and eigenstates are extracted
from the space spanned by a many-body basis. In our construction of this basis
we include cranking. We define a cranked Hamiltonian as
$H_{eff}^{\omega}\equiv H_{eff}-\omega I_{x}$ with some constraint on
$\expectationvalue{I_{x}}$. Here $I_{x}$ is the angular momentum around the
$x$-axis and $\omega$ is the angular frequency which acts like a Lagrange
multiplier. The basis states are obtained by solving the HFB equations for
$H_{eff}^{\omega}$ with different constraints on the $\beta$ and $\gamma$
deformation parameters; generating a many-body basis of HFB vacua spread out
over the $(\beta,\gamma)$-plane. In addition, each HFB vacuum is assigned
different values for $\expectationvalue{I_{x}}$ and for the pairing gaps,
which are used to determine $G$, of both neutrons and protons. Those values
are randomly drawn from a set of a few allowed values for each parameter.

In the framework of GCM, the matrix elements of $H_{eff}$ and the overlap,
$\mathcal{O}$, in the many-body basis $\\{\ket{\phi_{a}}\\}$ are needed. For
the overlap we use the new formula of Ref. [12] in our model.

At this stage, projections onto good particle numbers and angular momentum is
performed. Hence, the matrix elements to be calculated are

$H_{a^{\prime}K^{\prime},aK}=\matrixelement{\phi_{a^{\prime}}}{H_{eff}P^{I}_{K^{\prime}K}P^{N}P^{Z}}{\phi_{a}},\\\
\hskip
28.45274pt\mathcal{O}_{a^{\prime}K^{\prime},aK}=\matrixelement{\phi_{a^{\prime}}}{P^{I}_{K^{\prime}K}P^{N}P^{Z}}{\phi_{a}}.$
(2)

Here $P^{I}_{K^{\prime}K}$ is the projection operator for total angular
momentum $I$ with projection $K(K^{\prime})$ on the symmetry-axis for state
$\ket{\phi_{a}}(\ket{\phi_{a^{\prime}}})$; and $P^{N/Z}$ is the projection
operator for particle number of neutrons/protons. See Ref. [13] for the
explicit expressions of those projection operators and how to discretize them
in numerical computations.

In the end, the spectrum and wave functions, for a given angular momentum, are
obtained from the Hill-Wheeler equation

$\sum_{aK}H_{a^{\prime}K^{\prime},aK}h_{aK}^{n}=E_{n}\sum_{aK}\mathcal{O}_{a^{\prime}K^{\prime},aK}h_{aK}^{n}$
(3)

where $E_{n}$ is the energy and $\\{h_{aK}^{n}\\}$ are the expansion
coefficients for for the $n$:th state in the basis of projected HFB states.
Hence, the final state can be written

$\ket{IM,n}=\sum_{aK}h_{aK}^{n}P^{I}_{MK}P^{N}P^{Z}\ket{\phi_{a}}$ (4)

where $M$ is the projection of the angular momentum on the laboratory
$z$-axis.

In practice, the Hill-Wheeler equation is solved by first diagonalizing the
overlap matrix to obtain the so called natural basis. Then, in that basis, the
Hamiltonian is diagonalized. However, this can be numerically unstable because
of the overcompleteness of the many-body basis. Therefore we diagonalize the
Hamiltonian repeatedly. First in a small basis of just a few natural states
that have the largest norms. Then adding one state for each diagonalization;
sorted by decreasing norm. A converged result is indicated by the formation of
a plateau; where the addition of states no longer change the energy.

##  3 Results

As a demonstration of our method, and to give examples of the outputs from the
model, calculations have been performed for the nucleus ${}^{62}\text{Zn}$.
For this nucleus it exists a rich amount of data for which the theoretical
results can be compared to. The calculations has been done with three
different parametrizations of the Skyrme functional: UNE1, SLy4 and SKM*.

###  3.1 Mapping

For the HF calculations, 12 major shells in the spherical oscillator basis has
been used. The EDF results are from the code HFBTHO [14]; whereas the results
from $H_{eff}$ are obtained from an updated version of the code HOSPHE [15].
The curve $E(\beta_{2})$ has been fitted between $-0.15\leq\beta_{2}\leq
0.30$. In Fig. 1 the fits for the three parametrizations are presented.

Figure 1: Fits of the parameter $\chi$ for the three parametrizations of the
Skyrme functional given at the top of each plot. The blue curves are the HF
energies obtained from the EDF; whereas the red curves are HF energies
obtained from $H_{eff}$. The energy is the binding energy; which can be
compared to the experimental value: $E=-538.1$ MeV [16]. Since pairing
correlations are not included, we expect the HF energies to be higher than the
experimental one.

Even though the three parametrizations generate quite different $E(\beta_{2})$
curves, $H_{eff}$ can capture them all in a qualitative way.

###  3.2 The Many-Body Basis

For UNE1, the $(\beta,\gamma)$-plane was sampled with 250 points within
$-30^{\circ}\leq\gamma\leq 150^{\circ}$ and $\beta\leq 0.36$. Only states
below the excitation energy $E^{*}=16$ MeV were kept. This resulted in basis
size of 209 states.

For SLy4 and SKM*, the $(\beta,\gamma)$-plane was sampled with 400 points
within $-30^{\circ}\leq\gamma\leq 150^{\circ}$ and $\beta\leq$(0.60 / 0.48)
for (SLy4 / SKM*). Furthermore, a quadratic cut in the excitation energy was
chosen; only states below $E^{*}=[8+0.07\times I_{x}(I_{x}+1)]$ MeV were kept.
The number 0.07 was selected to fit with the experimental yrast band. This
resulted in a basis size of (151 / 145) states for (SLy4 / SKM*).

For the discretization of the projection onto good particle numbers, 10 points
was used. For the angular momentum, 28 points for each Euler angle was used
for UNE1; whereas 32 points was used for SLy4 and SKM*. All spins up to
$10\hbar$ was projected out.

###  3.3 Spectra

The ground state of 62Zn can be viewed as a core plus six valence nucleons,
four neutrons and two protons, in the $fp$-group above the 28 gap. At low
excitation energies, the valence nucleons are limited to this group. The
maximal spin that can be built in this configuration, by aligning all valence
nucleons, is $10\hbar$. Hence, one expects the yrast band to terminate at
$10\hbar$. At higher energies two neutrons, to ensure positive parity, can be
excited into the $g_{9/2}$ subshell. This configuration terminates at
$16\hbar$. In this work we focus on the low energy regime up to $I=10\hbar$.

The spectra is calculated in the way described in section 2.3. As a
visualization of the procedure, Fig. 2 shows the energy obtained, for the
yrast states of even angular momentum, with different sizes of the natural
basis. This is from the calculation built upon UNE1.

Figure 2: The energy obtained from diagonalizing the Hamiltonian in a basis
with different numbers of natural states. The curves are for the yrast states
of even spins. The crosses show where the solutions have been selected. This
is from the calculation built upon UNE1.

The region where the solution becomes unstable is distinct for all spins above
zero. The solutions are taken just before those regions; shown by the crosses.

Figure 3: The spectrum for 62Zn from experiment [16] and from our model built
upon the parametrizations UNE1, SLy4 and SKM*. For the experimental data, only
states with positive parity and $I\leq 10\hbar$ are shown. The blue lines
represent the bands which are found in experiment and their corresponding
positions in the theoretical spectra.

The final results for the spectra is shown in Fig. 3. Also experimental
values, taken from [16], are shown for comparison. First consider the UNE1
parametrization. The overall structure close to the yrast band is in good
agreement with experiments; even though the energies increases too fast with
angular momentum. It is interesting to compare the band structure. In
experiment, there is a clear yrast band up to $I=10\hbar$ and two bands on top
of each other just above the yrast one; ranging from $I=2\hbar$ to $I=7\hbar$.
This structure of states is reproduced by the calculation. Compared to our
previous work [9], we have included odd spin states in the model. The
agreement with experiment for those states indicates that our method holds
also for the odd spin case.

For the other parametrizations, SLy4 gives very similar spectra as UNE1.
Hence, for those parametrizations, the method is not too sensitive to the
underlying EDF it is founded on.

In contrast to that, SKM* does not succeed as well as the other
parametrizations to reproduce the structure above the yrast band. And the
increase in energy with angular momentum is even higher; in accordance with
the different energy landscape in Fig. 1.

##  4 Discussion

A beyond-mean-field method has been described. This method uses the results
from an EDF to constrain an effective Hamiltonian. This Hamiltonian is then
used in the framework of GCM, in order to include correlations associated with
vibrations and rotations, and projections are used to restore symmetries. The
form of the effective Hamiltonian allows for many-body calculations in a
straightforward way.

The results shown indicate, despite the simplicity of the interaction, that
the effective Hamiltonian capture enough physics in order to reproduce
experimental data. And also, because of this simplicity, the model can be a
valuable tool for the description of heavy and superheavy nuclei [17].

## References

## References

  * [1] Bender M, Heenen P H and Reinhard P G 2003 Rev. Mod. Phys. 75(1) 121–180
  * [2] Skyrme T 1958 Nuclear Physics 9 615–634 ISSN 0029-5582
  * [3] Carlsson B G, Toivanen J and von Barth U 2013 Phys. Rev. C 87(5) 054303
  * [4] Satuła W and Dobaczewski J 2014 Phys. Rev. C 90
  * [5] Duguet T and Bonche P 2003 Phys. Rev. C 67(5) 054308
  * [6] Waroquier M, Ryckebusch J, Moreau J, Heyde K, Blasi N, van der Werf S and Wenes G 1987 Phys. Rep. 148 249–306
  * [7] Bender M, Duguet T and Lacroix D 2009 Phys. Rev. C 79(4) 044319
  * [8] Idini A, Rotureau J, Boström J, Ljungberg J and Carlsson B G In this volume
  * [9] Ljungberg J, Carlsson B G, Rotureau J, Idini A and Ragnarsson I 2022 Phys. Rev. C 106(1) 014314
  * [10] Nilsson S and Ragnarsson I 1995 Shapes and Shells in Nuclear Structure (Cambridge: Cambridge University Press)
  * [11] Kumar K and Sørensen B 1970 Nucl. Phys. A 146 1–14
  * [12] Carlsson B G and Rotureau J 2021 Phys. Rev. Lett. 126(17) 172501
  * [13] Bally B and Bender M 2021 Phys. Rev. C 103(2) 024315
  * [14] Perez R N, Schunck N, Lasseri R D, Zhang C and Sarich J 2017 Comput. Phys. Commun. 220 363–375
  * [15] Carlsson B, Dobaczewski J, Toivanen J and Veselý P 2010 Comput. Phys. Commun. 181 1641–1657
  * [16] 2022 Evaluated nuclear structure data file
  * [17] Såmark-Roth A, Cox D M, Rudolph D, Sarmiento L G, Carlsson B G, Egido J L, Golubev P, Heery J, Yakushev A, Åberg S et al. 2021 Phys. Rev. Lett. 126(3) 032503


# Attack-Resilient Distributed Algorithms for Exponential Nash Equilibrium
Seeking

Zhi Feng and Guoqiang Hu  This work was supported by Singapore Ministry of
Education Academic Research Fund Tier 1 RG180/17(2017-T1-002-158). Z. Feng and
G. Hu are with the School of Electrical and Electronic Engineering, Nanyang
Technological University, Singapore 639798 (E-mail<EMAIL_ADDRESS>gqhu@ntu.edu.sg).

###### Abstract

This paper investigates a resilient distributed Nash equilibrium (NE) seeking
problem on a directed communication network subject to malicious cyber-
attacks. The considered attacks, named as Denial-of-Service (DoS) attacks, are
allowed to occur aperiodically, which refers to interruptions of communication
channels carried out by intelligent adversaries. In such an insecure network
environment, the existence of cyber-attacks may result in undesirable
performance degradations or even the failures of distributed algorithm to seek
the NE of noncooperative games. Hence, the aforementioned setting can improve
the practical relevance of the problem to be addressed and meanwhile, it poses
some technical challenges to the distributed algorithm design and exponential
convergence analysis. In contrast to the existing distributed NE seeking
results over a prefect communication network, an attack-resilient distributed
algorithm is presented such that the NE can be exactly reached with an
exponential convergence rate in the presence of DoS attacks. Inspired by the
previous works in [23, 21, 22, 24, 25, 26], an explicit analysis of the attack
frequency and duration is investigated to enable exponential NE seeking with
resilience against attacks. Examples and numerical simulation results are
given to show the effectiveness of the proposed design.

###### Index Terms:

Distributed algorithm, NE seeking, Directed graph, Exponential convergence,
Cyber-attack, Resilience.

##  I Introduction

Recently, distributed Nash equilibrium (NE) seeking of non-cooperative games
has been attracting increasingly attention due to its broad applications in
multi-robot systems [1], mobile sensor networks [3], smart grids [5], and so
on. In contrast to early works (e.g., [2, 4, 6]) with a complete information
setting, players in distributed NE seeking have limited local information,
i.e., each player needs to make the decision based on the local or relative
information, e.g., information from its neighbors, to optimize its own
objective function. Thus, the players involved in the game are required to
communicate with each other over the network to estimate other players’
actions.

Literature review: various distributed continuous- and discrete-time
algorithms are developed in [8, 9, 11, 10, 12, 13, 14, 15, 16, 17, 18] to
solve distributed NE seeking problems in games, in which each player cannot
observe all other players’ actions, but can exchange information between
neighbors over an undirected graph or weight-balanced digraph. Gradient-based
NE seeking algorithms are popular techniques to find the NE of games with
differentiable objective functions where each player modifies its current
action based on the gradient with respect to its own action. The distributed
NE seeking issue with continuous-time agent dynamics is studied in [8, 9, 11,
10, 12, 13, 14]. In particular, distributed NE seeking strategies are proposed
in [8] and [9] by combining the leader-follower based consensus algorithms and
gradient-play strategies over an undirected and connected graph. The authors
in [11] exploit some incremental passivity properties of pseudo-gradients to
illustrate that the estimates of the proposed augmented gradient dynamics
converge to the NE exponentially under the graph coupling condition. To remove
this condition, the two time-scale singular perturbation analysis is further
developed. Distributed NE seeking of aggregative games has been studied in
[12, 13, 14], where each player’ objective function relies on its own action
and an aggregate of all players’ actions. On the other hand, discrete-time
consensus algorithms are presented in [15] to search for the NE of congestion
games, where each player has the linear cost function. In [16], the gossip-
based algorithm is developed to seek the NE. An alternating direction method
of multipliers’ algorithm with the constant step-size is proposed in [17] by
exploiting the convex and smooth properties of pseudogradients. Recently, [18]
adopts consensus-based gradient-free NE seeking algorithms for limited cost
function knowledge. Distributed algorithms via diminishing and constant step-
sizes are studied, where the former ensures an almost sure NE convergence,
while the latter provides convergence to the neighborhood of the NE.

One observation in distributed NE seeking problems is that all players require
to find the NE via information exchange between neighbors through a
communication network. Unfortunately, due to malicious cyber-attacks such as
DoS attacks, deception attacks (false data injection, replay attack),
disclosure attacks [20, 27, 28, 29, 30, 21, 23, 22, 24, 25, 26], the secure
network environment is hardly guaranteed in practice. The network security
plays a fundamental yet critical role in successful information transmission.
The malicious attacks interrupt, incorrect, or tamper transmitting information
so that efficiency of distributed NE seeking algorithms is degraded
significantly, and it might even lead to the failure of seeking the NE under
malicious attacks. In light of wide applications of distributed algorithms in
cyber-physical systems (safety-critical), and inspired by studies of security
issues in many existing works (e.g., see [20, 27, 28, 29, 30, 21, 23, 22, 24,
25, 26]), it is highly desirable to determine how resilient distributed NE
seeking algorithms will be against malicious attacks based on the fact that
distributed algorithms are easily disrupted by malicious behaviors. Thus, the
main objective of this work is to address attack-resilient problems of
distributed NE seeking so as to provide certain safety and resilience
performances against attacks. So far, few efforts are made on resilient
distributed NE seeking.

This paper focuses on the attack-resilient research of distributed NE seeking
of a non-cooperative game. In particular, an attack-resilient distributed
algorithm is proposed to solve the NE seeking problem in the multi-player non-
cooperative game over a directed communication network under malicious DoS
attacks. The major contributions of this paper are as follows. Firstly, to the
best of our knowledge, this paper is the first work to solve this issue. The
proposed distributed NE seeking algorithm is capable of exactly seeking the NE
in the presence of DoS attacks. The exponential convergence of the proposed
algorithm is rigorously guaranteed, provided that the frequency and duration
of attacks satisfy certain bounded conditions. Moreover, compared with
existing NE seeking works requiring ideal communication in [8, 9, 11, 10, 12,
13, 14, 15, 16, 17, 18], we develop attack-resilient distributed NE seeking
algorithms to search for the NE in adversarial network environments. The
algorithms employ consensus-based pesudo-gradient strategies with a hybrid
system method to constrain DoS attacks. The explicit analysis is provided
according to the Lyapunov stability. In addition, different from the
distributed convex optimization works in [31, 32, 33] that consider faults on
nodes and require the removal of their states to be prior known, this paper
investigates the more practical DoS attack on communication network. The
attacks are time-sequence based and allowed to occur aperiodically. Another
contribution of this work is that unlike results in [8, 9, 11, 10, 12, 13, 14,
15, 16, 17, 18] and [31, 32, 33] over an undirected graph or weight-balanced
digraph, the strongly connected directed graph is allowed here, and in the
presence of malicious attacks, this directed graph can be disrupted or totally
paralyzed.

The paper is organized as follows. Section II gave mathematical preliminaries.
In Section III, the non-cooperative game and DoS attack model are described,
and the main objective is presented. The attack-resilient distributed
algorithm is proposed in Section IV, where the NE seeking results with an
exponential convergence analysis are provided. Examples and numerical
simulation results are given in Section V, followed by the conclusion in
Section VI.

##  II Preliminaries

###  II-A Notation

Denote $\mathbb{R}$, $\mathbb{R}^{n}$, and $\mathbb{R}^{n\times m}$ as the
sets of the real numbers, real $n$-dimensional vectors and real $n\times m$
matrices, respectively. Denote $\mathbb{R}_{\geq 0}$ as the set of nonnegative
numbers, while $\mathbb{N^{+}}$ represents the set of positive integers. Let
$0_{n}$ ($1_{n}$) be the $n\times 1$ vector with all zeros (ones) and $I_{n}$
be the identity matrix. Let col$(x_{1},...,x_{n})$ and
diag$\\{a_{1},...,a_{n}\\}$ be a column vector with entries $x_{i}$ and a
diagonal matrix with entries $a_{i}$, $i=1,2,\cdots,n$, respectively. The
symbols $\otimes$ and $\left\|\cdot\right\|$ represent the Kronecker product
and the Euclidean norm, respectively. Given a real symmetric matrix $M$, let
$M>0$ ($M\geq 0$) denote that $M$ is positive (or positive semi-definite), and
$\lambda_{\min}(M)$, $\lambda_{\max}(M)$ are its minimum and maximum
eigenvalues, respectively. For a function $f$, it is said to be
$\mathcal{C}^{m}$ if it is $m$th continuously differentiable. For two sets
${X}$ and ${Y}$, ${Y}\backslash{X}$ denotes the set of elements belonging to
$Y,$ but not to $X.$

###  II-B Convex analysis

A function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is convex if
$f(ax+(1-a)y)\leq af(x)+(1-a)f(y)$ for any scalar $a\in[0,1]$ and vectors
$x,y\in\mathbb{R}^{n}$. $f$ is locally Lipschitz on $\mathbb{R}^{n}$ if it is
locally Lipschitz at $x$ for $\forall x\in\mathbb{R}^{n}$. If $f$ is a
differentiable function, $\triangledown f$ denotes the gradient of $f$. A
vector-valued function (or mapping)
$F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is said to be
$\iota_{F}$-Lipschitz continuous if, for any $x,y\in\mathbb{R}^{n}$,
$\|F(x)-F(y)\|\leq\iota_{F}\|x-y\|$. Function
$F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is (strictly) monotone if, for any
$x,y\in\mathbb{R}^{n}$, $(x-y)^{T}(F(x)-F(y))(>)\geq 0$. Further, $F$ is a
$\varepsilon$-strongly monotone, if for any scalar $\varepsilon>0$, and
$x,y\in\mathbb{R}^{n}$, $(x-y)^{T}(F(x)-F(y))>\varepsilon\|x-y\|^{2}$.

###  II-C Graph Theory

Let $\mathcal{\bar{G}}$ $=$ $\left\\{\mathcal{V},\mathcal{\bar{E}}\right\\}$
denote a directed graph, where $\mathcal{V}$ $\in$ $\left\\{1,...,N\right\\}$
is a set of nodes and $\mathcal{\bar{E}}$ $\subseteq$ $\mathcal{V\times V}$ is
a set of edges. An edge $(i,j)\in\mathcal{\bar{E}}$ denotes that $i$th agent
receives the information from $j$th agent, but not vice versa.
$\mathcal{\bar{N}}_{i}$ $=$
$\left\\{j\in\mathcal{V\mid}(j,i)\in\mathcal{\bar{E}}\right\\}$ is a
neighborhood set of the agent $i$. A directed graph is strongly connected if
there exists a directed path connecting every pair of nodes. The matrix
$\bar{A}=\left[\bar{a}_{ij}\right]$ denotes the adjacency matrix of
$\mathcal{\bar{G}}$, where $\bar{a}_{ij}>0$ if and only if
$(j,i)\in\mathcal{\bar{E}}$, else $\bar{a}_{ij}=0$. A matrix
$\mathcal{\bar{L}}=[\bar{l}_{ij}]\in\mathbb{R}^{N\times N}$ is called the
Laplacian matrix of $\mathcal{\bar{G}}$, where
$\bar{l}_{ii}=\sum^{N}_{j=1}\bar{a}_{ij}$ and $\bar{l}_{ij}=-\bar{a}_{ij}$,
$i\neq j$. Denote $\mathcal{\bar{L}}\triangleq\bar{D}-\bar{A}$, where
$\bar{D}=\left[\bar{d}_{ii}\right]$ is the diagonal matrix with
$\bar{d}_{ii}=\sum\nolimits_{j=1}^{N}\bar{a}_{ij}$. Similarly, let
$\mathcal{G}$ denote a weighted digraph, where $a_{ij}=\omega_{i}\bar{a}_{ij}$
for a scalar $\omega_{i}>0$. With the same definition of $\mathcal{\bar{L}}$,
the Laplacian matrix of this new digraph becomes
$\mathcal{L}=\mathcal{W}\mathcal{\bar{L}}$, where
$\mathcal{W}=\text{diag}\\{\omega_{1},\cdots,\omega_{N}\\}$.

###### Assumption 1

The directed graph $\mathcal{\bar{G}}$ is strongly connected.

###### Lemma 1

[37] By Assumption 1, the Laplacian matrix $\mathcal{\bar{L}}$ of
$\mathcal{\bar{G}}$ has the following properties.

  1. 1.

$\mathcal{\bar{L}}$ has a simple zero eigenvalue associated with right
eigenvector $1_{N}$, and nonzero eigenvalues have positive real parts;

  2. 2.

$\omega=\text{col}(\omega_{1},\cdots,\omega_{N})$ with $\omega^{T}1_{N}=1$ is
a left eigenvector of $\mathcal{\bar{L}}$ associated with the zero eigenvalue.
Then, $\omega^{T}\mathcal{\bar{L}}=0^{T}_{N}$ and $\min_{\zeta^{T}x=0,\ x\neq
0}\frac{x^{T}\mathcal{\hat{L}}x}{x^{T}x}>\frac{\lambda_{2}(\mathcal{\hat{L}})}{N}$,
where $\mathcal{\hat{L}}=\frac{1}{2}(\mathcal{L}+\mathcal{L}^{T})$ and $\zeta$
is any vector with positive entries. Moreover, $\omega=1_{N}$ if and only if
$\mathcal{\bar{G}}$ is strongly connected and weight-balanced.

##  III Problem Formulation

###  III-A Non-cooperative Game over Networks

In this paper, we consider a multi-agent network consisting of $N$ players,
which form a N-player non-cooperative game defined as follows. For each agent
$i\in\mathcal{V}$, the $i$th player aims to minimize its cost function
$J_{i}(x_{i},x_{-i}):\mathbb{R}^{n_{i}}\rightarrow\mathbb{R}$ by choosing its
strategy $x_{i}\in\mathbb{R}^{n_{i}}$, and
$x_{-i}=\text{col}(x_{1},\cdots,x_{i-1},\cdots,x_{N})$ is the strategy profile
of the whole strategy profile except for player $i$. Let $x=\\\
(x_{i},x_{-i})$ represent all players’ action profile. Alternatively, let
$x\\\ =\text{col}(x_{1},\cdots,x_{N})\in\mathbb{R}^{n}$,
$n=\sum_{i\in\mathcal{V}}n_{i}$.

The concept of Nash equilibrium is given below.

###### Definition 1

(Nash equilibrium) A strategy profile $x^{*}=(x^{*}_{i},\\\
x^{*}_{-i})\in\mathbb{R}^{n}$ is said to be an Nash equilibrium of the game if

$J_{i}(x^{*}_{i},x^{*}_{-i})\leq J_{i}(x_{i},x^{*}_{-i}),\ \text{for}\ \forall
x_{i}\in\mathbb{R}^{n_{i}},\ i\in\mathcal{V}.$ (1)

Condition (1) means that all players simultaneously take their own best
(feasible) responses at the NE $x^{*}$, where no player can unilaterally
decrease its cost by changing its strategy.

###### Assumption 2

For each player $i$, $J_{i}(x_{i},x_{-i})$ is $\mathcal{C}^{2}$, strictly
convex, and radially unbounded in $x_{i}$ for each $x_{-i}$.

Under Assumption 2, it follows from [35] that an NE $x^{*}$ exists, and
satisfies $\triangledown_{i}J_{i}(x^{*}_{i},x^{*}_{-i})=0_{n_{i}}$, and
$\triangledown_{i}J_{i}(x_{i},x_{-i})=\partial J_{i}(x_{i},\\\
x_{-i})/\partial x_{i}\in\mathbb{R}^{n_{i}}$ represents the partial gradient
of player $i$’s cost with respect to its own action $x_{i}$. We define

$F(x)\triangleq\text{col}(\triangledown_{1}J_{1}(x_{1},x_{-1}),\cdots,\triangledown_{N}J_{N}(x_{N}\\\
,x_{-N})),$ (2)

where $F(x)\in\mathbb{R}^{n}$ denotes the pseudo-gradient (the stacked vector
of all players’ partial gradient). Thus, we have $F(x^{*})=0_{n}$.

###### Assumption 3

The pesudogradient $F$ is $\varepsilon$-strongly monotone and
$\iota_{F}$-Lipschitz continuous for certain constants
$\varepsilon,\iota_{F}>0$.

###### Remark 1

Assumptions 2 and 3 were widely used in existing works (e.g., [11, 10, 12, 13,
14, 15, 16, 17, 18]) to guarantee the unique NE $x^{*}$.

In a game with perfect information (a complete communication graph) requiring
the global knowledge on all players’ actions, a gradient-play algorithm
$\dot{x}_{i}=-\triangledown_{i}J_{i}(x_{i},x_{-i})$ can be used to seek the NE
asymptotically (Lemma 2, [11]). However, this design is impractical if the
communication network is not complete. Then, a distributed algorithm is
desirable to broadcast the local information among players. Unfortunately,
when distributed computations are required and executed over insure
communication networks under cyber-attacks, the exact NE may not be found.

###  III-B Malicious DoS Attack Model

As studied in the preliminary results [21, 23, 22, 24, 25, 26], the well-known
DoS attack refers to a class of malicious cyber-attacks where an attacker aims
to corrupt and interrupt certain or all components of communication channels
between the players. Without loss of generality, we suppose that each
adversary follows an independent attacking strategy with a varying active
period. Due to the limited energy, then it terminates attacking activities and
shifts to a sleep period to supply its energy for next attacks.

To be specific, suppose that there exists an $m\in\mathbb{N}$ and denote
$\\{a^{ij}_{m}\\}_{m\in\mathbb{N}}$ as an attack sequence when a DoS attack is
lunched at channel $(i,j)$. Let the duration of this attack be
$\tau^{ij}_{m}\geq 0$. Then, the $m$-th attack strategy is generated with
$\mathcal{A}^{(i,j)}_{m}={a^{ij}_{m}}\cup[a^{ij}_{m},a^{ij}_{m}+\tau^{ij}_{m})$
where $a^{ij}_{m+1}>a^{ij}_{m}+\tau^{ij}_{m}$ for $m\in\mathbb{N^{+}}.$ Given
$t\geq t_{0}$, the sets of time instants where communication between $(i,j)$
is denied, are described in the following form [24, 25, 26]:

$\Xi^{(i,j)}_{a}(t_{0},t)=(\cup^{\infty}_{m=1}\
\mathcal{A}^{(i,j)}_{m})\cap[t_{0},t],\ m\in\mathbb{N^{+}},$ (3)

which implies that the sets where communication between $(i,j)$ is allowed,
are: $\Xi^{(i,j)}_{s}(t_{0},t)=[t_{0},t]\setminus\Xi^{(i,j)}_{a}(t_{0},t)$.
Then, $|\Xi_{a}(t_{0},t)|\\\
=|{\cup}_{(i,j)\in\mathcal{E}}\Xi^{(i,j)}_{a}(t_{0},t)|$ is the total length
of attacker being active over $[t_{0},t]$, while
$|\Xi_{s}(t_{0},t)|=[t_{0},t]\setminus|\Xi_{a}(t_{0},t)|$ denotes the total
length of attacker being sleeping over $[t_{0},t]$.

###### Definition 2

(Attack Frequency) For any $T_{2}>T_{1}\geq t_{0}$, let $N_{a}(T_{1},T_{2})$
denote the total number of DoS attacks over $[T_{1},T_{2})$. Then,
$F_{a}(T_{1},T_{2})=\frac{N_{a}(T_{1},T_{2})}{T_{2}-T_{1}}$ denotes the attack
frequency over $[T_{1},T_{2})$ for $\forall T_{2}>T_{1}\geq t_{0}$, where
there exists scalars $N_{0},T_{f}>0$ so that $N_{a}(T_{1},T_{2})\leq
N_{0}+(T_{2}-T_{1})/T_{f}.$

###### Definition 3

(Attack Duration) For any $T_{2}>T_{1}\geq t_{0},$ denote
$|\Xi_{a}(T_{1},T_{2})|$ as the total time interval under attacks over
$[T_{1},T_{2}).$ The attack duration over $[T_{1},T_{2})$ is defined as: there
exist scalars $T_{0}\geq 0$, $T_{a}>1$ so that $|\Xi_{a}(T_{1},T_{2})|\leq
T_{0}+(T_{2}-T_{1})/T_{a}.$

###### Remark 2

As investigated in the pioneer works [21, 23, 22, 24, 25, 26] for multi-agent
systems under attacks, Definitions 2 and 3 are firstly introduced in [21] to
specify attack signals in terms of frequency and time-ratio constraints. The
$m$-th DoS attacks occurring at $a^{ij}_{m}$ between communication channels
$(i,j)$ with $\tau^{ij}_{m}$ are allowed to occur aperiodically and can
interrupt any communication channels synchronously or asynchronously. Fig. 1
provides an example for more details. In Definition 2, $1/T_{f}$ provides an
upper bound on the average DoS frequency, while $1/T_{a}$ in Definition 3
provides an upper bound on the average DoS duration. It requires attacks to
neither occur at an infinitely fast rate or be always activated.

Figure 1: An illustration of DoS attacks on communication channels $(i,j)$
occurring at $2$s, $9$s, $a^{ij}_{m}$s with their durations being $3$s, $1$s,
$\tau^{ij}_{m}$s, respectively. Thus, we obtain
$\Xi^{(i,j)}_{a}(3,10)=[3,5)\cup[9,10)$ for just an example.

###  III-C Main Objective

This work aims to study an attack-resilient issue of NE seeking under an
insecure communication network as follows.

###### Problem 1

(Attack-Resilient Distributed NE Seeking)
Consider a non-cooperative game consisting of $N$ players communicating over
an insecure communication network induced by DoS attacks. Design an attack-
resilient distributed NE algorithm so that all players can exactly reach the
NE $x^{*}$ with an exponential convergence rate and a resilient feature
against attacks.

$\hskip 10.00002pt\begin{split}&\text{minimize}\ J_{i}(x_{i}(t),x_{-i}(t)),\
x_{i}(t)\in\mathbb{R}^{n_{i}},\ i=1,\cdots,N,\\\ &\text{subject to:}\
\dot{x}_{i}(t)=u_{i}(t),\
t\in\Xi^{(i,j)}_{s}(t_{0},t)\cup\Xi^{(i,j)}_{a}(t_{0},t).\end{split}$ (4)

###### Remark 3

In contrast to existing works in [8, 9, 11, 10, 12, 13, 14, 15, 16, 17, 18],
solving Problem 1 is much more challenging at least from the following
aspects: (1) Player communication network: the games involved in an insecure
communication network may lead to the interruption of communication
transmission caused by DoS attacks, which makes existing NE seeking algorithms
[8, 9, 11, 10, 12, 13, 14, 15, 16, 17, 18] inapplicable. (2) Assumption: in
the absence of attacks, the graph is directed rather than being an undirected
graph or weight-balanced digraph in [8, 9, 11, 10, 12, 13, 14, 15, 16, 17,
18]. Further, this digraph under DoS attacks can be disconnected or totally
paralyzed, which brings nontrivial convergence analysis. (3) Design
requirement: propose an attack-resilient distributed NE seeking scheme with
the exponential convergence and resilience features against attacks. Due to
aforementioned challenges, existing NE seeking algorithms cannot be directly
applied.

###  III-D Motivating Example

In order to illustrate the influence of DoS attacks on distributed NE seeking,
the following motivation example is provided. We consider a classic example in
economy, namely the Nash-Cournot game (e.g., [38, 11]). This game includes
some firms involved in the production of a homogeneous commodity, where the
quantity produced by firm $i\in\mathcal{V}$ is denoted by $x_{i}$, and the
overall cost function of each firm $i$ is described by

$J_{i}(x_{i},x_{-i})=g_{i}(x_{i})-x_{i}f(x),\ i\in\mathcal{V},\\\ $ (5)

where $g_{i}(x_{i})=a_{i}+b_{i}(x_{i}-c_{i})+d_{i}x^{2}_{i}$ is the production
cost with $a_{i},b_{i},c_{i},d_{i}$ describing the characteristics of firm
$i$, and $f(x)=f_{0}-f_{1}\sum_{j=1}^{N}x_{j}$ is the commodity price with
constants $f_{0},f_{1}$. As studied in [11], these parameters are chosen as
$a_{i}=c_{i}=d_{i}=0$, $b_{i}=10+4(i-1)$, $f_{0}=720$, and $f_{1}=1$. By
certain calculations, the NE is $x^{*}=\text{col}(110,106,102,98,94,90)$.

Next, we consider the following different cases:

i) in the absence of attacks, the existing distributed NE seeking algorithm in
[11] is performed, and Fig. 2 shows that all players’ strategies can reach
consensus and converge to the NE.

---
Figure 2: Simulated results generated by the existing algorithm [11] in the
absence of DoS attacks. Left: plot of players’ strategies and the NE
$x^{*}_{i}$; Right: plot of an undirected graph in [11] and its associated
Laplacian matrix.

ii) in the presence of attacks, this algorithm is further performed under
graphs in Fig. 3(a) and the simulation result is shown in Fig. 3(b). As we
see, all players’ strategies neither reach consensus or converge to the NE. In
contrast, the simulation result generated by the proposed attack-resilient
algorithm is shown in 3(c), where all players’ strategies reach consensus and
converge to the NE.

---

(a) Various graphs under DoS attacks on communication channels

(b) Algorithm in [11]

(c) Proposed attack-resilient algorithm

Figure 3: Simulated results generated by different algorithms under DoS
attacks.

##  IV Attack-Resilient Distributed NE Seeking

In distributed NE seeking games, each player $i$ has no access to the full
information of all players’ strategies. Then, each agent $i$ shall estimate
all other players’ strategies. Inspired by [11], let each player combine its
gradient-play dynamics with an auxiliary dynamics, i.e., implement the
following dynamics:

$\left\\{\begin{array}[]{c}\dot{x}^{i}_{i}=u^{i}_{i},\
u^{i}_{i}=-\triangledown_{i}J_{i}(x^{i}_{i},\textbf{x}^{i}_{-i})+e^{i}_{i},\
i\in\mathcal{V},\\\ \dot{x}^{i}_{j}=u^{i}_{j},\ u^{i}_{j}=e^{i}_{j},\ \forall
j\in\mathcal{V},\ j\neq i,\end{array}\right.$ (6)

where player $i$ maintains an estimate vector
${\textbf{x}}^{i}=\text{col}(x^{i}_{1},\cdots,x^{i}_{i},\\\ \cdots,x^{i}_{N})$
in which $x^{i}_{j}$ is player $i$’s estimate of player $j$’s action,
$x^{i}_{i}=x_{i}$ is the player $i$’s actual action, $\textbf{x}^{i}_{-i}$ is
the player $i$’ estimate vector without its own action, $u^{i}_{i}=u_{i}$ is
the player $i$’s actual input, $u^{i}_{j}$ is the other players’ input, and
$e^{i}_{i},e^{i}_{j}$ are to be developed. In (6), each player $i$ updates
$x^{i}_{i}$ to reduce its own cost function and updates $x^{i}_{j}$ to reach
consensus with the other players. In addition, each player $i$ relies on its
local estimated action $\textbf{x}^{i}_{-i}$.

For each player $i$, (6) can be rewritten in a compact form

$\dot{\textbf{x}}^{i}=-\mathcal{R}^{T}_{i}\triangledown_{i}J_{i}(\textbf{x}^{i})+\textbf{e}^{i},\
i\in\mathcal{V},$ (7)

where
$\textbf{e}^{i}=\text{col}(e^{i}_{1},\cdots,e^{i}_{i},\cdots,e^{i}_{N})\in\mathbb{R}^{n}$
is a relative estimated error to be designed, and
$\mathcal{R}_{i}\in\mathbb{R}^{n_{i}\times n}$ used to align the gradient to
the action component, is a matrix given by

$\mathcal{R}_{i}=\left[0_{n_{i}\times n_{1}}\cdots 0_{n_{i}\times n_{i-1}}\
I_{n_{i}\times n_{i}}\ 0_{n_{i}\times n_{i+1}}\cdots 0_{n_{i}\times
n_{N}}\right].$ (8)

###  IV-A Attack-Resilient Distributed NE Seeking Algorithm Design

From the DoS attack model in the subsection III-B, we consider sequence-based
attacks where the $m$-th attack is lunched over a communication channel
$(i,j)$. Without loss of generality, suppose that there exists an infinite
sequence $k=0,1,2,\cdots$ for intervals $[t_{2k},t_{2(k+1)})$ such that in the
absence of DoS attacks, each player $i$ updates its input $\textbf{e}^{i}$
based on an original communication network during $[t_{2k},t_{2k+1})$, while
for each communication channel $(i,j)$ under DoS attacks during
$[t_{2k+1},t_{2k+2})$, the attackers interrupt its information transmission to
make original connected communication network disrupted or totally paralyzed.

To analyze the influence of DoS attacks, we propose an attack-resilient
distributed NE seeking algorithm as

$\dot{\textbf{x}}^{i}=-\mathcal{R}^{T}_{i}\triangledown_{i}J_{i}(\textbf{x}^{i})+\textbf{e}^{i},\
t\in[t_{2k},t_{2k+2}),\ k=0,1,\cdots,$ (9)

where $\textbf{e}^{i}$ denotes the relative estimated errors under attacks

$\textbf{e}^{i}=\left\\{\begin{array}[]{c}-\kappa\sum^{N}_{j=1}a_{ij}(\textbf{x}^{i}-\textbf{x}^{j}),\
t\in[t_{2k},t_{2k+1}),\\\
\sum^{N}_{j=1}a^{\Psi(t)}_{ij}(\textbf{x}^{j}-\textbf{x}^{i}),\
t\in[t_{2k+1},t_{2k+2}),\end{array}\right.$ (10)

where $\kappa$ is a positive constant gain to be specified later, and the
$a^{\Psi(t)}_{ij}$, $t\in[t_{2k+1},t_{2k+2})$, dependent on an attack flag
$\psi(i,j,t)$, is defined as $a^{\Psi(t)}_{ij}=0$ if $\psi(i,j,t)=1$ or $-1$;
otherwise $a^{\Psi(t)}_{ij}=1$ if $\psi(i,j,t)=0$. The expression of
$\psi(i,j,t)$ is presented below. Set the initial value $\psi(i,j,t)=0$ for
$\forall(i,j)\in\bar{\mathcal{E}}$, and then each player $i$ updates the
attack flag $\psi(j,i,t)$ as follows

  1. 1.

if player $i$ can receive information from player $j$ at $t$, then
$\psi(i,j,t)=0$ and it sends its information to player $j$;

  2. 2.

if player $i$ cannot receive information from player $j$ at $t$, then
$\psi(i,j,t)=1$ and it sends $\psi(i,j,t)=1$ to player $j$;

  3. 3.

if player $i$ receives $\psi(i,j,t)=1$ which means it knows the attacking of
channel $(i,j)$, then denote $\psi(i,j,t)=-1$ and send its information to
player $j$.

###### Remark 4

To facilitate understandings of (10), Fig. 4 shows the schematic of time
sequences with and without DoS attacks. Intuitively, not all communication
networks are secured anytime in practice, while it is reasonable to secure
some original network. In the presence of DoS attacks, $\textbf{e}^{i}$ in
(10) relies on $\psi(i,j,t)$. Each player updates this attack flag once an
attack signal over communication channels is detected by certain devices or
mechanisms. The attack detection design is beyond the scope of this work.

###### Remark 5

In the absence of DoS attacks, $\textbf{e}^{i}$ in (10) becomes
$\textbf{e}^{i}=-\kappa\sum^{N}_{j=1}a_{ij}(\textbf{x}^{i}-\textbf{x}^{j})$,
which is a modified version of the design in [11] that requires a restrictive
graph coupling condition under an undirected graph. In contrast, an adjustable
proportional control gain $\kappa$ is introduced to enable a natural trade-off
between the control effort and graph connectivity under Assumption 1.

Figure 4: Schematic of time sequences with and without DoS attacks.

Next, denote the following stacked vectors and matrices

x $\displaystyle=\text{col}(\textbf{x}^{1},\cdots,\textbf{x}^{N}),\
\mathcal{R}=\text{diag}\\{\mathcal{R}_{1},\cdots,\mathcal{R}_{N}\\},$ (11) e
$\displaystyle=\text{col}(\textbf{e}^{1},\cdots,\textbf{e}^{N}),\
\textbf{F({x})}=\text{col}(\triangledown_{1}J_{1}(\textbf{x}^{1}),\cdots,\triangledown_{N}J_{N}(\textbf{x}^{N})).$

Combining (9)-(11) gives rise to the closed-loop system under attacks in the
sense of a compact form

$\dot{\textbf{x}}=\left\\{\begin{array}[]{c}-\mathcal{R}^{T}\textbf{
F}(\textbf{x})-\kappa(\mathcal{L}\otimes I_{n})\textbf{x},\ t\in\
[t_{2k},t_{2k+1}),\\\ -\mathcal{R}^{T}\textbf{
F}(\textbf{x})-(\mathcal{L}^{\Psi(t)}\otimes I_{n})\textbf{x},\ t\in\
[t_{2k+1},t_{2k+2}),\end{array}\right.$ (12)

where $\mathcal{L}$ is the Laplacian matrix of the original strongly connected
digraph, while $\mathcal{L}^{\Psi(t)}$ is the Laplacian matrix of various
potential graphs under attacks, and its zero eigenvalues may not be simple as
those graphs can be unconnected under attacks.

###  IV-B Stability Analysis with An Exponential Convergence Rate

Before presenting the main result, we show that in the absence of attacks, the
equilibrium of the system occurs when all players can reach consensus at the
NE.

###### Proposition 1

Consider the game over the directed communication graph $\bar{\mathcal{G}}$.
Then, under Assumptions 1-3, $\tilde{\textbf{x}}=1_{N}\otimes x^{*}$ is the NE
of the game on networks in the absence of DoS attacks if
$\textbf{F}(\tilde{\textbf{x}})=0_{n}$ (or
$\triangledown_{i}J_{i}(\tilde{\textbf{x}}^{i})=0_{n_{i}}$). At the NE,
estimated vectors of all players reach consensus and equal to the NE $x^{*}$.
Thus, players’ action components coincide with optimal actions
($\tilde{\textbf{x}}^{i}_{i}=x^{*}_{i}$).

###### Proof:

In the absence of DoS attacks, let $\tilde{\textbf{x}}$ be an equilibrium of
the system. Then, we have $0_{Nn}=-\mathcal{R}^{T}\textbf{
F}(\tilde{\textbf{x}})-\kappa(\mathcal{L}\otimes I_{n})\tilde{\textbf{x}},$
which implies that multiplying both sides by $1^{T}_{N}\otimes I_{n}$ yields

$0_{n}=-(1^{T}_{N}\otimes I_{n})\mathcal{R}^{T}\textbf{
F}(\tilde{\textbf{x}})-\kappa(1^{T}_{N}\otimes I_{n})(\mathcal{L}\otimes
I_{n})\tilde{\textbf{x}}.$ (13)

Since $1^{T}_{N}\mathcal{L}=0^{T}_{Nn}$ under Assumption 1, we obtain
$0_{n}=(1^{T}_{N}\otimes I_{n})\mathcal{R}^{T}\textbf{F}(\tilde{\textbf{x}})$.
Then, it follows from the notations of $\mathcal{R}$ and F in (11) that
$\textbf{F}(\tilde{\textbf{x}})=0_{n}$. Then, submitting it into (13) gives
rise to $(\mathcal{L}\otimes I_{n})\tilde{\textbf{x}}=0_{Nn}$. Hence, there
exists certain $\theta\in\mathbb{R}^{n}$ such that
$\tilde{\textbf{x}}=1_{N}\otimes\theta$ under Assumption 1. Then, it has
$\textbf{F}(1_{N}\otimes\theta)=0_{n}$ for each player $i\in\mathcal{V}$.
Thus, $\triangledown_{i}J_{i}(\theta_{i},\theta_{-i})=0_{n_{i}}$. That is,
$\theta$ is a unique NE of the game and $\theta=x^{*}$. Thus,
$\tilde{\textbf{x}}=1_{N}\otimes x^{*}$ and for $i,j\in\mathcal{V}$, we have
$\tilde{\textbf{x}}^{i}=\tilde{\textbf{x}}^{j}=x^{*}$ (NE of the game). ∎

Notice that in the presence of DoS attacks, the system becomes
$\dot{\textbf{x}}=-\mathcal{R}^{T}\textbf{
F}(\textbf{x})-(\mathcal{L}^{\Psi(t)}\otimes I_{n})\textbf{x}$. Due to the
fact that the existence and uniqueness of the NE $x^{*}$ are guaranteed under
Assumptions 2 and 3, then, following a similar analysis above can give rise to
$0_{n}=(1^{T}_{N}\otimes I_{n})\mathcal{R}^{T}\textbf{
F}(\tilde{\textbf{x}})-(1^{T}_{N}\otimes I_{n})(\mathcal{L}^{\Psi(t)}\otimes
I_{n})\tilde{\textbf{x}}$. However, the presence of DoS attacks will make
$1^{T}_{N}\mathcal{L}^{\Psi(t)}=0^{T}_{Nn}$ not hold and then, consensus
estimates on the NE cannot be reached under attacks because there may not have
correct information exchange among all players. Under such a situation,
$x^{*}$ may not be the NE. Next, the main task is an explicit analysis of the
frequency and duration of attacks to guarantee
$\tilde{\textbf{x}}^{i}=\tilde{\textbf{x}}^{j}=x^{*}$.

Next, we present the main result on the resilient distributed NE seeking on
networks under DoS attacks.

###### Theorem 1

Under Assumptions 1-3, Problem 1 can be solvable for any $\textbf{x}_{i}(0)$
under the proposed resilient distributed optimization algorithm in (9)-(10)
provided that for
$\kappa>\frac{1}{\lambda_{2}(\hat{\mathcal{L}})}(\frac{\iota^{2}}{\varepsilon}+\iota)$
and positive scalars $\lambda_{a},\lambda_{b},u$ to be determined later, the
following two attack-related conditions are satisfied:

(1). There exists constants $\eta^{\ast}\in(0,\lambda_{a})$ and $\mu>1$ so
that $T_{f}$ in the attack frequency Definition 2 satisfies the condition:

$T_{f}>T^{*}_{f}=2\ln(\mu)/\eta^{\ast},$ (14)

(2). There exist constants $\lambda_{a}$, $\lambda_{b}>0$ such that $T_{a}$ in
the attack duration Definition 3 that satisfies the condition:

$\hskip
20.00003ptT_{a}>T^{*}_{a}=(\lambda_{a}+\lambda_{b})/(\lambda_{a}-\eta^{\ast}).$
(15)

Moreover, the estimated states can converge to the NE with an exponential
convergence rate, i.e.,

$||\textbf{x}(t)-\tilde{\textbf{x}}||^{2}\leq\varsigma
e^{-\eta(t-t_{0})}||\textbf{x}(t_{0})-\tilde{\textbf{x}}||^{2},\ \forall
t_{0}\geq 0,$ (16)

where $\varsigma$ is a positive scalar and
$\eta=\lambda_{a}-(\lambda_{a}+\lambda_{b})/T_{a}-\eta^{*}>0$.

###### Proof:

The proof includes four steps:
Step i): when communication networks do not suffer from DoS attacks during
$[t_{2k},t_{2k+1})$, we first show that the NE seeking can be achieved
exponentially under a strongly connected digraph.

Now, we first make a coordinate transformation as

$\displaystyle\overrightarrow{\textbf{x}}$
$\displaystyle=(1_{N}\otimes\mathcal{S})\textbf{x}\in\mathbb{R}^{Nn},\
\mathcal{S}=\frac{1}{N}(1^{T}_{N}\otimes I_{n}),$ (17)
$\displaystyle\overleftarrow{\textbf{x}}$ $\displaystyle=(\mathcal{T}\otimes
I_{n})\textbf{x}\in\mathbb{R}^{Nn},\
\mathcal{T}=I_{N}-\frac{1}{N}(1_{N}1^{T}_{N}).$ (18)

Then, it follows from (17) that the average estimate of $\textbf{x}^{i}$ can
be described by
$\bar{\textbf{x}}=\frac{1}{N}\sum^{N}_{i=1}\textbf{x}^{i}=\frac{1}{N}(1^{T}_{N}\otimes
I_{n})\textbf{x}=\mathcal{S}\textbf{x}$. Further, for any
$\textbf{x}\in\mathbb{R}^{Nn}$, it can be decomposed as
$\textbf{x}=\overrightarrow{\textbf{x}}+\overleftarrow{\textbf{x}}$ with
$(\overrightarrow{\textbf{x}})^{T}\overleftarrow{\textbf{x}}=0$ and
$(\mathcal{L}\otimes I_{n})\overrightarrow{\textbf{x}}=0_{Nn}$ under
Assumption 1.

For stability analysis, we select the Lyapunov function as

$\displaystyle V_{a}(\textbf{x})$
$\displaystyle=\frac{\alpha}{2}(\textbf{x}-\tilde{\textbf{x}})^{T}(\textbf{x}-\tilde{\textbf{x}})=\frac{\alpha}{2}(\overrightarrow{\textbf{x}}+\overleftarrow{\textbf{x}}-\tilde{\textbf{x}})^{T}(\overrightarrow{\textbf{x}}+\overleftarrow{\textbf{x}}-\tilde{\textbf{x}})$
$\displaystyle=\frac{1}{2}\left[\begin{array}[]{l}\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}\\\
\ \ \
\overleftarrow{\textbf{x}}\end{array}\right]^{T}P_{a}\left[\begin{array}[]{l}\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}\\\
\ \ \ \overleftarrow{\textbf{x}}\end{array}\right],$ (23)

where $P_{a}=\text{diag}\\{\alpha I_{Nn},\alpha I_{Nn}\\}$, $\alpha>0$ is an
adjustable constant, and the fact that
$(\overrightarrow{\textbf{x}})^{T}\overleftarrow{\textbf{x}}=0$ and
$1^{T}_{N}\mathcal{T}=0^{T}_{N}$ are used.

The time derivative of $V_{a}(\textbf{x})$ is expressed as

$\dot{V}_{a}(\textbf{x})=-\alpha(\textbf{x}-\tilde{\textbf{x}})^{T}[\mathcal{R}^{T}\textbf{F}(\textbf{x})+\kappa(\mathcal{L}\otimes
I_{n})\textbf{x}].$ (24)

As
$0_{Nn}=-\mathcal{R}^{T}\textbf{F}(\tilde{\textbf{x}})-\kappa(\mathcal{L}\otimes
I_{n})\tilde{\textbf{x}}$, (24) can be expressed as

$\dot{V}_{a}(\textbf{x})=-\alpha(\textbf{x}-\tilde{\textbf{x}})^{T}[\mathcal{R}^{T}(\textbf{F}(\textbf{x})-\textbf{F}(\tilde{\textbf{x}}))+\kappa(\mathcal{L}\otimes
I_{n})(\textbf{x}-\tilde{\textbf{x}})].$ (25)

In light of $\tilde{\textbf{x}}=1_{N}\otimes x^{*}$,
$\textbf{x}=\overrightarrow{\textbf{x}}+\overleftarrow{\textbf{x}}$, and
$(\mathcal{L}\otimes I_{n})\overrightarrow{\textbf{x}}=0_{Nn}$ under
Assumption 1, then the first term in (25) becomes

$\displaystyle-(\textbf{x}-\tilde{\textbf{x}})^{T}\mathcal{R}^{T}[\textbf{F}(\textbf{x})-\textbf{F}(\tilde{\textbf{x}})]$
(26)
$\displaystyle=-(\overleftarrow{\textbf{x}})^{T}\mathcal{R}^{T}[\textbf{F}(\textbf{x})-\textbf{F}(\overrightarrow{\textbf{x}})]-(\overleftarrow{\textbf{x}})^{T}\mathcal{R}^{T}[\textbf{F}(\overrightarrow{\textbf{x}})-\textbf{F}(\tilde{\textbf{x}})]$
$\displaystyle\ \ \
-(\overrightarrow{\textbf{x}}-\tilde{\textbf{x}})^{T}\mathcal{R}^{T}[\textbf{F}(\textbf{x})-\textbf{F}(\overrightarrow{\textbf{x}})]-(\overrightarrow{\textbf{x}}-\tilde{\textbf{x}})^{T}\mathcal{R}^{T}[\textbf{F}(\overrightarrow{\textbf{x}})-\textbf{F}(\tilde{\textbf{x}})].$

It follows from Assumption 3 that according to the $\iota_{F}$-Lipschitz
continuity of $F$, it yields that
$\|F(\textbf{x})-F(\overrightarrow{\textbf{x}})\|\leq\iota_{F}\|\overleftarrow{\textbf{x}}\|$.
Further,
$\|\textbf{F}(\textbf{x})-\textbf{F}(\overrightarrow{\textbf{x}})\|\leq\iota_{\textbf{F}}\|\overleftarrow{\textbf{x}}\|$
for certain scalar $\iota_{\textbf{F}}>0$. In addition, since
$\|\mathcal{R}^{T}\|=1$,
$\textbf{F}(\overrightarrow{\textbf{x}})=F(\bar{\textbf{x}})$, and
$\textbf{F}(\tilde{\textbf{x}})=F(x^{*})=0$,

$\displaystyle-(\overleftarrow{\textbf{x}})^{T}\mathcal{R}^{T}[\textbf{F}(\overrightarrow{\textbf{x}})-\textbf{F}(\tilde{\textbf{x}})]$
$\displaystyle=-(\overleftarrow{\textbf{x}})^{T}\mathcal{R}^{T}(F(\bar{\textbf{x}})-F(x^{*}))$
$\displaystyle\leq\iota_{F}\|\overleftarrow{\textbf{x}}\|\|\bar{\textbf{x}}-x^{*}\|,$
(27)
$\displaystyle-(\overrightarrow{\textbf{x}}-\tilde{\textbf{x}})^{T}\mathcal{R}^{T}[\textbf{F}(\textbf{x})-\textbf{F}(\overrightarrow{\textbf{x}})]$
$\displaystyle=-(\bar{\textbf{x}}-x^{*})^{T}(\textbf{F}(\textbf{x})-\textbf{F}(\overrightarrow{\textbf{x}}))$
$\displaystyle\leq\iota_{\textbf{F}}\|\bar{\textbf{x}}-x^{*}\|\|\overleftarrow{\textbf{x}}\|,$
(28)

where the fact that $\mathcal{R}\overrightarrow{\textbf{x}}=\bar{\textbf{x}}$
and $\mathcal{R}\tilde{\textbf{x}}=x^{*}$ is used, and exploiting the
$\varepsilon$-strong monotonicity of $F$, we can obtain

$\displaystyle-(\overrightarrow{\textbf{x}}-\tilde{\textbf{x}})^{T}\mathcal{R}^{T}[\textbf{F}(\overrightarrow{\textbf{x}})-\textbf{F}(\tilde{\textbf{x}})]$
$\displaystyle=-(\bar{\textbf{x}}-x^{*})^{T}(F(\bar{\textbf{x}})-F(x^{*}))\leq-\varepsilon\|\bar{\textbf{x}}-x^{*}\|^{2}.$
(29)

In addition, the second term in (25) can be rewritten as

$\displaystyle-(\textbf{x}-\tilde{\textbf{x}})^{T}(\mathcal{L}\otimes
I_{n})(\textbf{x}-\tilde{\textbf{x}})=-(\overrightarrow{\textbf{x}}+\overleftarrow{\textbf{x}})^{T}(\mathcal{L}\otimes
I_{n})(\overrightarrow{\textbf{x}}+\overleftarrow{\textbf{x}})$
$\displaystyle=-\frac{1}{2}\overleftarrow{\textbf{x}}^{T}[(\mathcal{L}^{T}+\mathcal{L})\otimes
I_{n}]\overleftarrow{\textbf{x}}\leq-\lambda_{2}(\hat{\mathcal{L}})\|\overleftarrow{\textbf{x}}\|^{2},$
(30)

where $\lambda_{2}(\hat{\mathcal{L}})$ is a minimal positive eigenvalue of
$\hat{\mathcal{L}}=\frac{1}{2}(\mathcal{L}^{T}+\mathcal{L})$.

Let $\iota=\max\\{\iota_{F},\iota_{\textbf{F}}\\}$. Substituting (27)-(30)
into (25) gives

$\displaystyle\dot{V}_{a}(\textbf{x})$
$\displaystyle\leq-\alpha\\{(\kappa\lambda_{2}(\hat{\mathcal{L}})-\iota)\|\overleftarrow{\textbf{x}}\|^{2}+\varepsilon\|\bar{\textbf{x}}-x^{*}\|^{2}-2\iota\|\overleftarrow{\textbf{x}}\|\|\bar{\textbf{x}}-x^{*}\|\\}$
$\displaystyle=-\alpha\left[\begin{array}[]{l}\|\bar{\textbf{x}}-x^{*}\|\\\ \
\
\|\overleftarrow{\textbf{x}}\|\end{array}\right]^{T}\left[\begin{array}[]{cc}\varepsilon&-\iota\\\
-\iota&\kappa\lambda_{2}(\hat{\mathcal{L}})-\iota\end{array}\right]\left[\begin{array}[]{l}\|\bar{\textbf{x}}-x^{*}\|\\\
\ \ \|\overleftarrow{\textbf{x}}\|\end{array}\right],$ (37)
$\displaystyle=-\frac{1}{2}\left[\begin{array}[]{l}\|\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}\|\\\
\ \
\|\overleftarrow{\textbf{x}}\|\end{array}\right]^{T}Q_{a}\left[\begin{array}[]{l}\|\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}\|\\\
\ \ \|\overleftarrow{\textbf{x}}\|\end{array}\right],$ (42)

where the fact that
$\|\bar{\textbf{x}}-x^{*}\|=\frac{1}{\sqrt{N}}\|\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}\|$
is used and $Q_{a}=2\alpha\\\
\left[\begin{array}[]{cc}\frac{\varepsilon}{N}&-\frac{\iota}{\sqrt{N}}\\\
-\frac{\iota}{\sqrt{N}}&\kappa\lambda_{2}(\hat{\mathcal{L}})-\iota\end{array}\right]>0$
if
$\kappa>\frac{1}{\lambda_{2}(\hat{\mathcal{L}})}(\frac{\iota^{2}}{\varepsilon}+\iota)$.

Thus, it concludes that there exists positive scalars $\alpha,\epsilon,\iota$,
and a matrix $P_{a}$, so that for the Lyapunov function in (23),

$V_{a}(\textbf{x})=\frac{1}{2}\chi^{T}P_{a}\chi\Rightarrow\dot{V}_{a}(\textbf{x})\leq-\lambda_{a}V_{a}(\textbf{x}),\
t\in[t_{2k},t_{2k+1}),$ (43)

where
$\chi=\text{col}(\overrightarrow{\textbf{x}}-\tilde{\textbf{x}},\overleftarrow{\textbf{x}})$
and $\lambda_{a}=\lambda_{min}(Q_{a})/\lambda_{max}(P_{a})>0$.

Step ii): when considering the presence of DoS attacks during
$[t_{2k+1},t_{2(k+1)})$, we select the Lyapunov function as

$\displaystyle V_{b}(\textbf{x})$
$\displaystyle=\frac{\beta}{2}(\textbf{x}-\tilde{\textbf{x}})^{T}(\textbf{x}-\tilde{\textbf{x}})=\frac{1}{2}\chi^{T}P_{b}\chi,$
(44)

where $P_{b}=\text{diag}\\{\beta I_{Nn},\beta I_{Nn}\\}$, $\beta\neq\alpha$ is
a positive constant.

In the presence of DoS attacks, $\dot{\textbf{x}}=-\mathcal{R}^{T}\textbf{
F}(\textbf{x})-(\mathcal{L}^{\Psi(t)}\otimes I_{n})\textbf{x}$ by (12). Since
$x^{*}$ is the unique NE, $0_{Nn}=-\mathcal{R}^{T}\textbf{
F}(\tilde{\textbf{x}})-(\mathcal{L}^{\Psi(t)}\otimes
I_{n})\tilde{\textbf{x}}$. By adding this equation, it follows from (26)-(42)
that the time derivative of $V_{b}(\textbf{x})$ can be described by

$\displaystyle\dot{V}_{b}(\textbf{x})$
$\displaystyle=-\beta(\textbf{x}-\tilde{\textbf{x}})^{T}[\mathcal{R}^{T}(\textbf{F}(\textbf{x})-\textbf{F}(\tilde{\textbf{x}}))+(\mathcal{L}^{\Psi(t)}\otimes
I_{n})(\textbf{x}-\tilde{\textbf{x}})]$
$\displaystyle\leq-\beta\left[\begin{array}[]{l}\|\bar{\textbf{x}}-x^{*}\|\\\
\ \
\|\overleftarrow{\textbf{x}}\|\end{array}\right]^{T}\left[\begin{array}[]{cc}\varepsilon&-\iota\\\
-\iota&-\iota\end{array}\right]\left[\begin{array}[]{l}\|\bar{\textbf{x}}-x^{*}\|\\\
\ \ \|\overleftarrow{\textbf{x}}\|\end{array}\right]$ (51) $\displaystyle\ \ \
-\beta(\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}+\overleftarrow{\textbf{x}})^{T}(\mathcal{L}^{\Psi(t)}\otimes
I_{n})(\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}+\overleftarrow{\textbf{x}})$
$\displaystyle\leq\beta\left[\begin{array}[]{l}\|\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}\|\\\
\ \
\|\overleftarrow{\textbf{x}}\|\end{array}\right]^{T}\left[\begin{array}[]{cc}-\frac{\varepsilon}{N}&\frac{\iota}{\sqrt{N}}\\\
\frac{\iota}{\sqrt{N}}&\iota\end{array}\right]\left[\begin{array}[]{l}\|\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}\|\\\
\ \ \|\overleftarrow{\textbf{x}}\|\end{array}\right]$ (58) $\displaystyle\ \ \
+\beta\|\mathcal{L}^{\Psi(t)}\|\left[\begin{array}[]{l}\|\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}\|\\\
\ \
\|\overleftarrow{\textbf{x}}\|\end{array}\right]^{T}\left[\begin{array}[]{l}\|\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}\|\\\
\ \ \|\overleftarrow{\textbf{x}}\|\end{array}\right]$ (63)
$\displaystyle=\frac{1}{2}\left[\begin{array}[]{l}\|\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}\|\\\
\ \
\|\overleftarrow{\textbf{x}}\|\end{array}\right]^{T}Q_{b}\left[\begin{array}[]{l}\|\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}\|\\\
\ \ \|\overleftarrow{\textbf{x}}\|\end{array}\right]$ (68)

where
$Q_{b}=2\beta\left[\begin{array}[]{cc}c-\frac{\varepsilon}{N}&c+\frac{\iota}{\sqrt{N}}\\\
c+\frac{\iota}{\sqrt{N}}&c+\iota\end{array}\right]$ with
$c=\|\mathcal{L}^{\Psi(t)}\|$, may have both positive and negative
eigenvalues.

Thus, it concludes that there exists positive scalars $\beta,\epsilon,\iota$,
and a matrix $P_{b}$, so that for the Lyapunov function in (44),

$V_{b}(\textbf{x})=\frac{1}{2}\chi^{T}P_{b}\chi\Rightarrow\dot{V}_{b}(\textbf{x})\leq\lambda_{b}V_{b}(\textbf{x}),\
t\in[t_{2k+1},t_{2k+2}),$ (69)

where $\lambda_{a}=\sigma_{\max}(Q_{b})/\lambda_{min}(P_{b})>0$ and
$\sigma_{\max}(Q_{b})$ denotes the maximum singular value of $Q_{b}$.

Step iii): we analyze the exponential convergence of the closed-loop system
from a switching perspective [21, 22, 24, 23, 25].

Let $\delta(t)\in\\{a,b\\}$ be a switching signal. Then, we can choose

$V(t)=\left\\{\begin{array}[]{c}V_{a}(\textbf{x}),\ \text{if}\
t\in[t_{2k},t_{2k+1}),\\\ V_{b}(\textbf{x}),\text{ if}\
t\in[t_{2k+1},t_{2k+2}),\end{array}\right.$ (70)

where $V_{a}(\textbf{x})$ and $V_{b}(\textbf{x})$ are defined in (23) and
(44), respectively.

Suppose that $V_{a}$ is activated in $[t_{2k},t_{2k+1})$, while $V_{b}$ is
activated in $[t_{2k+1},t_{2k+2})$. Then, by (43) and (69), we have

$V(t)\leq\left\\{\begin{array}[]{c}e^{-\lambda_{a}(t-t_{2k})}V_{a}(t_{2k}),\
\text{if}\ t\in\ [t_{2k},t_{2k+1}),\\\
e^{\lambda_{b}(t-t_{2k+1})}V_{b}(t_{2k+1}),\ \text{ if}\ t\in\
[t_{2k+1},t_{2k+2}).\end{array}\right.$ (71)

The closed-loop system is switched at $t=t^{+}_{2k}$ or $t=t^{+}_{2k+1}$. Let
$\mu=\max\\{\lambda_{max}(P_{a})/\lambda_{min}(P_{b}),\lambda_{max}(P_{b})/\lambda_{min}(P_{a})\\}>1$,
and next, we discuss the following two cases:

Case a): If $t\in[t_{2k},t_{2k+1})$, it follows from (71) that

$\displaystyle V(t)$ $\displaystyle\leq
e^{-\lambda_{a}(t-t_{2k})}V_{a}(t_{2k})\leq\mu
e^{-\lambda_{a}(t-t_{2k})}V_{b}(t_{2k}^{-})$ $\displaystyle\leq\mu
e^{-\lambda_{a}(t-t_{2k})}[e^{\lambda_{b}(t_{2k}-t_{2k-1})}V_{b}(t_{2k-1})]$
$\displaystyle\leq\mu
e^{-\lambda_{a}(t-t_{2k})}e^{\lambda_{b}(t_{2k}-t_{2k-1})}[\mu
V_{a}(t_{2k-1}^{-})]$
$\displaystyle=\mu^{2}e^{-\lambda_{a}(t-t_{2k})}e^{\lambda_{b}(t_{2k}-t_{2k-1})}V_{a}(t_{2k-1}^{-})$
$\displaystyle\leq\mu^{2}e^{-\lambda_{a}(t-t_{2k})}e^{\lambda_{b}(t_{2k}-t_{2k-1})}[e^{-\lambda_{a}(t_{2k-1}-t_{2k-2})}V_{a}(t_{2k-2})]$
$\displaystyle\leq\cdots$
$\displaystyle\leq\mu^{2k}e^{-\lambda_{a}|\Xi_{s}(t_{0},t)|}e^{\lambda_{b}|\Xi_{a}(t_{0},t)|}V_{a}(t_{0}).$
(72)

Case b): If $t\in[t_{2k+1},t_{2(k+1)})$, it follows from (71) that

$\displaystyle V(t)$ $\displaystyle\leq
e^{\lambda_{b}(t-t_{2k+1})}V_{b}(t_{2k+1})\leq\mu
e^{\lambda_{b}(t-t_{2k+1})}V_{a}(t_{2k+1}^{-})$ $\displaystyle\leq\mu
e^{\lambda_{b}(t-t_{2k+1})}[e^{-\lambda_{a}(t_{2k+1}-t_{2k})}V_{a}(t_{2k})$
$\displaystyle\leq\mu
e^{\lambda_{b}(t-t_{2k+1})}e^{-\lambda_{a}(t_{2k+1}-t_{2k})}[\mu
V_{b}(t^{-}_{2k-1})]$
$\displaystyle\leq\mu^{2}e^{\lambda_{b}(t-t_{2k+1})}e^{-\lambda_{a}(t_{2k+1}-t_{2k})}[e^{\lambda_{b}(t_{2k-1}-t_{2k-2})}V_{b}(t_{2k-2})]$
$\displaystyle\leq\cdots$
$\displaystyle\leq\mu^{2k+1}e^{-\lambda_{a}|\Xi_{s}(t_{0},t)|}e^{\lambda_{b}|\Xi_{a}(t_{0},t)|}V_{a}(t_{0}).$
(73)

Step iv): we consider bounds on attack frequency and duration.

Notice that $N_{a}(t_{0},t)=k$ for $t\in[t_{2k},t_{2k+1})$ and $k+1$ for
$t\in[t_{2k+1},t_{2(k+1)})$. Thus, for $\forall t\geq t_{0},$ by (72) and
(73),

$V(t)\leq\mu^{2N_{a}(t_{0},t)}e^{-\lambda_{a}|\Xi_{s}(t_{0},t)|}e^{\lambda_{b}|\Xi_{a}(t_{0},t)|}V(t_{0}).$
(74)

Notice that for all $t\geq t_{0},$
$|\Xi_{s}(t_{0},t)|=t-t_{0}-|\Xi_{a}(t_{0},t)|$ and $|\Xi_{a}(t_{0},t)|\leq
T_{0}+(t-t_{0})/T_{a}$ by Definition 3. Thus, we have

$\displaystyle-\lambda_{a}(t-t_{0}-|\Xi_{a}(t_{0},t)|)+\lambda_{b}|\Xi_{a}(t_{0},t)|$
$\displaystyle=-\lambda_{a}(t-t_{0})+(\lambda_{a}+\lambda_{b})|\Xi_{a}(t_{0},t)|$
$\displaystyle\leq-\lambda_{a}(t-t_{0})+(\lambda_{a}+\lambda_{b})[T_{0}+(t-t_{0})/T_{a}].$
(75)

Next, substituting (75) into (74) yields

$\displaystyle V(t)$
$\displaystyle\leq\mu^{2N_{a}(t_{0},t)}e^{-\lambda_{a}(t-t_{0}-|\Xi_{a}(t_{0},t)|)}e^{\lambda_{b}|\Xi_{a}(t_{0},t)|}V(t_{0})$
(76) $\displaystyle\leq
e^{(\lambda_{a}+\lambda_{b})T_{0}}e^{-\lambda_{a}(t-t_{0})}e^{\frac{(\lambda_{a}+\lambda_{b})}{\tau_{a}}(t-t_{0})}e^{2\ln(\mu)N_{a}(t_{0},t)}V(t_{0}).$

By exploiting the attack condition in (14), we can have

$2\ln(\mu)N_{a}(t_{0},t)\leq 2\ln(\mu)N_{0}+\eta^{*}(t-t_{0}).$ (77)

Let $\eta=\lambda_{a}-(\lambda_{a}+\lambda_{b})/T_{a}-\eta^{\ast}>0$. Based on
another attack condition in (15), and using (77), we can rewrite (76) as

$V(t)\leq e^{(\lambda_{a}+\lambda_{b})T_{0}+2\ln(\mu)N_{0}}\
e^{-\eta(t-t_{0})}V(t_{0}).$ (78)

Further, it follows from (23), (44), and (78) that

$||\chi(t)||^{2}\leq\varsigma e^{-\eta(t-t_{0})}||\chi(t_{0})||^{2},$ (79)

where
$\varsigma=e^{(\lambda_{a}+\lambda_{b})T_{0}+2\ln(\mu)}\varsigma_{a}/\varsigma_{b}$,
$\varsigma_{a}=\max\\{\lambda_{max}(P_{a}),\lambda_{max}(\\\ P_{b})\\}$, and
$\varsigma_{b}=\min\\{\lambda_{min}(P_{a}),\lambda_{min}(P_{b})\\}$.

Therefore, it follows from (79) that all estimate states
$\overrightarrow{\textbf{x}}-\tilde{\textbf{x}}$ and
$\overleftarrow{\textbf{x}}$ are bounded, and converge to zero exponentially.
Furthermore,
$\lim_{t\rightarrow\infty}(\overrightarrow{\textbf{x}}-\tilde{\textbf{x}})=0_{Nn}$
and $\lim_{t\rightarrow\infty}\overleftarrow{\textbf{x}}=0_{Nn}$. Then,
according to the coordinate transformation $\overrightarrow{\textbf{x}}$,
$\overleftarrow{\textbf{x}}$ in (17) and (18), and by using the fact that
$\textbf{x}=\overrightarrow{\textbf{x}}+\overleftarrow{\textbf{x}}$, we can
obtain that x exponentially converges to $\tilde{\textbf{x}}=1\otimes x^{*}$.
That is, Problem 1 is solved. ∎

###### Remark 6

Theorem 1 presented the main resilient distributed NE seeking result with an
exponential convergence rate
$\eta=\lambda_{a}-(\lambda_{a}+\lambda_{b})/T_{a}-\eta^{\ast}$. Here,
$(1-\frac{1}{T_{a}})\lambda_{a}$ is mainly used to measure the average rate of
exponential decay of the stable subsystems, while $\lambda_{b}/T_{a}$ isy used
to measure the exponential growth rate of unstable subsystems. The $\eta^{*}$
explains the exponential growth due to switchings. Note that the convergence
rate is not only affected by $\lambda_{a}$ and $\lambda_{b}$ that rely on the
communication topology, number of players, and control gains, but the attack
frequency and duration. Moreover, the larger values of attack frequency and
duration are, the more active that those attacks are allowed to be.

Notice that in the absence of DoS attacks, Theorem 1 can be reduced to the
following corollary:

Corollary 1: Under Assumptions 1-3, the following distributed NE seeking
algorithm enables all players’ estimated strategies to exponentially converge
to the NE provided
$\kappa>(\frac{\iota^{2}}{\varepsilon}+\iota)/\lambda_{2}(\hat{\mathcal{L}})$.

$\dot{\textbf{x}}^{i}=-\mathcal{R}^{T}_{i}\triangledown_{i}J_{i}(\textbf{x}^{i})+\textbf{e}^{i},\
\textbf{e}^{i}=-\kappa\sum_{j=1}^{N}a_{ij}(\textbf{e}^{i}-\textbf{e}^{j}),\
i\in\mathcal{V}.$

###### Remark 7

The corollary can cover some existing results (e.g., [8, 11]) as special
cases. Moreover, it can avoid restrictive graph coupling conditions in [11] by
adding a proportional gain $\kappa$, and remove the use of two-timescale
singular perturbation that yields semi-global convergence in [8, 11]. The
design does not require any initial requirements and allow for a general
directed graph.

##  V Numerical Simulations

In this section, numerical examples are presented to verify the effectiveness
of the proposed NE seeking design. The communication graph for all players in
examples is depicted in Fig. 5, in which the original strongly connected
digraph and the paralyzed graphs induced by attacks are presented,
respectively.

Figure 5: Communication graph for the players in the examples: (a) original
strongly connected digraph; and (b)-(d) paralyzed graphs under various DoS
attacks.

###### Example 1

(Energy Consumption Game)

In this example, we consider an energy consumption game of $N$ players for
Heating Ventilation and Air Conditioning (HVAC) system (see [8]), where the
cost function of each player $i$ can be modeled by the following function:

$J_{i}(x_{i},x_{-i})=a_{i}(x_{i}-b_{i})^{2}+\left(c\sum_{j=1}^{N}x_{j}+d\right)x_{i},\
i\in\mathcal{V},$

where $a_{i}>0,c>0$, $b_{i}$ and $d$ are constants for $i\in\mathcal{V}$. It
can be verified that Assumptions 2 and 3 are satisfied. Throughout this
simulation, let $a_{i}=1$, $c=0.1$, $d=10$ for each player. In the following
simulation, we investigate the effectiveness of the proposed resilient
distributed NE seeking algorithm in (9)-(10) from the perspective of network
under DoS attacks with several comparisons with the existing algorithm,
controller gain, attack frequency/duration, network topology and number of
players.

###  V-A Resilient Algorithm for Exponential Distributed NE Seeking

We consider five players ($N=5$) in the game over a strongly connected digraph
in the absence of attacks and three types of disconnected digraphs caused by
DoS attacks as shown in Fig. 5. Constants $b_{i}$ for $i=1,\cdots,5$, are set
to $10$, $15$, $20$, $25$, and $30$, respectively. By certain calculation
based on those parameters, the NE is
$x^{*}=\text{col}(2.0147,6.7766,11.5385,16.3004,21.0623)$ [8]. The initial
states are given by $x^{i}_{i}(0)=\text{col}(-2,-4,-6,-8,-10)$ and
$x^{i}_{j}(0)=\text{col}(15,10,5,0)$, $\forall i\neq j$, which are not close
to $x^{*}$. The control gain of algorithm in (10) is set as $\kappa=10$ and
the variable to balance the weight is
$\omega=\text{col}(\frac{1}{6},\frac{2}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6})$.

Next, we perform the proposed resilient NE seeking algorithm in (9)-(10), and
simulation results are provided as shown in Fig. 6. In particular, Fig. 6(a)
shows the occurrence of DoS attacks, where attack frequency and duration
conditions in Theorem 1 are satisfied. Fig. 6(b) illustrates all players’
estimate strategies on the NE $x^{*}$, while the relative errors of all
players’ actions $\|\textbf{x}-\textbf{x}^{*}\|/\\\ \|\textbf{x}^{*}\|$ are
depicted in Fig. 6(c). As observed, all players’ estimate strategies reach
consensus and converge to the NE exponentially.

Figure 6: Simulated results of the proposed attack-resilient NE seeking
algorithm in (9)-(10): (a) DoS attacks on various communication channels; (b)
all players’ estimated strategies $x^{i}_{j}(t),i,j\in\mathcal{V}$; and (c)
relative errors of all players’ actions.

###  V-B Algorithm Comparison

In order to make some comparisons, we perform the algorithm in [11] to further
illustrate the proposed algorithm’s effectiveness under DoS attacks. All
simulation environments are set the same as those in the subsection V-A. It
can be observed from Fig. 7 (a) that in the presence of attacks, the design in
[11] cannot guarantee the exact convergence of all players’ estimates to the
NE $x^{*}$. In contrast, Fig. 7 (b) shows the performance of the proposed
attack-resilient algorithm, which can verify the design’s effectiveness.

---

(a) Algorithm in [11]

(b) Algorithm in (9)-(10)

Figure 7: The plot of players’ strategies $x_{i}(t)$ produced by the proposed
algorithm and the algorithm in [16] in the presence of DoS attacks.

###  V-C Performance Analysis of Algorithm

1) Controller gain: we show the influence of control gain $\kappa$ on the
performance of algorithm, we conduct the proposed algorithm with the same
simulation setting in the subsection V-A, but under two different gains
$\kappa=1$ and $\kappa=5$, respectively. Fig. 8 shows the plots of players’
strategies and relative errors under different control gains, which implies
that the larger the gain $\kappa$, the better the convergence performance is,
which is as analyzed.

2) Number of Player: we increase the number of players to $N=3$, $5$, $10$ and
analyze its influence on the performance of the proposed algorithm. We set
$b_{i}=5i+5$ for each $i=1,2,\cdots,N$, and select a cycle directed graph as
the original communication graph. Fig. 9 depicts the performance of the
proposed algorithm under the different number of players. The algorithm is
scalable to various number of players, and the smaller the number of players,
the better the performance is as expected.

3) Network Topology: in this part, we investigate the influence of the
different original communication topologies on the performance of algorithm.
As analyzed in the theorem, the larger the $\lambda_{a}$, the better
convergence performance is. Intuitively, the larger the nonzero eigenvalue of
the Laplacian matrix, the larger the $\lambda_{a}$. Fig. 10 shows the
logarithmic curve of relative errors. As we see, the performance is better for
a graph with more links.

4) Attack Frequency and Duration: we illustrate the influence of attack
frequency and duration on the performance of algorithm. According to the
theorem, the frequency and duration of attacks have to be constrained to
guarantee the convergence of all players’ estimates to the NE. Simulated
result is shown in Fig. 11, where a comparison, when the attack frequency and
duration conditions in (14) and (15) hold and do not not, is provided. As can
be seen, the result is consistent with the analysis in Remark 6.

---

(a) $\kappa=1$

(b) $\kappa=5$

Figure 8: The plot of players’ strategies $x_{i}$ and their relative errors
$\|\textbf{x}-\textbf{x}^{*}\|/\|\textbf{x}^{*}\|$ produced by the proposed
algorithm in (9)-(10) with different controller gains.

---

(a)

(b)

(c)

Figure 9: Performance of the proposed algorithm (9)-(10) with different
players: (a) the cycle digraph; (b) the relative error; and (c) its
logarithmic curve.

---

(a)

(b)

Figure 10: The logarithmic curve of the relative error produced by the
proposed algorithm (9)-(10) with different network graphs.

---

(a) Small attack requency and duration

(b) Large attack requency and duration

Figure 11: The plot of players’ strategies $x_{i}(t)$ produced by the proposed
algorithm: (a) conditions (14)-(15) hold; (b) conditions (14)-(15) are not
satisfied.

###### Example 2

(General Non-Quadratic Game)

In this example, we investigate a more general non-quadratic game, where the
cost functions for each player $i$ are given by

$\displaystyle J_{1}(x_{1},x_{-1})$
$\displaystyle=\frac{x^{2}_{1}}{2}+x_{1}\sum^{5}_{j=2}x_{j},J_{2}(x_{2},x_{-2})=\frac{e^{\frac{x_{2}}{2}}}{2}+x_{2}x_{4},$
$\displaystyle J_{3}(x_{3},x_{-3})$
$\displaystyle=\frac{x^{2}_{3}}{2}+x_{1}^{3},\
J_{4}(x_{4},x_{-4})=\ln(e^{x_{4}})+x^{2}_{4}+x^{3}_{3},$ $\displaystyle
J_{5}(x_{5},x_{-5})$
$\displaystyle=x^{2}_{5}-5x_{5}+x^{3}_{1}x_{2}+x_{3}x^{4}_{4},\ i=1,\cdots,5.$
(80)

Next, we perform the proposed algorithm for this non-quadratic game with the
same simulation setting in the subsection V-A. By the calculation, the NE is
$x^{*}=\text{col}(-4.6589,4.1589,0,-2,2.5)$. The simulation result is shown in
Fig. 12, and as we can see, under attacks, all players’ estimates can reach
consensus and converge to the NE exponentially.

Figure 12: Simulated results of the proposed attack-resilient NE seeking
algorithm in (9)-(10): (a) DoS attacks on various communication channels; (b)
all players’ estimated strategies $x^{i}_{j}(t),\ i,j\in\mathcal{V}$; and (c)
relative errors of all players’ actions.

##  VI Conclusion

In this paper, an attack-resilient distributed algorithm has been presented
for exponential NE seeking of non-cooperative games, where all players’
strategies have been updated through a directed communication network subject
to malicious DoS attacks. Under such an adversary network environment, the
exponential convergence of the proposed distributed algorithm has been
established through the explicit analysis of the attack frequency and
duration. Moreover, in the absence of DoS attacks, the corollary has been
provided, which can cover many existing results as special cases. The
effectiveness of the developed approach has been illustrated by the numerical
examples. Further work may consider distributed NE seeking problems for
aggregative games with constraints.

## References

  * [1] S. Li, T. Basar, “Distributed algorithms for the computation of noncooperative equilibria,” _Automatica_ , 23, 523–533, 1987.
  * [2] J. S. Shamma, G. Arslan, “Dynamic fictitious play, dynamic gradient play, and distributed convergence to Nash equilibria,” _IEEE Trans. Autom. Control_ , 50(3): 312–327, 2005.
  * [3] M. S. Stankovic, K. H. Johansson, D. M. Stipanovic, “Distributed seeking of Nash equilibria with applications to mobile sensor networks,” _IEEE Trans. Autom. Control_ , 57(4): 904–919, 2012.
  * [4] P. Frihauf, M. Krstic, T. Basar, “Nash equilibrium seeking in noncooperative games,” _IEEE Trans. Autom. Control_ , 57(5): 1192–1207, 2012.
  * [5] S. Grammatico, “Dynamic control of agents playing aggregative games with coupling constraints,” _IEEE Trans. Autom. Control_ , 62(9): 4537–4548, 2017.
  * [6] G. Scutari, F. Facchinei, J. Pang, D. Pallomar, “Real and complex monotone communication games,” _IEEE Trans. Inf. Theory_ , 60(7): 400–409, 2014.
  * [7] J. R. Marden, “Learning in large-scale games and cooperative control,” _Ph.D. dissertation, University of California_ , Los Angeles, CA, USA, 2007.
  * [8] M. Ye, G. Hu, “Distributed Nash equilibrium seeking by a consensus based approach,” _IEEE Trans. Autom. Control_ , 62(9): 4811–4818, 2017.
  * [9] C. Sun, G. Hu, “Distributed Nash equilibrium seeking for a generalized convex game with nonsmooth objective functions and certain nonsmooth constraints,” _International Conference on Information and Automation_ , pp. 77–82, Fujian, China, Aug. 11–14, 2018.
  * [10] K. Lu, G. Jing, L. Wang, “Distributed algorithms for searching generalized Nash equilibrium of noncooperative Games,” _IEEE Trans. Cybern.,_ 49(6): 2362–2371, 2019.
  * [11] D. Gadjov, L. Pavel, “A passivity-based approach to Nash equilibrium seeking over networks,” _IEEE Trans. Autom. Control_ , 64: 1077–1092, 2019.
  * [12] S. Liang, P. Yi, Y. Hong, “Distributed Nash equilibrium seeking for aggregative games with coupled constraints,” _Automatica_ , 85: 179–185, 2017.
  * [13] Z. Deng, X. Nian, “Distributed generalized Nash equilibrium seeking algorithm design for aggregative games over weight-balanced digraphs,” _IEEE Trans. Neur. Net. Lear._ , 30(3): 695–706, 2019.
  * [14] Y. Zhang, X. Wang, H. Ji, “Distributed Nash equilibrium seeking in aggregative game with disturbance rejection,” _the 15th International Conference on Control and Automation_ , Scotland, July 16-19, pp. 1091–1095, 2019.
  * [15] X. Wang, N. Xiao, T. Wongpiromsarn, L. Xie, E. Frazzoli, D. Rus, “Distributed consensus in noncooperative congestion games: an application to road pricing,” _in IEEE Int. Conf. Control Autom._ , pp. 1668–1673, 2013.
  * [16] F. Salehisadaghiani, L. Pavel, “Distributed Nash equilibrium seeking: a gossip-based algorithm,” _Automatica_ , 72(1): 209–216, 2016.
  * [17] F. Salehisadaghiani, W. Shi, L. Pavel, “Distributed Nash equilibrium seeking under partial-decision information via the alternating direction method of multipliers,” _Automatica_ , 103: 27–35, 2019.
  * [18] Y. Pang, G. Hu, “Distributed Nash equilibrium seeking with limited cost function knowledge via a consensus-based gradient-free method,” _IEEE Trans. Autom. Control_ , DOI: 10.1109/TAC.2020.2995666, 2020.
  * [19] Y. Pang, G. Hu, “Randomized gradient-free distributed optimization methods for a multi-agent system with unknown cost function,” _IEEE Trans. Autom. Control_ , 65(1), pp. 333–340, 2020.
  * [20] Y. Mo, J. P. Hespanha, B. Sinopoli, “Resilient detection in the presence of integrity attacks,” _IEEE Trans. Signal Process_. 62(1): 31–43, 2014.
  * [21] Z. Feng, G. Hu, “Distributed coordinated control for multi-agent systems under two types of attacks with an application to power system,” _the 19 $\mathit{th}$ IFAC World Congress,_ August 24-29, South Africa, pages, 124–-130, 2014.
  * [22] Z. Feng, G. Hu, G. Wen, “Distributed consensus tracking for multi-agent systems under two types of attacks,” _Int. J. Robust. Nonlinear Control_ , 26(5) : 896–918, 2015.
  * [23] Z. Feng, G. Wen, G. Hu, “Distributed secure control for multi-agent systems under strategic attacks,” _IEEE Trans. Cybern.,_ 47(5): 1273–1284, 2017.
  * [24] Z. Feng, G. Hu, “Distributed secure average consensus for linear multi-agent systems under DoS attacks,” _American Control Conference,_ May 24-26, Seattle, WA, USA, 2017.
  * [25] Z. Feng, G. Hu, “Secure cooperative event-triggered control of linear multi-agent systems under DoS attacks,” _IEEE Trans. Contr. Syst. Tech._ , 28(3): 741–752 2020.
  * [26] W. Xu, G. Hu, D. W. C. Ho, Z. Feng, “Distributed secure cooperative control under Denial-of-Service attacks from multiple adversaries,” _IEEE Trans. Cybern.,_ 50(2): 3458–3467, 2020.
  * [27] S. Feng, P. Tesi, “Resilient control under denial-of-service: Robust design,” _Automatica_ , 79: 42–51, 2017.
  * [28] A. Lu, G. Yang, “Input-to-state stabilizing control for cyber-physical systems with multiple transmission channels under denial of service,” _IEEE Trans. Autom. Control_ , 63(6): 1813–1820, 2018.
  * [29] B. Chen, D. Ho, W. Zhang, L. Yu, “Distributed dimensionality reduction fusion estimation for cyber-physical systems under DoS attacks,” _IEEE Trans. Syst., Man, Cybern., Syst._ , 49(2): 455–468, 2019.
  * [30] D. Zhang, L. Liu, G. Feng, “Consensus of heterogeneous linear multi-agent systems subject to aperiodic sampled-data and DoS attack,” _IEEE Trans. Cybern._ , 49(4): 1501–1511, 2019.
  * [31] L. Su, N. Vaidya, “Multi-agent optimization in the presence of Byzantine adversaries: Fundamental limits,” _American Control Conference_ , July 6-8, Boston, MA, USA, pages 7183–-7188, 2016\.
  * [32] S. Sundaram, B. Gharesifard, “Distributed optimization under adversarial nodes,” _IEEE Trans. Autom. Control,_ 64(3): 1063–1076, 2019.
  * [33] C. Zhao, J. He, Q. Wang, “Resilient distributed optimization algorithm against adversarial attacks,” _IEEE Trans. Autom. Control,_ DOI: 10.1109/TAC. 2019.2954363, 2020.
  * [34] Z. Feng, G. Hu, C. G. Cassandras, “Finite-time distributed convex optimization for continuous-time multi-agent systems with disturbance rejection,” _IEEE Trans. Control Netw. Syst.,_ 7(2): 686–698, 2020.
  * [35] T. Başar, G. J. Olsder. Dynamic noncooperative game theory. SIAM, 1995.
  * [36] H. K. Khalil, Nonlinear systems. 3rd ed, Prentice-Hall, 2002.
  * [37] Z. Li, Z. Duan. Cooperative control of multi-agent systems: a consensus region approach. CRC Press, 2014.
  * [38] E. Ahmed, H. N. Agiza, “Dynamics of a Cournot game with ncompetitors,” _Chaos Solitons & Fractals_, 9(9): 1513–1517, 1998.