# Resilient Sparse Controller Design with
Guaranteed Disturbance Attenuation

MirSaleh Bahavarnia§ and Hossein K. Mousavi† $\S$ M. Bahavarnia is with the
Department of Electrical and Computer Engineering and the Institute for
Systems Research at the University of Maryland, College Park, MD 20742, USA
(e-mail: mbahavar@umd.edu).$\dagger$ H.K. Mousavi is with the Department of
Mechanical Engineering and Mechanics, Lehigh University, Bethlehem, PA 18015,
USA (e-mail: mousavi@lehigh.edu).

###### Abstract

We design resilient sparse state-feedback controllers for a linear time-
invariant (LTI) control system while attaining a pre-specified guarantee on
${\mathcal{H}}_{\infty}$ performance measure. We leverage a technique from
non-fragile control theory to identify a region of resilient state-feedback
controllers. Afterward, we explore the region to identify a sparse controller.
To this end, we use two different techniques: the greedy method of
sparsification, as well as the re-weighted $\ell_{1}$ norm minimization. Our
approach highlights a tradeoff between the sparsity of the feedback gain,
performance measure, and fragility of the design. To best of our knowledge,
this work is the first framework providing performance guarantees for sparse
feedback gain design.

##  I Introduction

During the last two decades, the robust resilient control has been an active
branch of research in control theory [1, 2, 3, 4]111This field of research is
alternatively called robust non-fragile control.. In this context, it is of
paramount interest to design controllers that are not only robust against the
exogenous processes such as disturbance or feedback noises but also resilient
to perturbations in the controller space. Designing controllers with this
property is undoubtedly necessary because in many real-world applications we
always face uncertainties and perturbations in the controller architecture.

