1

Local pinning of networks of multi-agent systems
with transmission and pinning delays

arXiv:1604.00459v1 [cs.SY] 2 Apr 2016

Wenlian Lu Senior Member, IEEE Fatihcan M. Atay

Abstract—We study the stability of networks of multi-agent
systems with local pinning strategies and two types of time
delays, namely the transmission delay in the network and
the pinning delay of the controllers. Sufficient conditions for
stability are derived under specific scenarios by computing
or estimating the dominant eigenvalue of the characteristic
equation. In addition, controlling the network by pinning a single
node is studied. Moreover, perturbation methods are employed
to derive conditions in the limit of small and large pinning
strengths. Numerical algorithms are proposed to verify stability,
and simulation examples are presented to confirm the efficiency
of analytic results.

I. I NTRODUCTION
Control problems in multi-agent systems have been attracting attention in diverse contexts [1]–[7]. In the consensus problem, for example, the objective is to make all agents converge
to some common state by designing proper algorithms [2]-[5],
such as the linear consensus protocol
ẋi = −

n
X

Lij xj (t), i = 1, . . . , n.

(1)

j=1

Here, xi ∈ R is the state of agent i and Lij are the components
of the Laplacian
P matrix L, satisfying Lij ≤ 0 for all i 6= j
and Lii = − j6=i Lij . The Laplacian is associated with the
underlying graph G, whose links can be directed and weighted.
It can be shown that, if the underlying graph has a spanning
tree, then all agents converge to a common number, which
depends on the initial values [1], [4], [5]. On the other hand,
if it is desired to steer the system to a prescribed consensus
value, auxiliary control strategies are necessary. Among these,
pinning control is particularly attractive because it is easily
realizable by controlling only a few agents, driving them to
the desired value s through feedback action:
ẋi = −

n
X
j=1

Lij xj (t) − δD (i)c(xi − s),

i = 1, . . . , n, (2)

where D denotes the subset of agents where feedback is
applied, with cardinality |D| = m, δD (i) is the indicator
W. L. Lu is with the Centre for Computational Systems Biology and School
of Mathematical Sciences, Fudan University, Shanghai, P. R. China (email:
w.l.lu@ieee.org); F. M. Atay is with the Max Planck Institute for Mathematics
in the Sciences, Leipzig, Germany (email: fatay@mis.mpg.de).
The authors thank the ZiF (Center for Interdisciplinary Research) of
Bielefeld University, where part of this research was conducted under the
program Discrete and Continuous Models in the Theory of Networks. FMA
acknowledges the support of the European Union (FP7/2007-2013) under
grant #318723 (MatheMACS). WLL acknowledges the support of the National
Natural Sciences Foundation of China under Grant (61273309) and the
Program for New Century Excellent Talents in University (NCET-13-0139).

function (1 if i ∈ D and 0 otherwise), and c > 0 is the
pinning strength. Eq. (2) provides the local strategy that pins a
few nodes to stabilize the whole network at a common desired
value. The following hypothesis is natural in pinning problems
and assumed in this paper.
(H) Each strongly connected component of G without
incoming links from the outside has at least one node in D.
The following result is proved in [8], [9].
Proposition 1: If (H) holds, then system (2) is asymptotically stable at xi = s ∀i.
In many networked systems, however, time delays inevitably
occur due to limited information transmission speed; so Proposition 1 does not apply. In this paper we consider systems with
both transmission and pinning delays,

ẋi = −

n
X

Lij (xj (t−τr )−xi (t))−cδD (i)(xi (t−τp )−s),

j=1,j6=i

(3)
for i = 1, . . . , n, where τr denotes the transmission delay
in the network and τp is the pinning delay of the controllers.
Several recent papers have addressed the stability of consensus
systems with various delays. It has been shown that consensus
can be achieved under transmission delays if the graph has
a spanning tree [13]-[15]. However, if a sufficiently large
delay is present also in the self-feedback of the node’s own
state, then consensus may be destroyed [16]; similar conclusions also hold in cases of time-varying topologies [17]–[19]
and heterogeneous delays [20]-[22]. The stability of pinning
networks with nonlinear node dynamics have been studied
in [6]–[12], [23]–[26]. However, the role of pinning delay
was considered in only a few papers [23]–[26], where it was
argued that stability can be guaranteed if the pinning delays are
sufficiently small. Precise conditions on the pinning delay for
stability, the relation to the network topology, and the selection
of pinned nodes have not yet been addressed.
In this paper, we study the stability of the model (3)
under both transmission and pinning delays. First, we derive
an estimate of the largest admissible pinning delay. Next,
we consider several specific scenarios and present numerical
algorithms to verify stability by calculating the dominant
eigenvalue of the system. Included among the scenarios are
the cases when only a single node is pinned in the absence
of transmission delay, or when the transmission and pinning
delays are identical. Finally, we use a perturbation approach
to estimate the dominant eigenvalue for very small and very
large pinning strengths.