On the other hand, sparse control design has amused the control theorists in
recent years [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. In sparse
control design, the main objective is to strike a balance between the
performance and structure of the controller. In many cases, the goal is to
push more elements of the feedback gains to be zero, while the performance of
the resulting design is satisfactory. This objective may arise from a need to
promote decentralization in large-scale systems wherein traditional (or
classic) centralized controllers may no longer be applicable. From another
point of view, sparse feedback gains reduce the computational burden in the
implementation of control protocols. The research in this area has found
applications, for instance in synchronization networks [18], dynamic mode
decomposition [19], and adaptive optics systems [20].

In this paper, we focus on a rather general idea to come up with a framework
to find controllers that are simultaneously robust, resilient, and sparse.
First, we use the tools from [21] to find a continuum of resilient controllers
with a guaranteed $\mathcal{H}_{\infty}$ performance measure. Then, we uncover
sparsified controllers that lie inside the specified geometric bounds. Our
approach allows us to tune the level of sparsity using a scalar parameter,
which highlights the tradeoff between the desired levels of sparsity and the
performance and non-fragility of the feedback controller design. To the best
of our knowledge, this paper is the first research paper that provides
performance guarantees on the designed sparse controller. We have included
several numerical experiments in this paper to support the theoretical
contributions. While our approach does not guarantee a certain level of
sparsity in the designed controller, the numerical examples create
meaningfully sparse solutions.

Notations: Throughout the paper, we adopt the following notations: all the
vectors and matrices are represented by lower and upper case letters,
respectively. The transpose of a matrix $X$ is denoted by $X^{T}$. The set of
$n\times n$ symmetric matrices is denoted by $\mathbb{S}_{n}$. The identity
matrix of appropriate size is denoted by $I$. The positive definite matrix
inequality operators are shown by $\prec$ and $\succ$. The $2$-norm of a
matrix is represented by $\|.\|_{2}$. The $\ell_{0}$ measure of a matrix is
denoted by $\|.\|_{0}$, which is equal to the number of nonzero elements of
the matrix. The $p$-norms of the matrices for $p\in\mathbb{Z}_{++}$ is denoted
by $\|.\|_{p}$ where $\mathbb{Z}_{++}$ represents the set of all positive
integers. The element-wise Hadamard matrix product operator is shown by
$\circ$. A normal random variable with mean $a$ and covariance $\Sigma$ is
denoted by $\mathcal{N}(a,\Sigma)$. The maximal eigenvalue of matrix
$M\in\mathbb{S}_{n}$ is denoted by $\lambda_{\max}(M)$. The maximum singular
value of a matrix $M$ is denoted by $\sigma_{\max}(M)$. The supremum of a set
is denoted by $\sup$.

##  II Problem Statement

We consider a plant with linear time-invariant (LTI) dynamics that are given
by

$\displaystyle\left\\{\begin{array}[]{l}\dot{x}(t)=Ax(t)+Bu(t)+B_{v}v(t)\vspace{2mm}\\\
y(t)=Cx(t)+D_{gu}u(t)+D_{gv}v(t)\end{array}\right.,$ (3)

where $x(t)\in\mathbb{R}^{n}$, $u(t)\in\mathbb{R}^{m}$,
$v(t)\in\mathbb{R}^{m_{v}}$, and $y(t)\in\mathbb{R}^{p}$ denote the state
vector, control input, disturbance, and output, respectively. One may choose
to control the system using a static state-feedback controller

$u(t)=Fx(t),$ (4)

for state-feedback controller222or alternatively, a feedback gain
$F\in\mathbb{R}^{m\times n}$, such that not only the closed-loop stability is
achieved, but also desired performance objectives and constraints are
satisfied. In this paper, we aim at finding state-feedback designs that are
robust to disturbances, sparse, and resilient to the uncertainties in the
implementation. In what follows, we elaborate and formalize these objectives.

Robustness: in plants with an external disturbance, the $\mathcal{H}_{\infty}$
norm is a common measure to describe the quality of the disturbance
attenuation in the output; that is

$\displaystyle\|G(s)\|_{\mathcal{H}_{\infty}}:=\sup\limits_{\omega}\leavevmode\nobreak\
\sigma_{\max}\left(G(j\omega)\right),$ (5)

where $G(s)\in\mathbb{R}^{p\times m_{v}}$ is the transfer matrix from the
disturbance $v$ to the output $y$[22]. A $\gamma$-level disturbance
attenuation concerns finding a controller $F$ such that

$\displaystyle\|G(s)\|_{\mathcal{H}_{\infty}}\leq\gamma.$ (6)

Sparsity: at the same time, we need controllers with the smallest number of
nonzero elements, (i.e., a sparse design); the direct measure to consider this
objective is the $\ell_{0}$ measure $\|F\|_{0}.$ A sparse controller design
may require lower levels of computation or communication in the control
system. Moreover, the privacy/security concern could become less of an issue
in systems wherein the information for control action should be shared over a
medium such as a cloud [23] or sent over communication channels.

Resilience: in practice the implementation of any feedback controller always
suffers from uncertainty to a certain extent. This might arise due to the
round-off or modeling errors, actuator degradation, etc. Hence, we are
generally interested in designs that are resilient (or non-fragile) to such
perturbations.

These objectives motivate us to define the following problem, wherein
$\Omega_{F}$ is a region in the controller space, inside which, certain
performance criteria hold. This set characterizes the structure of controller
perturbations.

######  Problem 1 (Robust Resilient Sparse State-Feedback)

For a desired performance level $\gamma>0$, find a state-feedback controller
$F$, which solves

$\displaystyle\begin{array}[]{ll}\underset{F}{\textrm{minimize}}&\leavevmode\nobreak\
\|F\|_{0}\vspace{2mm}\\\ {\textrm{subject to:}}&\leavevmode\nobreak\
\|G(s)\|_{\mathcal{H}_{\infty}}\leq\gamma,\textrm{ \leavevmode\nobreak\
\leavevmode\nobreak\ for any }\hat{F}\in\Omega_{F}.\end{array}$

Paper Objectives: even in the absence of resilience requirement, (i.e., the
set is a singleton $\Omega_{F}\equiv\\{F\\}$), it can be shown that Problem 1
is NP-hard [24], so we will not seek the exact solution to this problem.
Therefore, as a relaxation, we propose a two-level algorithm to find a
feedback gain $F$ such that

(i) the controller $F$ is (hopefully) sparse.
(ii) the $\gamma$-level performance of the system is guaranteed.
(iii) we may evaluate an explicitly defined region $\Omega_{F}$ in the
controller space such that $F\in{\Omega_{F}}$ and that for any feedback gain
$\hat{F}\in\Omega_{F}$, the stability and $\gamma$-level performance of the
controller are still guaranteed.

In this paper, we develop a general two-stage method in which, first the
region $\Omega_{F}$ is identified using linear matrix inequalities. Then, a
sparse controller is found inside this region in the second stage. To realize
the second step, we develop two different ways that the sparsification may be
conducted: first, we leverage re-weighted $\ell_{1}$ minimization techniques,
which also result in a set of linear matrix inequalities. Alternatively, we
show that we can efficiently sparsify the controller by a greedy method. In
both cases, we illustrate and control the tradeoff between resilience and
sparsity of the controllers using an appropriate design parameter.

##  III Generic Two-Level Algorithm

First, we recall the machinery developed by Peacelle et. al. [4], which
characterizes a quadratic region of controllers inside which every controller
provides us with $\gamma$-level performance attenuation. Then, we define a
generic two-level algorithm to show that how we can use the output of that
method to find sparse state-feedback controllers that enjoy the robustness and
resilience of the first method.

###  III-A Quadratically Resilient Robust Control

In the following theorem, we recall an important result from [21], wherein the
authors introduce a method to find a region in the controller space such that
any controller lying in the region induces the $\gamma$-level disturbance
attenuation.

######  Theorem 1 (see [21])

the linear time-invariant control system (3) is stabilizable by means of a
state-feedback controller $F$ such that the transfer matrix from the
disturbance to the output satisfies
$\|G(s)\|_{\mathcal{H}_{\infty}}\leq\gamma$ if and only if there exists a
Lyapunov matrix $P$ and three matrices $\hat{X}\in\mathbb{S}_{m}$,
$\hat{Y}\in\mathbb{R}^{m\times n}$, and $\hat{Z}\in\mathbb{S}_{n}$, which
satisfy the linear matrix inequalities given by

$\displaystyle\leavevmode\nobreak\ \hat{X}\preceq 0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ P\succ
0,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \hat{Z}\succ
0,$ (7) $\displaystyle\begin{bmatrix}Q_{11}&Q_{12}\\\
Q_{21}&Q_{22}\end{bmatrix}+\begin{bmatrix}B_{v}\\\
D_{gv}\end{bmatrix}\begin{bmatrix}B_{v}\\\ D_{gv}\end{bmatrix}^{T}\prec 0,$

where the blocks are defined to be

$\displaystyle Q_{11}$
$\displaystyle:=AP-B\hat{Y}+\left(AP-B\hat{Y}\right)^{T}-B\hat{X}B^{T}+\hat{Z},$
$\displaystyle Q_{12}$ $\displaystyle:=\left(CP-
D_{gu}\hat{X}B^{T}-D_{gu}\hat{Y}\right)^{T}=Q_{21}^{T},$ $\displaystyle
Q_{22}$ $\displaystyle:=-\gamma^{2}I-D_{gu}\hat{X}D_{gu}^{T}.$

Moreover, for any feasible solution, if we define

$\displaystyle X:=\hat{X}+\hat{Y}\hat{Z}^{-1}\hat{Y}^{T},\leavevmode\nobreak\
Y:=\hat{Y}\hat{Z}^{-1}P,\leavevmode\nobreak\ Z:=P\hat{Z}^{-1}P,$
$\displaystyle F_{o}:=-YZ^{-1}=-\hat{Y}P^{-1},\leavevmode\nobreak\
R:=YZY^{T}-X,$

then, for any controller $F$ satisfying the quadratic inequality

$\displaystyle(F-F_{o})\,Z\,(F-F_{o})^{T}\,\preceq\,{R},$ (8)

the closed-loop systems is asymptotically stable and the performance measure
of the system satisfies $\|G(s)\|_{\mathcal{H}_{\infty}}\leq\gamma$.

We observe that the feedback gain $F_{o}$ is the center of the ellipsoidal
region defined by (8). In fact, the quadratic constraint (8) implies that all
controllers that are close enough to control gain $F_{o}$ will inherit the
stability and performance guarantee, (i.e., bound on $\mathcal{H}_{\infty}$
performance measure) from central controller $F_{o}$.

###  III-B Generic Search for Sparse Feedback Controllers

Here, we go over the idea for finding sparse controllers. Suppose that we have
found a quadratic resilient region; i.e., the interior of the matrix ellipsoid
centered at $F_{o}$ defined by quadratic inequality (8). Now, we wish to find
a sparse controller $F$ inside this region. As long as we do not leave this
region, we would benefit from the performance guarantee of
$\left\|G(s)\right\|_{\mathcal{H}_{\infty}}\leq\gamma$. Moreover, if we
maintain a ”safe distance” from the boundary of this matrix ellipsoid, we can
preserve the resilience of the design to certain extent. To do so, instead of
searching for sparse controllers in the whole region defined by (8), we
introduce a design parameter $\theta\in[0,1]$ and use it to control the
breadth of the search space around the initial controller $F_{o}$.
Mathematically, compared to the region defined by (8), we consider a possibly
shrunk region of controllers

$\displaystyle(F-F_{o})Z(F-F_{o})^{T}\preceq\theta{R},$ (9)

for design parameter $\theta\in[0,1]$. We are interested in rewriting this
inequality in a form that is linear in its argument $F$. This is particularly
useful for arriving at semi-definite programs [25]. To achieve this, one way
is to apply the Schur complement [26] on matrix inequality (9) to obtain an
equivalent linear matrix inequality

$\displaystyle\begin{bmatrix}\theta{R}&F-F_{o}\\\
(F-F_{o})^{T}&Z^{-1}\end{bmatrix}\succeq 0.$ (10)

Given this search space, we define the following generic problem for finding
sparse feedback gains.

###### Problem 2

Given matrix ellipsoid (8), find a sparse controller from the optimization
problem

$\displaystyle\begin{array}[]{ll}\underset{F}{\textrm{minimize}}&\|F\|_{0}\vspace{2mm}\\\
{\textrm{subject to:}}&\begin{bmatrix}\theta{R}&F-F_{o}\\\
(F-F_{o})^{T}&Z^{-1}\end{bmatrix}\succeq 0.\end{array}$

Depending on the value of $\theta$, Problem 2 will imply either of the
following designs with guaranteed $\mathcal{H}_{\infty}$ performance:

  1. 1.

$\theta=0$: quadratically resilient design, i.e., $F=F_{o}$.

  2. 2.

$\theta\in(0,1)$: resilient sparse design.

  3. 3.

$\theta=1$: sparse design.

We merge these two steps to form the following _generic_ two-stage algorithm:

Algorithm 1:
Generic Two-Stage Algorithm of Resilient Sparse State-Feedback Controller with
Guaranteed $\mathcal{H}_{\infty}$ Performance

  1. 1.

Find $Z,\leavevmode\nobreak\ R$, $F_{o}$ from linear matrix inequalities (7).

  2. 2.

Solve Problem 2 for a sparse controller $F$.

We can call this method the generic algorithm for finding a sparse controller
because at this point we have not stated how to perform the second step in
Problem 2. In the next section, we preserve the constraint of Problem 2 and
consider tractable relaxations for the objective function, which hopefully
give us controllers that are sparse.

##  IV Realization of Sparsity

We consider two methods to fill in the blanks of Algorithm $1$, i.e., methods
that let us find sparse feedback gains.

###  IV-A Re-Weighted $\ell_{1}$ Regularization

Instead of dealing directly with Problem 2, we can consider a weighted
$\ell_{1}$ regularization and define a problem:

###### Problem 3

Given the matrix ellipsoid (8), find a sparse controller from the optimization
problem,

$\displaystyle\begin{array}[]{ll}\underset{F}{\textrm{minimize}}&\|W\circ
F\|_{1}\vspace{2mm}\\\ {\textrm{subject
to:}}&\begin{bmatrix}\theta{R}&F-F_{o}\\\
(F-F_{o})^{T}&Z^{-1}\end{bmatrix}\succeq 0,\end{array}$

where $W\in\mathbb{R}^{m\times n}$ represents the element-wise positive weight
matrix.

In fact, the $\ell_{1}$ term in Problem 3 is expected to act as a _proxy_ to
the sparsity [7]. Once we relax Problem 2 and replace it with Problem 3, we
get Algorithm $2$.

Algorithm 2:
Resilient Sparse State-Feedback Controller Design with $\mathcal{H}_{\infty}$
Performance using $\ell_{1}$ Regularization

  1. 1.

Find $Z,\leavevmode\nobreak\ R$, $F_{o}$ from the linear matrix inequalities
(7).

  2. 2.

Solve Problem 3 for a sparse controller $F$.

Regularization with Re-weighting: The choice of the weight matrix $W$ plays a
significant role in the properties of the proposed method. When an appropriate
weight matrix is not available, the re-weighted $\ell_{1}$ norm technique can
also be utilized. In this method, the weight assigned to each element is
updated iteratively. This element-wise update is inversely proportional to the
absolute value of the corresponding element in $F^{(k)}$ that is recovered
from the past iteration as

$\displaystyle W_{ij}^{(k+1)}\,{\leftarrow}\,\frac{1}{|F_{ij}^{(k)}|+\zeta},$
(11)

where the constant $\zeta>0$ is added to the denominator of update law (11)
for numerical stability [27]. The value of $\zeta$ is usually chosen to be
relatively small. The stopping criterion is implemented by examining the ratio

$\epsilon^{(k+1)}\,{\leftarrow}\,\frac{\|F^{(k+1)}-F^{(k)}\|_{2}}{\|F^{(k+1)}\|_{2}},$
(12)

for the $k$’th iteration. For a desired level of precision $\epsilon_{d}$,
once $\epsilon^{(k+1)}\leq\epsilon_{d}$ holds, we terminate the iterations. In
the last step of the algorithm, we truncate the elements with negligible
magnitude, e.g., those smaller than a certain threshold (for instance,
$5\times 10^{-5}$) in resulting controller $F$.

###### Remark 1

Our simulation results are obtained by incorporating update law (11) for the
first few iterations. Also, $W\in\mathbb{R}^{m\times n}$ is initially set to
the matrix of all ones.

###  IV-B Greedy Sparsification

Instead of application of the regularization in the second step, we may
alternatively consider the following greedy algorithm. We start with
$F^{(0)}=F_{o}$, and iteratively for $k=1,2,\dots$, the next controller
$F^{(k)}$ is derived by setting exactly one nonzero element of $F^{(k-1)}$
equal to zero. This nonzero element is chosen such that: the _distance_ of
updated constraint (10) from the boundary of positive-definiteness has a
minimal decrease. We define our proxy to distance as follows: Let us define
the matrix

$\displaystyle E^{(k)}:=\begin{bmatrix}\theta{R}&F^{(k)}-F_{o}\\\
(F^{(k)}-F_{o})^{T}&Z^{-1}\end{bmatrix},$ (13)

Also, we define $\Delta({i,j})\in{R^{m\times n}}$ to be a matrix of all zeros
expect with a $1$ at row $i$ and column $j$. Now, if we choose to push the
nonzero element $F_{ij}^{(k-1)}$ to $0$, then one inspects that matrix
$E^{(k)}$ is updated according to

$\displaystyle
E^{(k)}=E^{(k-1)}-F_{ij}^{(k-1)}\begin{bmatrix}0&\Delta({i,j})\\\
\Delta({i,j})^{T}&0\end{bmatrix}.$ (14)

Now, we consider the smallest eigenvalue of matrix $E^{(k)}$ as a measure of
distance to the boundary of constraint (10). This is equal to the inverse of
the largest eigenvalue of $(E^{(k)})^{-1}$. Hence, we will keep track of
matrix $(E^{(k)})^{-1}$. It is straight-forward to show that (14) is a rank-
two update to matrix $E^{(k-1)}$ according to

$\displaystyle E^{(k)}=E^{(k-1)}+V(i,j)H^{(k-1)}(i,j)V^{T}(i,j),$ (15)

where we have used the matrix of eigenvectors

$\displaystyle V(i,j):=[v_{1}(i,j)\,|\,v_{2}(i,j)],$ (16)

where $v_{1}(i,j)\in\mathbb{R}^{m+n}$ is a vector with $1/\sqrt{2}$ at
locations $i$ and $j+m$ and zero elsewhere. Also,
$v_{2}(i,j)\in\mathbb{R}^{m+n}$ a vector with $1/\sqrt{2}$ at location $i$,
and $-1/\sqrt{2}$ at location $j+m$, and zero elsewhere. Moreover,
$H^{(k-1)}(i,j)$ is given by

$\displaystyle
H^{(k-1)}(i,j):=\mathrm{diag}\left(-F_{ij}^{(k-1)},F_{ij}^{(k-1)}\right)$ (17)

We can see that $v_{1}$ is perpendicular to $v_{2}$. Let us derive the update
law for $(E^{(k)})^{-1}$ as well. The rank-two update derived for $E^{(k)}$
can be used in Woodbury matrix identity [28] to get

$\displaystyle(E^{(k)})^{-1}=(E^{(k-1)})^{-1}-(E^{(k-1)})^{-1}V(i,j)\times$
$\displaystyle\left(\left(H^{(k-1)}(i,j)\right)^{-1}+V(i,j)^{T}(E^{(k-1)})^{-1}V(i,j)\right)^{-1}\times$
$\displaystyle V(i,j)^{T}(E^{(k-1)})^{-1}.$ (18)

Then, we compute the maximum eigenvalue corresponding to each candidate for
update using the Power method [29]. Finally, at iteration $k$ we choose the
$i_{k}$ and $j_{k}$ using

$\displaystyle(i_{k},j_{k})\,\leftarrow\leavevmode\nobreak\
\underset{i,j}{\textrm{argmin}}\left\\{\lambda_{\max}\left((E^{(k)})^{-1}\right):\leavevmode\nobreak\
F^{(k-1)}_{ij}\neq 0\right\\}$ (19)

and set the element at row $i_{k}$ and column $j_{k}$, (i.e, element
$F_{i_{k}j_{k}}$) to zero. Using this method as the sparsification tool in
Algorithm $1$, the overall procedure would look like this:

Algorithm 3:
Greedy Sparsification inside Matrix Ellipsoid

  1. 1.

Start with $F^{(0)}=F_{o}$

  2. 2.

For $k=1,2,\dots$, choose $i$ and $j$ from

$(i_{k},j_{k})\,\leftarrow\leavevmode\nobreak\
\underset{i,j}{\textrm{argmin}}\left\\{\lambda_{\max}\left((E^{(k)})^{-1}\right):\leavevmode\nobreak\
F^{(k-1)}_{ij}\neq 0\right\\},$

and find $F^{(k)}$ with setting $F_{i_{k}j_{k}}^{(k-1)}$ to zero in
$F^{(k-1)}$. Stop when no such update exists.

##  V Numerical Examples

We examine the effectiveness of our proposed two-level algorithm. We define
the following measures that are useful for assessment of the output of the
sparsification. For a sparse controller $F$ that is built by searching around
the pre-designed $F_{o}$, we define the density level of $F$ relative to
$F_{o}$ (in percent) as

$\sigma_{d}:=100\times\dfrac{\|F\|_{0}}{\|F_{o}\|_{0}},$ (20)

and relative $\mathcal{H}_{\infty}$ performance loss (in percent) as

$\sigma_{p}:=100\times\frac{\|G_{F}(s)\|_{\mathcal{H}_{\infty}}-\|G_{F_{o}}(s)\|_{\mathcal{H}_{\infty}}}{\|G_{F_{o}}(s)\|_{\mathcal{H}_{\infty}}},$
(21)

where $G_{F}$ is the closed-loop transfer matrix from the disturbance to the
output upon application of controller $F$.

###  V-A Sample Feedback Control Designs

We consider two classes of control systems and examine the effectiveness of
our proposed procedure as follows:

####  V-A1 Randomly Generated System

we randomly create the system matrix $A\in\mathbb{R}^{n\times n}$ and input
matrix $B\in\mathbb{R}^{n\times m}$, where $n=m=30$ and each element of these
matrices is independently sampled from the standard normal distribution. We
set the rest of the matrices to be $C=I$, $D_{gv}=I$, $D_{gu}=I$, and
$B_{v}=B$. First, we use the weighted $\ell_{1}$ regularization method with
performance demand of $\gamma=2$ and design parameter $\theta=0.5$ and find a
value of feedback gain, which is denoted by $F_{s}\in\mathbb{R}^{m\times n}$.
On the other hand, we use the greedy method with the same value of $\theta$ to
find another feedback gain, namely, $F_{g}\in\mathbb{R}^{m\times n}$. The
density level and performance measure of these controllers are shown in Table
I. It turns out that sparsity level $\sigma_{d}$ in $F_{s}$ and $F_{g}$ have
decreased to about $65\%$ and $78\%$, while the relative performance
degradation compared to the design with controller $F_{o}$ is only $7.4\%$ and
$4.9\%$, respectively. The sparsity patterns of these feedback gain designs
have been illustrated in Fig. 1.

####  V-A2 Spatially-Decaying Interacting Subsystems:

We consider the matrix $A$ to be constructed as follows: We consider $n$
agents that are located in the square-shaped region $[0,1]\times[0,1]$, with a
random position $p_{l}\in[0,1]\times[0,1]$ for the $l$’th agent that was
uniformly picked. Then, we suppose that the element in row $i$ and column $j$
of $A$, namely, $A_{ij}$ is zero if $\|p_{i}-p_{j}\|_{2}>r$, and otherwise its
value is given by

$A_{ij}=c_{ij}\mathrm{e}^{-\alpha\|p_{i}-p_{j}\|_{2}^{\beta}},$

where ${c}_{ij}\sim\mathcal{N}(0,1)$, $\alpha$ determines the bandwidth of the
system, $\beta$ specifies the decay rate of the interaction between two
spatially located agents, and $r$ is the radius of the connectivity disk
around each agent. We consider $n=30$ with system parameters $\alpha=\beta=1$
and ${r=0.25}$. We set $B=I$, (i.e., each agent has an actuator that controls
its state), and $C=D_{gv}=D_{gu}=I$, $B_{v}=4B$. For the sparsification
parameters, we set $\gamma=5$, $\theta=0.5$, and similar to the previous
example, we obtain two sparse controllers $F_{s}$ and $F_{g}$ using the
weighted $\ell_{1}$ regularization and greedy method, respectively. Their
corresponding sparsity patterns are visualized in Fig. 2. According to Table
II, the full density level of $F_{o}$ has decreased to almost $37.5\%$ and
$36.3\%$, while the performance loss is around $1.2\%$ and $0.4\%$, in the
case of $F_{s}$ and $F_{g}$, respectively. In Fig. 3, we have illustrated the
links corresponding to connectivity architecture of the system as well as the
links that correspond to the information structure of sparse controller
$F_{s}$. Higher levels of achieved sparsity in this example are supposedly due
to the diagonal input matrix $B$.

Figure 1: The sparsity patterns of $F_{s}$ and $F_{g}$ for a uniform random system. Each blue dot corresponds to a nonzero element of the matrix. $F$ | $F_{o}$ | $F_{s}$ | $F_{g}$
---|---|---|---
$\|F\|_{0}$ | $900$ | $587$ | $701$
$\|G_{F}(s)\|_{\mathcal{H}_{\infty}}$ | $1.0001$ | $1.0740$ | $1.0487$
Table I: $\mathcal{H}_{\infty}$ norm and number of nonzero elements in $F_{o}$, $F_{s}$, and $F_{g}$ (for a randomly generated system). Figure 2: The sparsity patterns of feedback gains $F_{s}$ and $F_{g}$ for the exponentially decaying random system. Each blue dot corresponds to a nonzero element of the matrix. Figure 3: The links in the sparse controller $F_{s}$ and the system with spatially decaying interconnections. The links that exist in both structures have been highlighted as well.  $F$ | $F_{o}$ | $F_{s}$ | $F_{g}$
---|---|---|---
$\|F\|_{0}$ | $900$ | $337$ | $326$
$\|G_{F}(s)\|_{\mathcal{H}_{\infty}}$ | $3.4574$ | $3.4993$ | $3.4723$
Table II: $\mathcal{H}_{\infty}$ norm and number of nonzero elements in
$F_{o}$, $F_{s}$, and $F_{g}$ (for a sub-exponentially decaying randomly
generated system). Figure 4: Relative density levels and relative performance
losses of sparse controllers $F_{s}$ and $F_{g}$ vs the parameter $\theta$.
Recall that $F_{s}$ and $F_{g}$ are outputs of sparsification using Algorithms
2 and 3, respectively. Figure 5: The tradeoff curves between the density
levels and the relative performance degradation for sparse set of controllers
$F_{s}$ and $F_{g}$.

###  V-B Interplay between Resilience, Performance, and Sparsity

#### Randomly Generated System

Reconsidering the control system given in the previous subsection, we vary
parameter $\theta$ from $0$ to $1$ and look at the resulting density levels
and performance losses of the sparse controllers. In Fig. 4, we observe that
as the parameter $\theta$ increases, the density level of the outcome is
enhanced, while the performance measure is deteriorating. We observe that
these data represent a tradeoff between the relative density level of the
feedback gain and the relative performance loss upon sparsification, which is
further illustrated in Fig. 5. In these examples, the greedy method turns out
to be more conservative, while trading more performance for lower levels of
sparsity.

If we do the same experiment on the control system with spatially decaying
parameter, we arrive at similar trends and tradeoffs, which we omit for
brevity.

###  V-C Study of the Resilience

We create a random sample of the system with sub-exponentially decaying
interactions with $n=20$, $\alpha=1$, $\beta=0.5$, and $r=0.6$. Then, we study
the resilience of the designed sparse feedback control laws using re-weighted
$\ell_{1}$ regularization. We vary the parameter $\theta$ to find a sparse
feedback gain. Then, given the sparse controller $F_{s}$, we consider random
perturbations of each nonzero element of the controller that have the form

$[F_{s}]_{ij}\leftarrow[F_{s}]_{ij}+0.5\,\eta_{ij},$

where $\eta_{ij}\sim\mathcal{N}(0,1)$, independent of other elements. For each
design, (i.e., each value of the sparse controller corresponding to the value
of $\theta$), we find $5,000$ random perturbations. Then, for each randomly
perturbed feedback gain, we evaluate the $\mathcal{H}_{\infty}$ performance
measure of the closed-loop system.

The empirical cumulative density function (CDF) of the degradation in the
$\mathcal{H}_{\infty}$ performance measure (relative the sparse design) under
these random perturbations is illustrated in Fig. 6. We observe that for the
larger values of $\theta$, the sparse control design is more fragile.

Figure 6: The empirical CDF of the relative degradation in the
$\mathcal{H}_{\infty}$ performance measure $\sigma_{p}$ for 5,000 random
perturbations of the controller with different tuning parameter $\theta$.

##  VI Final Remarks

1) The general methodology followed in this paper can be summarized as
follows: first, obtain a region in the controller space inside which a
specific control objective holds. Then, explore the geometry of the region in
the controller space to find a sparse controller that inherits the performance
guarantee. While in this paper we have addressed $\mathcal{H}_{\infty}$
control design using the quadratic inequalities, as mentioned in [21], one can
do similar developments for the $\mathcal{H}_{2}$ norm performance measure as
well.

2) Because some of the matrix inequalities involved in the first method are
strict, solving the linear matrix inequalities in the first stage is subject
to a great deal of flexibility. For instance, from extensive numerical
simulations, we have learned that limiting the condition number of matrix $P$
in (7) is practical for finding a ”rich” enough region of controllers that
will be used in the sparsification.

3) Although we have limited the value of $\theta$ to $[0,1]$, some numerical
examples suggest that even for values of $\theta$ greater than $1$ the second
step for finding sparse controllers may result in satisfactory results. This
is due to the conservative nature of the quadratic regions of controllers
[21].

4) The exact running time analysis for the sparsification methods may not be
derivable. However, we expect that the greedy method becomes less
computationally expensive compared to $\ell_{1}$ regularization as the size of
the system increases.

5) Our extensive numerical studies suggest that as the designed controller
become sparser, the fragility of the design would increase as well.

6) The distance measure used for greedy sparsification, (i.e., smallest
eigenvalue) can be also changed to a number of other measures. For instance,
we can replace the protocol (19) to be

$(i_{k},j_{k})\,\leftarrow\leavevmode\nobreak\
\underset{i,j}{\textrm{argmin}}\left\\{\sum_{i=1}^{m+n}\lambda_{i}\left((E^{(k)})^{-1}\right):\leavevmode\nobreak\
F^{(k-1)}_{ij}\neq 0\right\\},$

where $\lambda_{i}$ denotes the $i$’th eigenvalue of the matrix.
Alternatively, one can set

$(i_{k},j_{k})\,\leftarrow\leavevmode\nobreak\
\underset{i,j}{\textrm{argmin}}\left\\{\log\left(\det\left((E^{(k)})^{-1}\right)\right):\leavevmode\nobreak\
F^{(k-1)}_{ij}\neq 0\right\\}.$