2

II. N OTATION

AND

P RELIMINARIES

A directed graph G = {V, E} consists of a node set V =
{v1 , . . . , vn } and a link set E ⊆ V × V. A (directed) path
of length l from node vj to vi , denoted (vr1 , . . . , vrl+1 ), is a
sequence of l+1 distinct vertices with vr1 = vi and vrl+1 = vj
such that (vrk , vrk+1 ) ∈ E for k = 1, . . . , l. The graph is called
strongly connected if there is a directed path from any node
to any other node, and it is said to have a spanning tree if
there is a node vp ∈ V such that for any other node j there is
a path from vp to vj .
We denote the imaginary unit by j and the n × n identity
matrix by In . For a matrix L, Lij denotes its (i, j)th element
and L⊤ its transpose. The Laplacian matrix L is associated
with the graph G in the sense that there is a link from vj to vi
in G if and only if Lij 6= 0. We denote the eigenvalues of L by
{θ1 , . . . , θn }. Recall that zero is always an eigenvalue, with the
corresponding eigenvector [1, . . . , 1]⊤ , and Re(θi ) > 0 for all
nonzero eigenvalues θi . Furthermore, if the graph G is strongly
connected (or equivalently, if L is irreducible), then zero is
a simple eigenvalue of L. The diagonal element Lii is the
weighted in-degree of node i. Let K = diag{L11 , . . . , Lnn }
be the diagonal matrix of in-degrees and A = K − L. Let
yi = xi − s, y = [y1 , . . . , yn ]⊤ , and D = diag{d1 , · · · , dn }
with di = δD (i). System (3) can be rewritten as
ẏ = −Ky + Ay(t − τr ) − cDy(t − τp ).

(4)

Considering solutions in the form y(t) = exp(λt)ξ with λ ∈ C
and ξ ∈ Cn , the characteristic equation of (4) is obtained as

characteristic root λ∗ of (5) such that Re(λ∗ ) ≥ 0. Applying
the Gershgorin disc theorem to (5), we have
X
X
|λ∗ + Lii + cdi | ≤
|Lij || exp(−λ∗ τr )| ≤
|Lij | = Lii
j6=i

j6=i

(8)

for some i, which implies
[Re(λ∗ ) + Lii + cdi ]2 + [Im(λ∗ )]2 ≤ L2ii .
Since Lii , c, di ≥ 0, it must be the case that Re(λ∗ ) =
Im(λ∗ ) = 0; i.e., λ∗ = 0. Then exp(−τr λ∗ ) = 1, and
since τp = 0, (5) gives det(λ∗ In + L + cD) = 0. This,
however, contradicts Corollary 1. Therefore, when τp = 0,
all characteristic roots of (5) have negative real parts.
We now let τp ≥ 0. Suppose (5) has a purely imaginary
root λ = jω, ω ∈ R. By (8), we have, for some index q,
X
|jω + Lqq + cdq exp(−jωτp )| ≤
|Lqj || exp(−jωτr )
j6=q

=

X
j6=q

|Lqj | = Lqq

implying
q
[Lqq + cdq cos(ωτp )]2 + [ω − cdq sin(ωτp )]2 ≤ Lqq .
Thus,

(cdq )2 + ω 2 + 2cdq (Lqq cos(ωτp ) − ω sin(ωτp )) ≤ 0. (9)

We claim that q must be a pinned node. For if dq = 0, then
ω must be zero, which implies that zero is a characteristic
χ(λ) := det [λIn + K − A exp(−λτr ) + cD exp(−λτp )] = 0. root of (5), contradicting Corollary 1. Therefore d = 1. In
q
(5) the notation of (6), the inequality (9) can then be written
The asymptotic stability of (4) is equivalent to all characteristic as F (w, c, L , τ ) ≤ 0. By (7), however, we have that
qq p
roots λ of (5) having negative real parts. The root having the F (w, c, L , τ ) > 0 for all p ∈ D, ω ∈ R and τ < τ ∗ .
qq p
p
p
largest real part will be termed as the dominant root or the We conclude that (5) does not have purely imaginary roots
dominant eigenvalue. For the undelayed case, Proposition 1 for τ < τ ∗ . Thus, by [27, Theorem 2.1], all characteristic
p
p
can be equivalently stated as follows.
roots of (5) have strictly negative real parts for τp < τp∗ .
Corollary 1: If (H) holds, then all eigenvalues of L + cD
Remark 1: Proposition 2 provides an estimate for the
have negative real parts.
largest admissible pinning delay for which system (4) is stable.
We also state an easy observation for later use:
This estimate needs only the knowledge of the set of pinned
Lemma 1: For any two column vectors u, v ∈ Rn , det(In + nodes and their weighted in-degrees.
uv ⊤ ) = 1 + v ⊤ u.
IV. P INNING
III. E STIMATION OF