These measures capture alternative aspects of distance of the matrix from the
boundary of positive-definite matrices and we can derive simple update laws
for these measures using update formula (IV-B) (e.g. see [30]).

## References

  * [1] R. H. Takahashi, D. A. Dutra, R. M. Palhares, and P. L. Peres, “On robust non-fragile static state-feedback controller synthesis,” in _Proceedings of the 39th IEEE Conference on Decision and Control_ , vol. 5, 2000, pp. 4909–4914.
  * [2] D. Famularo, P. Dorato, C. T. Abdallah, W. M. Haddad, and A. Jadbabaie, “Robust non-fragile lq controllers: the static state feedback case,” _International Journal of control_ , vol. 73, no. 2, pp. 159–165, 2000.
  * [3] J. H. Park, “Robust non-fragile control for uncertain discrete-delay large-scale systems with a class of controller gain variations,” _Applied Mathematics and Computation_ , vol. 149, no. 1, pp. 147–164, 2004\.
  * [4] D. Peaucelle and D. Arzelier, “Ellipsoidal sets for resilient and robust static output-feedback,” _IEEE Transactions on Automatic Control_ , vol. 50, no. 6, pp. 899–904, 2005.
  * [5] M. Rotkowitz and S. Lall, “A characterization of convex problems in decentralized control,” _IEEE Transactions on Automatic Control_ , vol. 51, no. 2, pp. 274–286, 2006.
  * [6] N. Motee and A. Jadbabaie, “Optimal control of spatially distributed systems,” _IEEE Transactions on Automatic Control_ , vol. 53, no. 7, pp. 1616–1629, 2008.
  * [7] F. Lin, M. Fardad, and M. R. Jovanović, “Design of optimal sparse feedback gains via the alternating direction method of multipliers,” _IEEE Transactions on Automatic Control_ , vol. 58, no. 9, pp. 2426–2431, 2013.
  * [8] F. Dörfler, M. R. Jovanović, M. Chertkov, and F. Bullo, “Sparsity-promoting optimal wide-area control of power networks,” _IEEE Transactions on Power Systems_ , vol. 29, no. 5, pp. 2281–2291, 2014\.
  * [9] N. K. Dhingra and M. R. Jovanović, “A method of multipliers algorithm for sparsity-promoting optimal control,” in _2016 American Control Conference (ACC)_. IEEE, 2016, pp. 1942–1947.
  * [10] R. Arastoo, M. Bahavarnia, M. V. Kothare, and N. Motee, “Closed-loop feedback sparsification under parametric uncertainties,” in _IEEE 55th Conference on Decision and Control (CDC)_ , 2016, pp. 123–128.
  * [11] M. Bahavarnia, C. Somarakis, and N. Motee, “State feedback controller sparsification via a notion of non-fragility,” in _IEEE 56th Annual Conference on Decision and Control (CDC)_ , 2017, pp. 4205–4210.
  * [12] M. Bahavarnia, “State-feedback controller sparsification via quasi-norms,” in _American Control Conference (ACC)_. IEEE, 2019, pp. 748–753.
  * [13] S. Sojoudi, J. Lavaei, and A. G. Aghdam, _Structurally Constrained Controllers: Analysis and Synthesis_. Springer Science & Business Media, 2011.
  * [14] J. Lavaei, “Optimal decentralized control problem as a rank-constrained optimization,” in _2013 51st Annual Allerton Conference on Communication, Control, and Computing (Allerton)_. IEEE, 2013, pp. 39–45.
  * [15] M. S. Sadabadi and A. Karimi, “Fixed-structure sparse control of interconnected systems with polytopic uncertainty,” _IFAC Proceedings Volumes_ , vol. 47, no. 3, pp. 2588–2593, 2014.
  * [16] Y. Wang, J. Lopez, and M. Sznaier, “Sparse static output feedback controller design via convex optimization,” in _53rd IEEE Conference on Decision and Control_. IEEE, 2014, pp. 376–381.
  * [17] S. Fattahi and J. Lavaei, “On the convexity of optimal decentralized control problem and sparsity path,” in _2017 American Control Conference (ACC)_. IEEE, 2017, pp. 3359–3366.
  * [18] X. Wu and M. R. Jovanović, “Sparsity-promoting optimal control of systems with symmetries, consensus and synchronization networks,” _Systems & Control Letters_, vol. 103, pp. 1–8, 2017.
  * [19] M. R. Jovanović, P. J. Schmid, and J. W. Nichols, “Sparsity-promoting dynamic mode decomposition,” _Physics of Fluids_ , vol. 26, no. 2, p. 024103, 2014.
  * [20] C. Yu and M. Verhaegen, “Structured modeling and control of adaptive optics systems,” _IEEE Transactions on Control Systems Technology_ , vol. 26, no. 2, pp. 664–674, 2018.
  * [21] D. Peaucelle, D. Arzelier, and C. Farges, “Lmi results for resilient state-feedback with $\mathcal{H}_{\infty}$ performance,” in _43rd IEEE Conference on Decision and Control (CDC)_ , vol. 1, 2004, pp. 400–405.
  * [22] K. Zhou, J. C. Doyle, K. Glover _et al._ , _Robust and optimal control_ , vol. 40.
  * [23] M. Mital, A. K. Pani, S. Damodaran, and R. Ramesh, “Cloud based management and control system for smart communities: A practical case study,” _Computers in Industry_ , vol. 74, pp. 162–172, 2015.
  * [24] V. Blondel and J. N. Tsitsiklis, “Np-hardness of some linear control design problems,” _SIAM Journal on Control and Optimization_ , vol. 35, no. 6, pp. 2118–2127, 1997.
  * [25] M. Grant, S. Boyd, and Y. Ye, “Cvx: Matlab software for disciplined convex programming,” 2008.
  * [26] F. Zhang, _The Schur complement and its applications_. Springer Science & Business Media, 2006, vol. 4.
  * [27] E. Candes, M. Wakin, and S. Boyd, “Enhancing sparsity by reweighted $\ell_{1}$ minimization,” _Journal of Fourier Analysis and Applications_ , vol. 14, pp. 877–905, July 2008.
  * [28] M. A. Woodbury, “Inverting modified matrices,” _Memorandum report_ , vol. 42, no. 106, p. 336, 1950.
  * [29] Y. Saad, _Numerical methods for large eigenvalue problems: revised edition_. Siam, 2011, vol. 66.
  * [30] H. K. Mousavi, Q. Sun, and N. Motee, “Space-time sampling for network observability,” _arXiv preprint arXiv:1811.01303_ , 2018.


# Effects of group size and noise on cooperation in population evolution of
dynamic groups

Hong-Bin Zhang<EMAIL_ADDRESS>Deng-Ping Tang  School of Computer Science
and Technology, Soochow University, Suzhou, 215006, China

###### Abstract

In a large population, the agents temporally form one by one group of the
Public Goods Game (PGG), and size of one group is randomly distributed at
$g\in[g_{l},g_{u}]$. Players in it have two strategies to be chosen to
cooperate, or to defect for playing the PGG. Based on this structure we
investigate the evolution of cooperation in public goods games as a function
of the noise level underlying strategy adoptions with the group size
distribution. In the process of the evolution, the payoffs of
cooperators(defectors) have related to the strategy selection implemented by
the death-birth process. We consider a particular case of time scales $s$
between the interactions of the agents and the strategy update, and if $s$=1,
the time scale of one strategy update is same to that of one interaction of
agents engaging in a PGG, where the cooperative behavior can be analytically
solved. We find the cooperative effect is affected by the noise level and the
group distribution, consequently it induces different dynamic behavior. In
this model, when at the little value of the noise ($K=0.1$), if the system has
different initial concentrations of cooperation, a mixed state with
coexistence of cooperators and defectors, and AllC state are observed. While
at other values of the noise level, the critical parameter of the
multiplication factor $r$ is obtained, and above the critical value of $r_{c}$
it induces a totally cooperative behavior, which is fitted well with
analytical expressions. The larger uncertainty with the increase of the noise
level is, the more he hesitate and unsure how to choose which strategy for a
medium value of the multiplication factor $r$. So, with a greater value of
noise strength, a larger value of $r$ induces a cooperative behavior; however
a small value of $r$ leads to the mixed states of cooperators and defectors,
and between the two players probably randomly choose to cooperate or to defect
in a sphere of unstable solution, where a hysteresis effect occurs for a
bigger size distribution of group and greater uncertainty.