THE LARGEST ADMISSIBLE PINNING
DELAY

We first show that the system (4) is stable for all values of
the pinning delay τp smaller than a certain value τp∗ .
Proposition 2: Assume condition (H). Let
F (w, c, l, τ ) = c2 + ω 2 + 2c [l cos(ωτ ) − ω sin(ωτ )]

(6)

and define




τp∗ = sup τ : min min F (w, c, Lii , τ ) > 0 .
τ >0

ω∈R i∈D

(7)

If τp < τp∗ , then system (4) is stable for all τr ≥ 0.
Proof: First, we take τp = 0 and prove stability for
all τr ≥ 0. Assume for contradiction that there exists some

A SINGLE NODE

We now consider the possibility of controlling the network
using a single node, say, the qth one. Then D = uq u⊤
q , where
uq denotes the qth standard basis vector, whose qth component
is one and other components zero. If λIn + K − A exp(−λτr )
is nonsingular, the characteristic equation (5) becomes


χ(λ) = det λIn + K − A exp(−λτr ) + cuq u⊤
q exp(−λτp )

= det(λIn + K − A exp(−λτr ))


−1
det In + cuq u⊤
exp(−λτp )
q (λIn + K − A exp(−λτr ))
= det(λIn + K − A exp(−λτr ))

−1
(1 + cu⊤
uq exp(−λτp ))
q (λIn + K − A exp(−λτr ))
(10)

using Lemma 1. Then we have the following result.

3

Proposition 3: Assume (H). If all solutions λ of the equation
−1
1 + cu⊤
uq exp(−λτp ) = 0 (11)
q (λIn + K − A exp(−λτr ))

satisfy Re(λ) < 0, then system (4) is stable.
Proof: As in the first part of the proof of Proposition 2,
the equation det[λIn +K −A exp(−λτr )] = 0 has no solutions
with Re(λ) ≥ 0. Hence, if all solutions λ of (11) have negative
real parts, then all roots of (5) have negative real parts.
We consider two specific cases to obtain more information
about the solutions of (11). First, we consider the absence
of transmission delays, i.e., τr = 0. Suppose for simplicity that L is diagonalizable and has only real eigenvalues:
L = Q−1 JQ for some nonsingular Q and a real diagonal
matrix J = diag{θ1 , . . . , θn } of eigenvalues of L. The column
vectors of Q−1 (resp, the row vectors of Q) are the right (resp.,
left) eigenvectors of L. Then, (11) can be written as
1 + cζ ⊤ (λIn + J)−1 ξ exp(−λτp ) = 0,
⊤

u⊤
q Q

where ζ =
is the qth left eigenvector and ξ = Q
is the qth right eigenvector of L. We expand (12) as
1+c

n
X
ξi ζi exp(−λτp )

λ + θi

i=1

=0

(12)
−1

uq

(13)