##  1 Introduction

The rule of natural selection that is proposed by Darwin favors the selfish
individuals who obtain the maximum benefit at the expense of the others[1, 2,
3, 4]. For example, in a competing population, the individuals are divided
into the cooperators and defectors, those cooperators are willing to
contribute money so as to get benefits for the whole population while the
defectors (also called free riders) don’t pay anything but share the
collections with the formers. Obviously, to be the defectors are rational and
they will be winners since they get more benefit than the cooperators.
Therefore, the defective behavior will be favored[3, 5]. So the selection of
cooperation in nature is not an obvious behavior. However, the phenomenon of
cooperation is ubiquitous in nature and society.

The evolutionary game theory is one such powerful tool and is widely used to
discover or understand the competition behavior in nature[6]. In order to get
a complete understanding of local interactions of multi-player systems, the
Public Goods game (PGG) has been widely used as a paradigm [7]. In the PGG,
$g$ players forming a group can take one of the two strategies to play the
game. One strategy is to cooperate ($C$) and the other is to defect ($D$). A
cooperator pays a cost $c$ so as to get a benefit $rc$ for the group, where
$r$ is the multiplication factor. Assume that $n$ cooperators are in the
group, then a total gained benefit $nrc$ will be shared by all the players in
the group. Therefore, the defectors have higher net payoffs since they do not
pay extra costs. When formed groups of $g$ players enter into the PGG in a
large population, for $1<r<g$ [8], the situation of the cooperators are worse
off, and a defective population will evolve to a dominate position.

The network effect that the competition environment of the population has been
localized will lead to the emergence of cooperation among selfish
individuals[2, 9, 10, 11, 12, 15, 16]. Szolnoki proposed a model that
individuals participate a PGG on regular graphs of different group size, which
in turn introduce effective overlapping triangles to the interaction structure
and study the impact of this fact on the noise $K$ dependence of the evolution
of cooperation in spatial public goods games[17]. Stivala, Zhang and Chen
assumed that the players have not perfect rationality, embedding the noise in
the evolution of cooperation can lead to the survival of cooperators[18, 19,
20]. Players selected from a well-mixed population randomly forming group
engage in the PGG of constant group size, the defectors getting better payoffs
on average[21]. Hauert et al. considered a coevolving model with population
dynamics and game dynamics [8]. The population dynamics involves the death-
birth processes of cooperators and defectors besides the empty space due to
the death of individuals. One empty space has been filled with one cooperator
(defector) depending on the probability of reproduction of one cooperator
(defector). The game dynamics involves the interactions of dynamic formed
group of $g$ size and thereafter causes the varied average payoffs of
cooperators and defectors in the population. The average payoffs will in turn
determine the fractions of cooperators and defectors together with the
fraction of empty space. The temporal treatment of interactions between
individuals will also greatly affect the cooperation of the whole population.
The parameter $s$ is considered by Roca to reflect the difference between the
evolutionary time scales of the interactions among individuals and the
strategy selection implemented by the death-birth process [25]. In finite
populations, the time scale parameter $s$ significantly affects the fixation
probability of totally cooperative population when the interaction of randomly
pairing individuals is considered. When the individuals are temporally formed
to play the $N$-person snowdrift game[26] with a group size distributions, a
totally cooperative population could occur under the strategy selection
implemented by the imitation process [27]. Furthermore, when the individuals
could choose a time slot to play the prisoner’s dilemma game [1], the
cooperative behavior of the whole population with spatial structure will be
affected by the various time slots [28, 29].

In what follows, we will investigate the effect of formed group of $g$ size
with the noise $K$ dependence on the evolution of cooperation by a fixed time
slot $s=1$. We considered that the sizes of the formed groups are varied but
the total number of players in the population remains unchanged. The
probability of the offspring players inheriting the strategies of their
parents. This reproduction probability is characterized by the Fermi function
with a noise effect. The death-birth process is applied in our model. It is
related to the total payoffs of cooperators(defectors). Compared with previous
essays when the noise $K$ is introduced, in addition to the simulation results
we can obtain, an analytical method related to the noise $K$ is also provided.
The details are as follows.

##  2 Model

We consider a large population of $M$ players. Some players in the population
will be randomly selected to form size-distributed groups so as to play the
PGG from time to time. For example, investors hunt for yield by engaging into
different competing groups. Unlike the structured population with fixed local
opponents, the dynamic interaction between agents will happen in the whole
system. In the PGG, the cooperators and the defectors will obtain their
payoffs respectively. Suppose that there are $g$ agents and $n$ cooperators in
the group playing the PGG, the total payoff to the cooperators in the group is
given by

$V_{c}(g,n)=(\frac{{rn}}{g}-1)n,$ (1)

Meanwhile, the total payoff to the defectors in the group is given by:

$V_{d}(g,n)=\frac{{rn}}{g}(g-n),$ (2)

Where $r$ ($>1$) is the multiplication factor. With simplicity and without
loss of generality, the cost that a cooperator contributes is taken as $c=1$.
In this paper, the payoffs that the two type agents accumulated after they
play the games will depend on the selection process the death-birth process to
simulate the strategy selection. The payoffs that the two type agents
accumulated decide the reproduction probabilities of their offspring. The
detailed evolutionary death-birth process is as follows. In an evolutionary
time step, there are two processes: at the first one, the initial payoffs of
cooperators and defectors are preseted as $V_{c}$=$V_{d}$=0, and then randomly
$g$ agents are selected from $M$ players to play the PGG. After the game, the
payoffs that the two type agents (cooperators and defectors) are denoted with
$V_{c}$ and $V_{d}$ respectively. The sizes $g$ of the dynamically formed
groups follow a random distribution with $g_{l}$ and $g_{u}$ being the lower
and upper limits of the group sizes. The second process is the death-birth
process. A randomly selected agent from the whole population is removed and
the two-type agents compete to birth an offspring who inherits his father’s
strategy. The reproduction probability $w_{c}$ of the cooperators is given by
$w_{c}=1/[1+e^{(\tilde{V}_{d}-\tilde{V}_{c})/K}]$ where $K$ is the noise
value. The learning mechanism helps the offsprings with a higher payoff in the
evolution of population. The reproduction probability of the defectors is
$w_{d}$=1-$w_{c}$. Then randomly make a agent die, and give birth to a new
offspring with certain probabilities. The C(D) feature of a new offspring
subjects to the reproduction probability $w_{c}$($w_{d}$). So far, after the
death-birth process, one evolutionary time step is over.

In this paper, $M=10^{5}$ players are selected to proceed the Monte-Carlo
simulation in which one Monte-Carlo step (MCS) consists of $10^{5}$ time
steps. When the system is in equilibrium state, we get the average cooperative
frequency $f_{c}$ by taking an average of $10^{4}$ MCS.

##  3 Theory and Results

The frequency of cooperation $f_{c}$ is $M_{c}/M$, where $M_{c}$ is the number
of cooperators at time step $t$. We define that $q(g)$ be the size of the
dynamically formed groups following random distribution. It should satisfies
$\sum\nolimits_{g=g_{l}}^{g_{u}}{q(g)}=1$ with $g\in[g_{l},g_{u}]$. During the
interactions among agents, suppose that the probability for the randomly
formed group in the $i$ round having group size $g_{i}$ and $n_{i}$
cooperators in it is given by

$p(g_{i},n_{i})=q(g_{i})\left({\begin{array}[]{*{20}c}{g_{i}}\\\ {n_{i}}\\\
\end{array}}\right)f_{c}^{n_{i}}(1-f_{c})^{g_{i}-n_{i}},$ (3)

where the first factors $q(g_{i})$ of Eq. (3) from the right side indicates
the probability of having a dynamic formed group of size $g_{i}$. The
probability of formed group of $g_{i}$ agents involving $n_{i}$ cooperators is
represented with the binomial distribution approximately. It shows that the
probability of having $n_{i}$ cooperators in the formed group of randomly
selected $g_{i}$ agents depends: (1) linearly on the probability of having
$n_{i}$ cooperators and $g_{i}-n_{i}$ defectors, i.e, $f_{c}^{n_{i}}$
$(1-f_{c})^{g_{i}-n_{i}}$. (2) on multiple possible expression of
$f_{c}^{n_{i}}$ $(1-f_{c})^{g_{i}-n_{i}}$, we pick up all non-repetitive
processes denoted as a binomial coefficient
$C_{g_{i}}^{n_{i}}$=$\tbinom{g_{i}}{n_{i}}$. After one round game of PGG, the
accumulated profits of C and D agents is given by:

$V_{tot}(g_{i},n_{i})=V_{c}(g_{i},n_{i})+V_{d}(g_{i},n_{i})=(r-1)n_{i}.$ (4)

In general, the dynamic equation can be obtained:

$\frac{{df_{c}}}{{dt}}=\frac{1}{M}(\sum\limits_{g=g_{l}}^{g_{u}}{\sum\limits_{n=0}^{g}{p(g,n)w_{c}(g,n)}}-f_{c}).$
(5)

We will set the right side equal to zero, a static solution of $f_{c}$ can be
obtained. If $g_{u}$ ($g_{u}<5$) has a relatively small value, analytical
solution may be obtained. For example, when $g\in[1,3]$ and $q(g)$ follows an
uniform distribution, the static solution of Eq.(5) can be obtained by solving
the following univariate cubic equation:

$af_{c}^{3}+bf_{c}^{2}+cf_{c}+\frac{1}{2}=0.$ (6)

where

$a=-\frac{1}{6}+\frac{1}{3(1+e^{\frac{3-3r}{K}})}+\frac{1}{1+e^{\frac{6-2r}{3K}}}+\frac{1}{1+e^{\frac{3+r}{3K}}}.$
(7)
$b=\frac{1}{3}(2-\frac{2}{1+e^{\frac{1}{K}}}+\frac{1}{1+e^{\frac{2-2r}{K}}}+\frac{3}{1+e^{\frac{6-2r}{3K}}}-\frac{6}{1+e^{\frac{3+r}{3K}}}).$
(8)
$c=-2+\frac{2}{3(1+e^{\frac{1}{K}})}+\frac{1}{3(1+e^{\frac{1-r}{K}})}+\frac{1}{1+e^{\frac{3+r}{3K}}}.$
(9)

As we can see that the theoretical curves fit well to all the simulation data
(see Figure 1)

Figure 1: (color online) The dependence of the frequency of cooperation
$f_{c}$ on the multiplication factor $r$ at different noise strength $K$. The
size of formed group $g\in[1,3]$ is set with the initial frequency of
cooperation $f_{c,i}=0.5$. The symbols are obtained by simulations and the
curves are calculated by theory.

The noise level $K$ will affect the frequency of cooperation. In Figure 1, we
show the frequency of cooperation as a function of $K$ when $g\in[1,3]$. For a
small value of $r$, we observe that a totally cooperative state with a smaller
value of $K$. In previous works, it showed the truth that cooperation was more
difficult to sustain for lower enhancement factors, e.g.,$1<r<g$ [8], on the
other hand, we argue that a lower value of the noise with a smaller size
distribution of dynamic formed group facilitates the domination of cooperators
based on payoffs of all the whole group. For a given value of $K$=0.1, when
the size distribution of formed group $g$ has been enlarged, e.g.,for
$g_{u}>5$, it benefits the the survival of cooperators(see Figure 2).
Moreover, we observe the non-trivial result that for different initial
concentrations of cooperation when $g\in[1,5]$, $r=3.5$ and $K$=0.1, it takes
$10^{5}$ time steps to evolve to the mixed state of C and D, or to the AllC
state.

Figure 2: (color online) After $20$ Monte-Carlo steps, the system converge to
one stable mixed state of C and D, or AllC state when $g\in[1,5]$, $r=3.5$ and
$K$=0.1. The symbols are obtained by simulations and the curves are calculated
by theory.

When the value of $K$ is increased, e.g., for $K=1.0$ and $g\in[1,5]$, it
shows that hysteresis phenomenon occurs in Figure 3 when the frequency of
cooperation changes with the multiplication factor $r$. Eventually, the
simulation has evolved to two stable solutions, and the cut-off line reacts
from one stable solution to another. We explore an explanation that one
defector(cooperator)gradually invade the territory of opposite side on two
stable branches, which are depicted with solid lines. What we are depicting
two branches, one is an upwards branch, the other is a downwards branch.
Between two branches, there is an unstable solution by a dashed line. In
Figure 3, the upwards branch reveals that the frequency of cooperation goes up
with the increase of $r$, as the other goes down with the decrease of $r$. In
the process evolving to upwards, or to downwards branch, what we set the
current initial frequency of cooperation depends on $f_{c}$ in the equilibrium
state with the value of $r$ in the previous process. In the upwards branch,
e.g., for the minimum value of $r$=1.1, the initial value of $f_{c}$ is close
to 0.00001. In the another branch, for a bigger value of $r$=9.0, the initial
value of $f_{c}$ is set to 0.99999.

Figure 3: (color online) The dependence of the frequency of cooperation on the
multiplication factor at different initial frequency of cooperation. We take a
random size distribution with $g\in[1,5]$ and $K=1.0$.

In Figure 4, for a smaller value of $K$, the system is in a cutoff period.
Here, we investigate how one C-agent invades the area of all defectors by the
description of square dots moving upward (upwards branch), and how one D-agent
bursts into the territory of all cooperators from another side by those
triangulation dots moving downward (downwards branch). By a curve passing
through the triangulation dots, or square dots, the theoretical solution is
obtained.

Figure 4: (color online) When $K=0.1$ and $g\in[1,5]$, we can obtain results
similar to Figure 3. The symbols describe the process that
(upwards(downwards)) branches are shaped by one C(D)-agent invading the
territory of Defectors(Cooperators).

More detailed explantation we make can be seen in the Figure 5. It gives a
phase diagram of a single D-agent invading the region of AllC (see triangular
point of Figure 4). We divide the phase diagram into three regions, e.g.,
region (I),region (II) and region(III). For a lower value of $K$ and $r$, an
AllC population can be certainly invaded by one D-agent in the region (I). The
region (II) is a part that one D-agent has chances to invade C-agents
population, and also for a C-agent to invade D-agents population. However, the
region(III) is that AllC population cannot be invaded by one D-agent. In the
region(III), if the system is in a state where there are one defector and
$M$-1 cooperators, and when the probability of increasing one cooperator is
higher than that of a defector, the last defector will be extinct. At the
moment, when the size distribution $g\in[g_{l},g_{u}]$ and $K$ will be given,
the outcome of a formed group will have $0$ or $1$ defector. The probability
$p_{1d}$ for the group of size $g$ having one defector is given by
$p_{1d}=g/M$, and the probability $p_{0d}$ for having zero defector is given
by $p_{0d}=1-g/M$. The probability $\beta_{c}$ for increasing one cooperator
is given by

$\displaystyle\beta_{c}$
$\displaystyle=\sum\limits_{g_{j}=1}^{g_{u}}\frac{1}{g_{u}}*{\left(\left(1-\frac{g_{j}}{M}\right)*{{\frac{1}{1+e^{\frac{0-\left(\frac{r*g_{j}}{g_{j}}-1\right)*g_{j}}{K}}}+\frac{g_{j}}{M}*\frac{1}{1+e^{\frac{\frac{r*(g_{j}-1)}{g_{j}}*1-\left(\frac{r*(g_{j}-1)}{g_{j}}-1\right)*(g_{j}-1)}{K}}}}}\right)}$
(10) $\displaystyle-\frac{M-1}{M},$