in terms of the components ξi , ζi of ξ and ζ, respectively.
Consider the smallest value of τp for which there exists
a purely imaginary solution, λ = jω. Then, the real and
imaginary parts of (13) give
(
1 + a(ω) cos(ωτp ) − b(ω) sin(ωτp ) = 0
b(ω) cos(ωτp ) + a(ω) sin(ωτp )
=0
where
a(ω) = c

n
X
ξi ζi θi
,
ω 2 + θi2
i=1

b(ω) = c

X ξi ζi ω
.
ω 2 + θi2
i

(14)

Rearranging gives cos(ωτp ) = −a(ω)/(a2 (ω) + b2 (ω)) and
sin(ωτp ) = b(ω)/(a2 (ω) + b2 (ω)). This implies a(ω)2 +
b2 (ω) = 1 and
cos(ωτp ) = −a(ω),

sin(ωτp ) = b(ω).

(15)

We then have the following result.
Proposition 4: Suppose τr = 0, L is diagonalizable, irreducible, and all its eigenvalues are real. Let the eigenvalues
{θi } of L be sorted so that θq = 0, and let ζ = [ζ1 , . . . , ζn ],
P
n
k=1 ζk = 1, be the left eigenvector of L corresponding to
the zero eigenvalue. Let Z denote the set of positive solutions
of the equation
a2 (ω) + b2 (ω) = 1
(16)
2

with respect to the variable ω , where a(ω) and b(ω) are given
by (14). Define
√
arccos(−a( max Z))
M
√
.
(17)
τp =
max Z
Then system (4) is stable for τp < τpM .
Proof: Eq. (10) implies that any purely imaginary solution
jω of (5) should also be a solution of (13). Then ω must be

a real solution of (16). By the definition
√ of Z, the solution
set of (16) with respect to ω is {± z : z ∈ Z}. By
the assumption of irreducibility,
θi > 0 for all i 6= q and
√
ζi , ξi > 0 ∀i. If ω = z, then the smallest√positive
√ solution
of (15) with respect to τ√
p is arccos(−a( z))/ z. If, on
the other hand, ω = − z, noting that a(ω) > 0 and
b(ω) ≤ 0, √the smallest
positive solution of (15) is again
√
arccos(−a( z))/ z. Therefore, given ω 2 ∈ Z, the smallest
nonnegative solution√of (15)
√ with respect to τp should be in
z))/
z : z ∈ Z}. Since the mapping
the set {arccos(−a(
√
√
z 7→ arccos(−a( z))/ z is a decreasing function of z > 0,
the quantity τpM defined in (17) is the smallest nonnegative
solution of (15) with respect to τp , given ω 2 ∈ Z. Hence, for
τp < τpM (13) does not have any purely imaginary solutions.
Since for τp = 0 all characteristic roots of (5) have negative
real parts, we conclude that all roots have negative real parts
for τp < τpM .
Remark 2: By derivation, Eq. (13) is independent of the
ordering of the eigenvalues or the eigenvectors in J. Therefore,
the bound τpM for allowable pinning delays given in Proposition 4 does not depend on the choice of the pinned node.
Proposition 4 suggests an algorithm to calculate τpM :
1) Find the largest positive solution ω 2 of the equation
n
X
X ξi ξj ζi ζj (θi θj + ω 2 )
(ξk ζk )2
1
= 2 . (18)
+2
2
2
2
2
2
2
ω + θk
(ω + θi )(ω + θj )
c
i>j
k=1

2) Calculate (17).
We illustrate this approach in an Erdős-Renyi (E-R) random
network of n = 100 nodes with linking probability 0.03,
where the first node is pinned. The left and right eigenvectors of L associated
√ with the zero eigenvalue are given by
ζ = [1, . . . , 1]/ n. Figure 1 shows the parameter region
{(c, τp ) : τp < τpM }, illustrating the inverse dependence
of τpM on c. Note that τp > τpM does not necessarily
imply instability, since Proposition 4 gives only a sufficient
condition. Nevertheless, the curve shown in Fig. 1(a) turns
out to be a good approximation of the boundary of the exact
stability region. To illustrate, we take two parameter points
very close (±10% of the τpM ) to the curve but on different
sides of it, as indicated by blue and red stars in Fig. 1(a). We
simulate (3) at the corresponding parameter values, with the
same Laplacian as above and τr = 0. As seen in Fig. 1(b)–(c),
the two points indeed yield different stability properties.
The other situation we consider is the homogeneous case
when L is diagonalisable and normalised, i.e., Lii = l ∀i for
some l > 0, and τr = τp . Then (11) becomes
−1
1+cu⊤
uq exp(−λτp ) = 0. (19)
q ((λ+l)In −A exp(−λτr ))

Let L = QJQ−1 ; thus A = Q(lIn − J)Q−1 . Then, by the
same algebra as above, (19) becomes
1+c

n
X

k=1

ζk ξk exp(−λτp )
= 0.
(λ + l) + (θk − l) exp(−λτp )

(20)

We have the following result.
Proposition 5: Suppose that τr = τp , L is diagonalizable,
irreducible, normalised (Lii = l ∀i), and all its eigenvalues
{θi } are real. Denote θq = 0 and let ζ = [ζ1 , . . . , ζn ] be the

4

Let λi (c) denote the characteristic roots of (5) and φ̃i (c)
and ψ̃ i (c) denote the right and left eigenvectors of Σ(c, λi (c)),
regarded as functions of c, with λi (0) = σi , φ̃i (0) = φi and
ψ̃ i (0) = ψ i . Using a perturbation expansion [29], [30],

1.2

1

τp

0.8

0.6

0.4

Stable

0.2
4

5

6

7
c

8

9

λi (c)
ψ̃ i (c)

10

(a) The stability region {(c, τp ) : τp < τpM }.

where o(c) denotes terms that satisfy limc→0 |o(c)|/c = 0.
Thus,

3
2

yi, i=1,...,100

= σi + λ1i c + o(c), φ̃i (c) = φi + φi,1 c + o(c)
= ψ i + ψ i,1 c + o(c)

1

[−K + A exp(−λi (c)τr ) − cD exp(−λi (c)τp )]φ̃i (c)

0

= λi (c)φ̃i (c).

−1
−2
−3
0

10

20

30

40

50
Time

60

70

80

90

When c is sufficiently small, the dominant eigenvalue is λ1 (c),
since σ1 = 0 is the dominant eigenvalue when c = 0. Hence,
we consider i = 1. Then exp(−λ1 (c)τ ) = 1 − cλ11 τ + o(c).
Comparing the first-order terms in c on both sides, (−Aλ11 τr −
⊤
D)φi − Lφi,1 = λ11 φi . Multiplying both sides with ψ 1 and
⊤
noting that ψ 1 φ1 = 1,

100

(b) c = 4.48 and τp = 0.7724
3
2

yi,i=1,...,100

1
0

⊤

−1

λ11 = −

−2
−3
0

10

20

30

40

50
Time

60

70

80

90

100

(c) c = 4.48 and τp = 0.9441
Fig. 1. (a) The stability region {(c, τp ) : τp < τpM } in the parameter
plane (c, τp ), where the dashed line depicts τpM as a function of c. Direct
simulation verifies that the system is indeed stable for the parameter values
c = 4.48 and τp = 0.7724 (b), and unstable for the slightly different values
c = 4.48 and τp = 0.9441 (c), corresponding to the blue and red stars,
respectively, in subfigure (a).

left eigenvector of L corresponding to the eigenvalue 0, with
P
i ζi = 1. Let S denote the set of all the branches of the
solutions of the equation
1+c

n
X

k=1

ζk ξk
=0
exp(−lτp )s/τp + (θk − l)

(21)

with respect to the variable s. Then system (4) is stable
whenever the real parts of the numbers { Wτ(s)
− l : s ∈ S}
p
are all negative, where W is the Lambert W function [28].
Proposition 5 can be proved by transforming (20) into (21)
with s = τp (λ + l) exp(τp (λ + l)) and using Proposition 3.
V. S MALL

AND LARGE PINNING STRENGTHS

In this section, we consider the extreme situations when
the pinning strength c is very small or very large. We will
employ the perturbation approach in [29], [30] to approximate
the eigenvalues and eigenvectors in terms of c.
The characteristic roots λ of (5) are eigenvalues of the
matrix Σ(c, λ) = −K+A exp(−λτr )−cD exp(−λτp ). Hence,
when c = 0, the characteristic roots of (5) equal to the
eigenvalues {σi } of Σ(0, λ). Under the condition (H), there
is a single eigenvalue σ1 = 0. We denote the right and left
⊤
eigenvectors of Σ(0, σi ) by φi and ψ i respectively, with
⊤
ψ i φi = 1. It can be seen that ψ 1 and φ1 (associated with
σ1 = 0) are, respectively, the right and left eigenvectors of L
associated with the zero Laplacian eigenvalue.

ψ 1 Dφ1
1 + τr (ψ 1 ⊤ Aφ1 )

.

(22)

Hence, we have the following result.
Proposition 6: Suppose that the underlying graph is
strongly connected and at least one node is pinned. Then, for
sufficiently small c, all characteristic roots of (5) have negative
real parts and the dominant root is given by
⊤

λ1 (c) = −

ψ 1 Dφ1
1 + τr (ψ 1 ⊤ Kφ1 )

c + o(c).

(23)

Proof: Since the graph is strongly connected, L has a
simple zero eigenvalue. When c = 0, the dominant root of (5)
is σ1 = λ1 (0). Since the roots of (5) depend analytically on c,
they are given by λ1 (c) for all sufficiently small c. Substituting
⊤
(22) into λ1 (c) and noting that ψ 1 (−K + A)φ1 = 0
completes the proof.
In order to understand the meaning of (23), consider the
special case of an undirected graph with binary adjacency matrix A. Then, with φ1 = [1, . . . , 1]⊤ and ψ 1 = [1, . . . , 1]⊤ /n,
Pn
⊤
we have ψ 1 Kφ1 = i=1 Lii /n, which equals the average
Pn
⊤
degree of the graph. In addition, ψ 1 Dφ1 = i=1 δD (i)/n,
which is the fraction of pinned agents. Then, (23) yields the
approximation
λ1 (c) ≈ −

Pinning Fraction
c
1 + τr × Mean Degree

(24)

for small c, which uses only the pinning fraction and the
mean degree of the graph. Since the real part of the dominant
characteristic value measures the exponential convergence of
the system, Proposition 6 implies that, for sufficiently small
c, the convergence rate is improved if the number of pinned
nodes is increased, the transmission delay is reduced, or the
mean degree is decreased. If the graph is directed, a similar
1
statement can be obtained
Pn by1 taking the components of ψ as
1
1
weights: ψ Dφ = j=1 ψj δD (j).
To illustrate this result, we employ a numerical method to
calculate the real part of λ1 (c), namely, by simulating the

5

system (4) and expressing its exponential convergence rate
in terms of its largest Lyapunov exponent. In detail, letting
τm = max{τr , τp }, we partition time into disjoint intervals
of length τm , tk = kτm , and define ηk (θ) = y(tk + θ)
for θ ∈ [0, τm ]. Then, the largest Lyapunov exponent, which
equals to the largest real part of solutions of (5), is numerically
calculated via [31]

is bounded as c grows large. Thus, we assume that τpc := τp c
remains bounded as c → ∞.
When ǫ = 0, (27) becomes approximately ẋ = −Dx(t −
τp∞ ), where τp∞ can be any value between limc→∞ τpc and
limc→∞ τpc . In terms of components, ẋi = −xi (t − τp∞ ) if
i ∈ D, and 0 otherwise. The characteristic equation (27) with
ǫ = 0 can be written as

N
(µ + exp(−µτp∞ ))m µn−m = 0
(28)
1
1 X
kηk k
Re(λ1,sim ) = lim
log
log kηN k = lim
,
N →∞ N τm
N →∞ N τm
kηk−1where
k
m = |D|. It is known that Re(µ) < 0 for all roots of
k=1
(25) the function µ 7→ µ + exp(−µτp∞ ) if and only if τp∞ < π .
2
where k · k stands for the function norm. The latter is Therefore, we impose the condition: τp c < π .
2
numerically calculated by approximating ηk (·) with a finiteThus, the largest real part of the solutions of (28) is zero,
dimensional vector ϕk obtained by evaluating ηk at a finite and is obtained for the solution µ = 0. The corresponding
number of equally spaced points and using the vector norm eigenspace has dimension n − m and has the form
kϕk k. The estimate (25) can then be compared with the
ES = {u = [u1 , . . . , un ]⊤ ∈ Rn : ui = 0, ∀ i ∈ D}.
analytical estimate for Re(λ1 ) obtained from (23):
⊤

Re(λ1,est ) = −

ψ 1 Dφ1
1 + τr (ψ

1⊤

Kφ1 )

c.

(26)

For simulations, we generate an undirected E-R random
graph of n = 100 nodes with linking probability p = 0.03
and randomly select a given fraction f of them as the pinned
nodes. The pinning delay is taken as τp = 0.1. Figure 2 shows
that the simulated value of Re(λ1 ) decreases almost linearly
with respect to c and f , and increases with respect to τr and the
mean degree. The simulation results are in a good agreement
with the theoretical results. The error between Re(λ1,est ) and
Re(λ1,sim ) depends on the values of λ11 and c. It can be seen
that the error will increase as c or λ11 (or equivalently, f )
increases, or else as the mean degree or τr decreases.

Re(λ1)

1

Re(λ )

−0.1

−0.15
−0.2

−0.04
−0.06
−0.08
−0.1

−0.25 −3
10

−2

−1

10

0

10

c

−0.12
0

10

0.2

−0.05

−0.005

−0.1
−0.15
−0.2
4

6

8

10

12

14

Mean degress

(c)

0.4

0.6

Pinning Fraction

0.8

1

(b)
0

Re(λ1)

1

Re(λ )

(a)
0

16

18

µ1i ǫ + o(ǫ), ξ˜i (ǫ) = ξ i + ξ i,1 ǫ + o(ǫ)

ζ̃ i (ǫ) =

ζ i + ζ i,1 ǫ + o(ǫ)

(29)

as ǫ → 0. Thus, from (27),
h

i
τr 
−ǫK + ǫA exp −µi (ǫ)
− D exp(−µi (ǫ)τP c ) ξ˜i (ǫ)
ǫ
= µi (ǫ)ξ˜i (ǫ).

−0.015
0.2

0.4

[−K + A exp(−µ1i τr )]ξ i − Dξ i,1 = µ1i ξ i .
We write

K1
K=
0

−0.01

−0.02

µi (ǫ) =

Since exp(−µi (ǫ)τ ) = 1 − ǫµ1i τ + o(ǫ), by comparing the
coefficients of order 1, we have

−0.02

−0.05

Without loss of generality, we assume D = {1, . . . , m}. Thus,
we consider perturbation in terms of ǫ near zero eigenvalues µi
⊤
and its corresponding right and left vectors, ξ i , ζ i ∈ ES such
i ⊤ i
j ⊤ i
that (ζ ) ξ = 1 and (ζ ) ξ = 0 if i 6= j, i, j = m+1, . . . , n.
Let µi (ǫ) stand for the perturbed solution of (27), ξ̃ i (ǫ)
and ζ̃ i (ǫ) be the corresponding right and left eigenvectors,
respectively. By a perturbation expansion,

τ

0.6

0.8

1

r

(d)

Fig. 2. Variation of Re(λ1 ) with system parameters. The estimate (26) (plotted with +) shows good agreement with the values obtained via simulation
and (25) (plotted with ⋄). The parameters that are kept fixed are: (a) f = 0.3,
τr = 0.1, mean degree = 3.4; (b) c = 0.1, τr = 0.1, mean degree = 3.4;
(c) c = 0.1, f = 0.3, τr = 0.1; (d) c = 0.1, f = 0.3, mean degree = 3.4.

Next, we consider the case of large c. Letting ǫ = 1/c and
µ = λ/c, (5) is rewritten as
det [µIn + ǫK − ǫA exp(−µτr /ǫ) + D exp(−µτp /ǫ)] = 0.
(27)
By the foregoing results, one can see that when ǫ is sufficiently
small, equivalently, c is sufficiently large, the largest admissible pinning delay for (4) approaches zero. It is therefore
natural to assume that τp depends on c in such a way that τp c

0
K2
⊤



, A=

⊤



A11
A21

A12
A22
⊤



, D=



(30)

Im
0

0
0



,

⊤

and ξ i = [ξ1i , ξ2i ]⊤ , ξ i,1 = [ξ1i,1 , ξ2i,1 ]⊤ , with K1 , A11 ,
ξ1i = 0 and ξ1i,1 corresponding to the pinned subset D of
dimension m. Then (30) becomes
(
[−K2 + A22 exp(−µ1i τr )]ξ2i = µ1i ξ2i
(31)
= 0.
exp(−µ1i τr )A12 ξ2i − ξ1i,1
We have the following result.
Proposition 7: Suppose that the underlying graph is
strongly connected and at least one node is pinned. Fix τr ≥ 0,
and suppose τp c < π2 as c → ∞. Then the dominant root of
(27) has the form
λ(c) = µ1∗ + o(1) as c → ∞,

(32)

where µ1∗ is the dominant eigenvalue of the delay-differential
equation
ẏ = −K2 y(t) + A22 y(t − τr ).
(33)

6

−0.6

1

Re(λ )

−0.65
−0.7
−0.75
−0.8
−0.85
−0.9 1
10

2

c

10

Fig. 3. Variation of Re(λ1 ) with large values of c, calculated for f = 0.3,
τr = τp = 0.1, and mean degree = 3.4. The estimation Re(λ1,est ) is plotted
by the blue solid line and the real values by the dash line with red ⋄.

Furthermore, Re(λ(c)) < 0 for all sufficiently large c.
Proof: The condition τpc < π/2 implies that, when
ǫ = 0, the dominant root of the characteristic equation
(27) is zero and corresponds to the eigenspace ES. So, for
sufficiently small ǫ, the dominant root of equation (27) and
the corresponding eigenvector have the form (29), where µ1i
satisfies the first equation in (31), i.e., is an eigenvalue of (33).
Since λ(ǫ) = µ/ǫ, (32) follows. Moreover, since −K2 +A22 is
diagonally dominant, one can see that Re(µ1i ) < 0 under condition (H). Therefore, for sufficiently large c, all characteristic
values of system (3) have negative real parts.
We note that µ1∗ depends only on the coupling structure
of the uncoupled nodes. To illustrate this result, we
consider examples with a similar setup as in Sec. V.
We take an E-R graph with n = 100 nodes and linking
probability p = 0.03, and pin m = 30 nodes. We set
τr = 0.1 and τp = 1c . The real part of the dominant
characteristic root of (5) is numerically calculated via
the largest Lyapunov exponent, using formula (25). Its
theoretical
estimation comes from Theorem 7: Re(λ1,est ) =

max Re(µ1 ) : det µ1 Im + K2 − A22 exp(−µ1 τr ) = 0 ,
where the largest real part of µ1 is similarly calculated from
the largest Lyapunov exponent of (33). Fig. 3 shows that as c
grows large, the real part of the dominant root of (5) obtained
from simulations approach the theoretical result Re(λ1,est ),
thus verifying Proposition 7.
We have shown in this paper that the stability of the multiagent systems with a local pinning strategy and transmission
delay may be destroyed by sufficiently large pinning delays.
Using theoretical and numerical methods, we have obtained
an upper-bound for the delay value such that the system is
stable for any pinning delay less than this bound. In this case,
the exponential convergence rate of the multi-agent, which
equals the smallest nonzero real part of the eigenvalues of the
characteristic equation, measures the control performance.
R EFERENCES
[1] M. H. DeGroot, Reaching a consensus, J. Amer. Statist. Assoc., 69
(1974), 118–121.
[2] C. Reynolds, Flocks, herds, and schools: A distributed behavioral model,
Comput. Graph., 21:4 (1987), 25–34.
[3] T. Vicsek, A. Czirök, E. Ben-Jacob, I. Cohen and O. Shochet, Novel
type of phase transition in a system of self-driven particles, Phys. Rev.
Lett., 75 (1995), 1226–1229.
[4] A. Jadbabaie, J. Lin, and A. S. Morse, Coordination of groups of mobile
autonomous agents using nearest neighbor rules, IEEE Trans Automat.
Control, 48 (2003), 988–1001.
[5] J. A. Fax and R. M. Murray, Information flow and cooperative control of
vehicle formations, IEEE Trans. Autom. Control, 49 (2004), 1465–1476.

[6] X. F. Wang, G. Chen, Pinning control of scale-free dynamical network,
Physica A, 310(2002), 521–531.
[7] X. Li, X. F. Wang, G. Chen, Pinning a complex dynamical network to
its equilibrium, IEEE Trans. Circuits Syst. I, 51:10 (2004), 2074–2087.
[8] T. P. Chen, X. W. Liu, W. L. Lu. Pinning complex networks by a single
controller. IEEE Trans. Circuits Syst. I, 54:6 (2007), 1317–1326.
[9] W. L. Lu, X. Li, Z. H. Rong. Global stabilization of complex networks
with digraph topologies via a local pinning algorithm. Automatica, 46
(2010), 116–121.
[10] Q. Song, F. Liu, J. Cao, W. Yu, M-Matrix Strategies for PinningControlled Leader-Following Consensus in Multiagent Systems With
Nonlinear Dynamics, IEEE Trans. Cybern., 43:6 (2013), 1688–1697.
[11] Q. Song, F. Liu, J. Cao, W. Yu, Pinning-Controllability Analysis of
Complex Networks: An M-Matrix Approach, IEEE Trans. Circuits Syst.I, 59:11(2012) 2692–270.
[12] Q. Song, J. Cao, On Pinning Synchronization Of Directed And Undirected Complex Dynamical Networks, IEEE Trans. Circuits Syst.-I,
57:3(2010), 672–680.
[13] R. Olfati-Saber and R. M. Murray, Consensus problems in networks of
agents with switching topology and time-delays, IEEE Trans. Autom.
Control, 49 (2004), 1520–1533.
[14] F. M. Atay, Consensus in networks under transmission delays and the
normalized Laplacian. Phil. Trans. Roy. Soc. A, 371 (2013), 20120460.
[15] F. M. Atay, On the duality between consensus problems and Markov
processes, with application to delay systems. Markov Processes and
Related Fields (in press).
[16] P.-A. Bliman and G. Ferrari-Trecate, Average consensus problems in networks of agents with delayed communications, Automatica, 44 (2008),
1985–1995.
[17] L. Moreau, Stability of multiagent systems with time-dependent communication links, IEEE Trans. Autom. Control, 50 (2005), 169–182.
[18] F. Xiao and L. Wang, Asynchronous consensus in continuous-time multiagent systems with switching topology and time-varying delays, IEEE
Trans. Autom. Control, 53 (2008), 1804–1816.
[19] W. L. Lu, F. M. Atay, J. Jost, Consensus and synchronization in discretetime networks of multi-agents with stochastically switching topologies
and time delays, Networks and Heterogeneous Media, 6:2(2011), 329–
349.
[20] L. Xiang, Z. Chen, Z. Liu et al. Pinning control of complex dynamical
networks with heterogeneous delays, Computers & Mathematics with
Applications, 56:5 (2008), 1423–1433.
[21] M. Ulrich, A. Papachristodoulou, F. Allgöwer, Generalized Nyquist
consensus condition for linear multi-agent systems with heterogeneous
delays, Proc. IFAC Workshop on Estimation and Control of Networked
Systems, 24–29, 2009.
[22] P. Andrey, H. Werner, Robust stability of a multi-agent system under arbitrary and time-varying communication topologies and communication
delays, IEEE Trans. Autom. Control, 57:9(2012), 2343–2347.
[23] Z.X. Liu, Z.Q. Chen, Z.Z. Yuan, Pinning control of weighted general
complex dynamical networks with time delay, Physica A, 375:1 (2007),
345–354.
[24] J. Zhao, J.-A. Lu, and Q. Zhang, Pinning a complex delayed dynamical
network to a homogeneous trajectory, IEEE Trans. Circuits Syst.-II, 56:6
(2009), 514–518.
[25] W. Gu, Lag synchronization of complex networks via pinning control,
Nonlinear Analysis: Real World Applications, 12:5 (2011), 2579–2585.
[26] Q. Song, J. Cao, F. Liu, W. Yu, Pinning-controlled synchronization of
hybrid-coupled complex dynamical networks with mixed time-delays,
Inter. J. Robust Nonlin. Control, 22(2012), 690–706
[27] S. Ruan, J. Wei. On the zeros of transcendental functions with applications to stability of delay differential equations with two delays. Dynam.
Contin. Discrete Impuls. Systems, Ser. A 10 (2003), 863–874.
[28] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, D. E. Knuth.
On the Lambert W function. Adv. Comput. Math. 5 (1996), 329–359.
[29] L. N. Trefethen. Numerical Linear Algebra. SIAM (Philadelphia, PA),
p. 287 1997.
[30] R.-C. Li. Matrix Perturbation Theory. In Hogben, Leslie. Handbook of
Linear Algebra (2nd ed.), 2014.
[31] J. D. Farmer, Chaotic attractors of an infinite-dimensional dynamical
system, Physica D, 4:3 (1982), 366–393.