The above equation is obtained by directly expanding the right-hand-side of
Eq. (5) and multiply $M$. Accordingly, the the probability $\beta_{d}$ for
increasing one defector is given by

$\displaystyle\beta_{d}$
$\displaystyle=\sum\limits_{g_{j}=1}^{g_{u}}\frac{1}{g_{u}}*{\left(\left(1-\frac{g_{j}}{M}\right)*{{\frac{1}{1+e^{\frac{\left(\frac{r*g_{j}}{g_{j}}-1\right)*g_{j}-0}{K}}}+\frac{g_{j}}{M}*\frac{1}{1+e^{\frac{\left(\frac{r*(g_{j}-1)}{g_{j}}-1\right)*(g_{j}-1)-\frac{r*(g_{j}-1)}{g_{j}}*1}{K}}}}}\right)}$
(11) $\displaystyle-\frac{1}{M},$

Compared Eq. (10) with Eq. (11), we get $\beta_{c}+\beta_{d}$=0. Obviously
when $\beta_{c}>0$, $\beta_{c}>\beta_{d}$ and the AllC state will happen. For
$g\in[1,5]$, a more simpler form is given as follows when $r>1$ and $K=0.1$,

$\beta_{c}=-\beta_{d}=\frac{1}{5}(Q1+Q2+Q3+Q4+Q5)-\frac{1}{M}$ (12)

where

$Q_{1}=\frac{1}{2M}$ (13)

$Q_{2}=\frac{2}{M}\frac{1}{1+e^{\frac{-1}{K}}}$ (14)

$Q_{3}=\frac{3}{M}\frac{1}{1+e^{\frac{(\frac{2}{3}*r)-2}{K}}}$ (15)

$Q_{4}=\frac{4}{M}\frac{1}{1+e^{\frac{(\frac{3}{2}*r)-3}{K}}}$ (16)

$Q_{5}=\frac{5}{M}\frac{1}{1+e^{\frac{(\frac{12}{5}*r)-4}{K}}}$ (17)

Based on Eqs. (12)-(17), the critical value of $r_{invade}$ for invading the
AllC population can be obtained by setting $\beta_{c}=0$. In the simulations,
we take $M=10^{5}$ and $K=0.1$, then the theoretical calculation of
$r_{invade}$ is given by 2.75, which is close to that obtained in the
simulations.

Figure 5: (color online) Phase diagram of a single D-agent invading the region
of AllC. The phase indicates three states: high probability of invasion,
moderate probability of invasion, and impossible invasion with the description
of region (I),region (II) and region(III) respectively. All points can be
obtained when $M=10^{5}$, $K=0.1$ and $g\in[1,5]$ by going through average 107
MCS steps.

##  4 Summary

Summarizing, in a large population, we investigate the evolution of
cooperation in the public goods game on dynamic randomly formed group with the
group size $g$ and the noise level $K$. The payoffs of cooperators and
defectors depends on the strategy update in the death-birth process. We show
that varied noise strength $K$ at $g\in[1,3]$ induce different dynamic
behavior of cooperation. At $K=0.5$ or $K=0.1$ when $g\in[1,3]$, approaching
to totally cooperative state has a lower critical value of
$r_{c}$($r_{c}=3.8$, or $r_{c}=3.0$, see Figure 1) compared with that of
$K=1.0$. For a greater value of the noise $K$ ($K$=1.0) and $g\in[1,3]$, ALLC
cooperative population also emerges when the multiplication factor is above a
critical value $r_{c}$($r_{c}=5.6$). Moreover, for $K$=0.1, $g\in[1,5]$ and
$r=3.5$, starting with different initial concentrations of cooperation, AllC
state and mixed state of C and D can be observed. For a greater value of the
noise $K=1.0$, we study the evolution of cooperation dependence on different
multiplication factor $r$, it can be observed that there exist a sphere of
unstable solutions between the two stable extremes. The phenomenon of
hysteresis effect can be given a comprehensible explanation for one C(D)-agent
to invade the territory of the D(C) population with the change of the
multiplication factor $r$. For $K$=1.0, on the downwards branch, a one by one
D-agent invade the region of C-population with decrease of the multiplication
factor $r$, and another phase of mixed state of C and D is eventually reached;
however on the upwards branch, similarly one by one C-agent invade the sphere
of D-population with the increase of $r$, finally approaching to the AllC
state. With a greater value of the noise strength and the mediate values of
multiplication factor $r$, one agent is in a dilemma about whether to select
the strategies of cooperation or defection, thus randomly selection leads to
unstable solutions(see Figure 3). Furthermore, for $K=0.1$ and $g\in[1,5]$,
the system enters into a cutoff period, and for a small value of $r$, one
defective agent easily invades the sphere of cooperators; but when above the
critical value of $r_{invade}$, e.g, for $r_{invade}=2.75$, he is hard to
invade the territory of cooperators. With a greater noise level $K$, e.g.,
$K=1.0$, the uncertainty makes the value of $r_{invade}$ more larger. Also, we
find that the stationary frequency of cooperation in this model could be
solved analytically and fits well with the simulated results. All the
analytical results are verified by computer simulations.

## References

## References

  * [1] R. Axelrod, W. D. Hamilton, Science 211 (1981) 1390.
  * [2] J. Maynard Smith, Evolution and the Theory of Games, Cambridge Univ. Press, Cambridge, 1982.
  * [3] J. Hofbauer, K. Sigmund, Evolutionary Games and Population Dynamics, Cambridge Univ. Press, 1998.
  * [4] M. A. Nowak, Nature 359 (1992) 6398.
  * [5] G. Szabo, C. Hauert, Phys.Rev.Lett 89 (2002) 11.
  * [6] M. A. Nowak, K.Sigmund, Nature, 437 (2005)7063
  * [7] F.C. Santos, M.D. Santos,and J.M. Pacheco, Nature 454 (2008)7201.
  * [8] C. Hauert, M. Holmes, M. Doebeli, Proc. R. Soc. B 273 (2006) 2565. Evolutionary Dynamics, Harvard Univ. Press, 2006.
  * [9] W. D. Hamilton, J. Theor. Biol. 7 (1964) 1.
  * [10] D. S. Wilson, E. Sober, Behav. Brain Sci. 17 (1994) 585.
  * [11] R. Trivers, Q. Rev. Biol. 46 (1971) 35.
  * [12] M. A. Nowak, K. Sigmund, 393 (1998) 573.
  * [13] M. A. Nowak, R. M. May, Nature 359 (1992) 826.
  * [14] A. Traulsen, M. A. Nowak, Proc. Natl. Acad. Sci. USA 103 (2006) 10952.
  * [15] C. Hauert, M. Doebeli, Nature 428 (2004) 643.
  * [16] C. Hauert, S. De Monte, J. Hofbauer, K. Sigmund, Science 296 (2002) 1129.
  * [17] A. Szolnoki, M. Perc, and G. Szabó,Phys. Rev. E 80(2009)056109
  * [18] A.Stivala, Y.Kashima,and M.Kirley, Phys. Rev. E, 94 (2016) 032303.
  * [19] G. Q. Zhang, Q. B. Sun, and L. Wang, Chaos, Solitons & Fractals, 51 (2013) 7.
  * [20] X. J. Chen, L. Wang, EPL, 90 (2010) 3.
  * [21] X.Y Xu, H.B Zhang Physica A 625(2023) 129005
  * [22] A. Szolnoki, M. Perc, Phys. Rev. E 81 (2010) 057101.
  * [23] M. A. Janssen, R. L. Goldstone, J. Theor. Biol. 243 (2006) 134.
  * [24] X. Wang, R. Ding, J. Zhao, C. Gu, Appl. Math. Comp. 413 (2022) 126624.
  * [25] C. P. Roca, J. A. Cuesta, A. Sánchez, Phys. Rev. Lett. 97 (2006) 158701.
  * [26] D. F. Zheng, H. P. Yin, C. H. Chan, and P. M. Hui, EPL 80 (2007) 18002.
  * [27] C. Xu, P. M. Hui, Phys. Rev. E 105 (2022) 054309.
  * [28] T. Johnson, O. Smirnov, Communication Biology 4 (2021) 1283.
  * [29] T. Johnson, Chaos, Solitons & Fractals: X, 9 (2022) 100081.