Didier Dubois
email: dubois@irit.fr

Henri Prade
email: prade@irit.fr

Formal Concept Analysis from the Standpoint of Possibility Theory

Formal concept analysis (FCA) and possibility theory (PoTh) have been developed independently. They address different concerns in information processing: while FCA exploits relations linking objects and properties, and has applications in data mining and clustering, PoTh deals with the modeling of (graded) epistemic uncertainty. However, making a formal parallel between FCA and PoTh is fruitful. The four set-functions at work in PoTh have meaningful counterparts in FCA; this leads to consider operators neglected in FCA, and thus new fixed point equations. One of these pairs of equations, paralleling the one defining formal concepts in FCA, defines independent sub-contexts of objects and properties that have nothing in common. The similarity of the structures underlying FCA and PoTh is still more striking, using a cube of opposition (a device extending the traditional square of opposition in logic). Beyond the parallel between FCA and PoTh, this invited contribution, which largely relies on several past publications by the authors, also addresses issues pertaining to the possible meanings, degree of satisfaction vs. degree of certainty, of graded object-property links, which calls for distinct manners of handling the degrees. Other lines of interest for further research are briefly mentioned.

Introduction

Formal concept analysis (FCA) and possibility theory (PoTh) are two theoretical frameworks that are addressing different concerns in the processing of information. Namely FCA builds concepts from a relation linking objects to the properties they satisfy, which has applications in data mining, clustering and related fields, while PoTh deals with the modeling of (graded) epistemic uncertainty. This difference of focus explains why the two settings have been developed completely independently for a very long time. However, it is possible to build a formal analogy between FCA and PoTh. Both theories heavily rely on the comparison of sets, in terms of containment or overlap. The four set-functions at work in PoTh actually determine all possible relative positions of two sets. Then the FCA operator defining the set of objects sharing a set of properties, which is at the basis of the definition of formal concepts, appears to be the counterpart of the set function expressing strong (or guaranteed) possibility in PoTh. Then, it suggests that the three other set functions existing in PoTh should also make sense in FCA, which leads to consider their FCA counterparts and new fixed point equations in terms of the new operators. One of these pairs of equations, paralleling the one defining formal concepts, define independent sub-contexts of objects and properties that have nothing in common.

The parallel of FCA with PoTh can still be made more striking using a cube of opposition (a device extending the traditional square of opposition existing in logic, and exhibiting a structure at work in many theories aiming at representing some aspects of the handling of information).

In this survey we shall indicate various issues pertaining to FCA that could be worth studying in the future. For instance, the object-property links in formal contexts of FCA may be a matter of degree. These degrees may refer to very different notions, such as the degree of satisfaction of a gradual property, the degree of certainty that an object has, or not, a property, or still the typicality of an object with respect to a set of properties. These different intended semantics call for distinct manners of handling the degrees, as advocated in the presentation.

Lastly, other examples of lines of interest for further research, such as the extension of the parallel of FCA with PoTh to conceptual pattern structures, or the applications to the fusion of conflicting pieces of information, to the clustering of sets of objects on the basis of approximate concepts, or to the building of conceptual analogical proportions, are briefly mentioned.

Possibility Theory and Formal Concept Analysis -A Parallel

Formal concept analysis [START_REF] Barbut | Ordre et Classification. Algèbre et Combinatoire[END_REF][START_REF] Ganter | Formal Concept Analysis[END_REF][START_REF] Wille | Restructuring lattice theory: an approach based on hierarchies of concepts[END_REF] associates objects with the set of their properties, through a formal context which is a binary relation R on the Cartesian product of the set of objects O and the set of properties P. Thus, knowing only that an object x has some property y, the set R t (y) = {x ∈ O|(x, y) ∈ R} is the set of the possible objects corresponding to the elementary piece of knowledge "the object has property y"(in the context R). This suggests a possibilistic reading of formal concept analysis and leads to considering the formal counterpart to possibility theory set-functions in this framework. After introducing some notations, we first provide a short refresher on possibility theory [START_REF] Dubois | Possibility Theory: An Approach to Computerized Processing of Uncertainty[END_REF][START_REF] Dubois | Possibility theory[END_REF][START_REF] Zadeh | Fuzzy sets as a basis for a theory of possibility[END_REF].

Describing Objects

An object, or item, is denoted by x, or x i in case we consider several ones at the same time. A subset of objects is denoted by a capital letter X, and we write X = {x 1 , . . . , x i , . . . , x m }. A set of objects associated with their respective sets of properties defines a formal context R ⊆ O × P [START_REF] Ganter | Formal Concept Analysis[END_REF]. An object x is associated with its description, denoted ∂(x). In the following, we only consider simple descriptions, expressible in terms of a subset Y of properties y j , namely, Y = {y 1 , . . . , y j , . . . , y n }. In such a case, we write

∂(x) = Y .
Besides, a useful kind of structured description is in terms of attributes. Let a, and A = {a 1 , . . . , a k , . . . , a r }, respectively denote an attribute, and a set of attributes. The value of attribute a for x is denoted a(x) = u, where u belongs to the attribute domain U a . In this case, we shall write ∂(x) = (a 1 (x), . . . , a k (x), . . . , a r (x)) = (u 1 , . . . , u k , . . . , u r ). This corresponds to a completely informed situation where all the considered attribute values are known for x. When it is not the case, the precise value a k (x) will be replaced by the possibility distribution π a k (x) . Such a possibility distribution [START_REF] Zadeh | Fuzzy sets as a basis for a theory of possibility[END_REF] is a mapping from U a k to [0, 1], or more generally any linearly ordered scale. Then π a k (x) (u) ∈ [0, [START_REF] Amgoud | A formal concept view of abstract argumentation[END_REF] estimates to what extent it is possible that the value of a k for x is u. 0 means impossibility; several distinct values may be fully possible (i.e. at degree 1). The characteristic function of an ordinary subset is a particular case of a possibility distribution. Precise information corresponds to the characteristic function of singletons. An elementary property y can be viewed as a subset of a single attribute domain, i.e. y ⊆ U . Note that while a set of properties Y is conjunctive (in the sense that an object possesses all properties in Y ), each property y corresponds to a subset of some attribute domain U that is disjunctive [START_REF] Dubois | Gradualness, uncertainty and bipolarity: making sense of fuzzy sets[END_REF]: it is a set of mutually exclusive values, since object x having property y possesses a single attribute value a(x) = u in U .

Taking inspiration from the existence of four set functions in possibility theory [START_REF] Dubois | Possibility theory: qualitative and quantitative aspects[END_REF], new operators have been suggested in the setting of formal concept analysis [START_REF] Dubois | A possibility-theoretic view of formal concept analysis[END_REF]. These set functions are now recalled, emphasizing the symmetrical roles played by the object x and the attribute value u, a point of view unusual in possibility theory, but echoing the symmetrical role played by objects and properties in formal concept analysis. See [START_REF] Dubois | Possibility theory: qualitative and quantitative aspects[END_REF][START_REF] Dubois | Possibility theory[END_REF] for more complete introductions and surveys on possibility theory.

Possibility Theory

Let π a(x) (u) denote the possibility that object x has value u ∈ U (for attribute a). For simplicity, we only consider the single-attribute case here. We assume that π a is bi-normalized: ∀x ∃u π a(x) (u) = 1 and ∀u ∃x π a(x) (u) = 1. This means that for any object x, there is some fully possible value for attribute a, and that for any value u there is an object x that takes this value. Let X be a set of objects, and y ⊆ U be a property. Then, one can define (i) the possibility measures [START_REF] Zadeh | Fuzzy sets as a basis for a theory of possibility[END_REF], denoted by Π:

Π(X) = max x∈X π a(x) (u) and Π(y) = max u∈y π a(x) (u).

Π(X)

estimates to what extent it is possible that there is an object in X having value u, while Π(y) is the possibility that object x has property y. Π is an indicator of non-empty intersection of the fuzzy set induced by the possibility distribution with an ordinary subset. They are measures of "weak, or potential possibility". Clearly, Π is max-decomposable with respect to set union.

(ii) the dual measures of necessity N (or "strong or actual necessity") [START_REF] Dubois | Fuzzy Sets and Systems: Theory and Applications[END_REF]:

N (X) = min x ∈X 1 -π a(x) (u) and N (y) = min u ∈y 1 -π a(x) (u) N (X)
estimates to what extent it is certain (necessarily true) that an object has value u is in X, while N (y) is the certainty that object x has property y. Note that N (y) = 1 -Π(y) where y = U \ y. N may be viewed as a degree of inclusion of the fuzzy set induced by the possibility distribution into an ordinary subset. N is min-decomposable with respect to set intersection.

(iii) the measures of "strong (or actual, or guaranteed) possibility" [START_REF] Dubois | Possibility theory as a basis for preference propagation in automated reasoning[END_REF] ∆(X) = min x∈X π a(x) (u) and ∆(y) = min u∈y π a(x) (u) ∆(X) estimates to what extent it is possible that all objects in X have value u, while ∆(y) estimates the possibility that object x takes any value in y. ∆ may be viewed as a degree of inclusion of an ordinary subset into the fuzzy set induced by the possibility distribution. ∆ is min-decomposable with respect to set union.

(iv) the dual measures of "weak (or potential) necessity or certainty" [START_REF] Dubois | Possibility theory as a basis for preference propagation in automated reasoning[END_REF]

∇(X) = 1 -min x ∈X π a(x) (u) and ∇(y) = 1 -min u ∈y π a(x) (u)
∇(X) estimates to what extent there exists at least one object outside X that has a low degree of possibility of having value u, while ∇(y) measures to what extent x has a low possibility value outside y. Note that ∇(y) = 1 -∆(y). ∇ is an indicator of non-full coverage of the considered universe by the fuzzy set induced by the possibility distribution together with an ordinary subset. ∇ is max-decomposable with respect to set intersection.

Formal Context Setting

The classical setting of formal concept analysis defined from a formal context relies on a single operator that associates a subset of objects with the set of properties shared by them (and the dual operator). In [START_REF] Dubois | A possibility-theoretic view of formal concept analysis[END_REF], this framework has been enlarged with the introduction of three other operators. We now recall the four operators which are counterparts to the possibility theory set functions in the setting of a formal context.

Namely, let R be the formal context. Then R(x) = {y ∈ P|(x, y) ∈ R} is the set of properties of object x, and R t (y) = {x ∈ O|(x, y) ∈ R} is the set of objects having properties y. Then, four remarkable sets can be associated with a subset X of objects (the notations have been chosen here in order to emphasize the parallel with possibility theory):

-the set R Π (X) of properties that are possessed by at least one object in X:

R Π (X) = {y ∈ P|R t (y) ∩ X = ∅} = x∈X R(x). Clearly, we have R Π (X 1 ∪ X 2 ) = R Π (X 1 ) ∪ R Π (X 2 ).
the set R N (X) of properties s. t. any object that satisfies one of them is necessarily in X:

R N (X) = {y ∈ P|R t (y) ⊆ X} = x ∈X R(x).
In other words, having any property in R N (X) is a sufficient condition for belonging to X. Moreover, we have

R N (X) = R Π (X) = P \ R Π (X), and R N (X 1 ∩ X 2 ) = R N (X 1 ) ∩ R N (X 2 ). -the set R △ (X) of properties shared by all objects in X: R △ (X) = {y ∈ P|R t (y) ⊇ X} = x∈X R(x).
In other words, satisfying all properties in R △ (X) is a necessary condition for an object for belonging to X. R △ (X) is a partial conceptual characterization of objects in X: objects in X have all the properties of R △ (X) and may have some others (that are not shared by all objects in X). It is worth noticing that R Π (X) provides a negative conceptual characterization of objects in X since it gathers all the properties that are never satisfied by any object in X.

Moreover, we have R △ (X 1 ∪ X 2 ) = R △ (X 1 ) ∩ R △ (X 2 ). Besides, as can be seen, R N (X) ∩ R △ (X)
is the set of properties possessed by all objects in X and only by them.

-the set R ∇ (X) of properties that are not satisfied by at least one object in X.

R ▽ (X) = {y ∈ P|R t (y) ∪ X = O} = x ∈X R(x). Note that R ▽ (X) = R △ (X) = P \ R △ (X).
In other words, in context R, for any property in R ▽ (X), there exists at least one object outside X that misses it. Moreover, we have

R ▽ (X 1 ∩ X 2 ) = R ▽ (X 1 ) ∪ R ▽ (X 2 ).
Note that R Π (X) and R N (X) become larger when X increases, while R △ (X) and R ▽ (X) get smaller. The four subsets R Π (X), R N (X), R △ (X), and R ▽ (X) have been considered by different authors (with different notations) without any reference to possibility theory. Düntsch et al. [START_REF] Düntsch | Approximation operators in qualitative data analysis[END_REF][START_REF] Düntsch | Mixing modal and sufficiency operators[END_REF] calls R △ a sufficiency operator, and its representation capabilities are studied in the theory of Boolean algebras. Taking inspiration as the previous authors from rough sets [START_REF] Pawlak | Rough Sets. Theoretical Aspects of. Reasoning about Data[END_REF], Yao [START_REF] Yao | Rough set approximations in formal concept analysis[END_REF][START_REF] Yao | A comparative study of formal concept analysis and rough set theory in data analysis[END_REF] also considers these four subsets. In both cases, the four operators were introduced. See also [START_REF] Georgescu | Non-dual fuzzy connections[END_REF][START_REF] Popescu | A general approach to fuzzy concepts[END_REF].

The Cube of Opposition in FCA

Before being able to present the structures of opposition relating the four operators introduced in the previous section, we need to start with a refresher on the Aristotelian square of opposition [START_REF] Parsons | The traditional square of opposition[END_REF]. The traditional square involves four logically related statements exhibiting universal or existential quantifications: it has been noticed that a statement A of the form "every x is p"is negated by the statement O "some x is not p", while a statement like E "no x is p"is clearly in even stronger opposition to the first statement A. These three statements, together with the negation of the last one, namely I "some x is p", give birth to the Aristotelian square of opposition in terms of quantifiers A: ∀x p(x), E: ∀x ¬p(x), I: ∃x p(x), O: ∃x ¬p(x), pictured in Fig. 1. Such a square is usually denoted by the letters A, I (affirmative half) and E, O (negative half). The names of the vertices come from a traditional Latin reading: AffIrmo, nEgO). Recently, it has been noticed that such a square can be generated by a binary relation and a subset that can be composed together [START_REF] Ciucci | The structure of oppositions in rough set theory and formal concept analysis -toward a new bridge between the two settings[END_REF]. Indeed, let R be a binary relation on a Cartesian product X × Y (nothing forbids Y = X in the construction we are going to describe). We assume R = ∅. Let R t denote the transposed relation ((y, x) ∈ R t iff (x, y) ∈ R). Moreover, we assume that ∀x, R(x) = ∅, which means that the relation R is serial, namely ∀x, ∃y such that (x, y) ∈ R; this is also referred to in the following as the X -normalization condition. In the same way R t is also supposed to be serial, i.e., ∀y, R t (y) = ∅ (Ynormalization). We further assume that the complementary relation R ((x, y) ∈ R iff (x, y) ∈ R), and its transpose are also serial, i.e. ∀x, R(x) = Y and ∀y, R t (y) = X . These conditions enforce a non trivial relation between X and Y. In the following, set complementations will be denoted by means of overbars.

Let S be a subset of Y. We assume S = ∅ and S = Y. The relation R and the subset S give birth to the following subset of X, namely the (left) image of

S by R R(S) = {x ∈ X | ∃s ∈ S, (x, s) ∈ R} = {x ∈ X | S ∩ R(x) = ∅}.
Similarly, we consider R(S), R(S), and

R(S) = {x ∈ X | ∀s ∈ S, (x, s) ∈ R} = {x ∈ X | R(x) ⊆ S}.
The four subsets thus defined can be nicely organized into a square of opposition. See Fig. 2. Indeed, it can be checked that the set counterparts of the relations existing between the logical statements of the traditional square of oppositions still hold here. Namely, R(S) and R(S) are complements of each other, as R(S) and R(S); we have R(S) ⊆ R(S) and R(S) ⊆ R(S), thanks to Let us now consider the complementary relation R. We further assume that R = ∅ (i.e., R = X × Y). Moreover we have also assumed the X -normalization of R, i.e. ∀x, ∃y (x, y) ∈ R. In the same way as previously, we get four other subsets of

X -normalization condition; R(S) ∩ R(S) = ∅ ; R(S) ∪ R(S) = X .
X from R. Namely, R(S) = {x ∈ X | ∃s ∈ S, (x, s) ∈ R} = {x ∈ X | S∪R(x) = X }; R(S); R(S); R(S) = {x ∈ X | ∀s ∈ S, (x, s) ∈ R} = {x ∈ X | S ⊆ R(x)}.
This generates a second square of opposition denoted by aeoi.

As can be seen, when R is a formal context (i.e., X = O, Y = P), we have R Π (S) = R(S), R N (S) = R(S), R ∆ (S) = R(S), R ∇ (S) = R(S). The eight subsets involving R and its complement can be organized into a cube of opposition as in Fig. 3. The four formal concept analysis operators correspond to the left side facet of the cube of oppositions. The full cube is then obtained by introducing their complements, giving birth to the right side facet. Since R Π (S) = R N (S), and R ∆ (S) = R ∇ (S), the classical square of oppositions AEOI is given by the four corners R N (S), R N (S), R Π (S), and R Π (S), and the second square aeoi on the back of the cube is given by R ∆ (S), R ∆ (S), R ∇ (S), and R ∇ (S).

Moreover, in the side facets, all edges are uni-directed, including the diagonal ones, and express inclusions. Indeed, as already established in [START_REF] Dubois | A possibility-theoretic view of formal concept analysis[END_REF], under the Xand Y-normalization hypotheses, the following inclusion relation holds:

R N (S) ∪ R ∆ (S) ⊆ R Π (S) ∩ R ∇ (S).
R N (S), R ∆ (S), R Π (S), and R ∇ (S) constitute, four distinct pieces of information [START_REF] Dubois | A possibility-theoretic view of formal concept analysis[END_REF], which are only (weakly) related by the above relation.

Lastly, it can be checked that we also have R ∆ (S) ∩ R N (S) = ∅ and R ∆ (S) ∩ R N (S) = ∅ on the one hand, and R ∇ (S) ∪ R Π (S) = X and R ∇ (S) ∪ R Π (S) = X on the other hand. These are the relations that holds on the top and on the bottom facets of the cube respectively. The cube of oppositions not only underlie FCA (and PoTh) [START_REF] Dubois | From Blanché's hexagonal organization of concepts to formal concept analysis and possibility theory[END_REF], but also is a setting of interest for building bridges with rough set theory [START_REF] Pawlak | Rough Sets. Theoretical Aspects of. Reasoning about Data[END_REF] (see [START_REF] Ciucci | The structure of oppositions in rough set theory and formal concept analysis -toward a new bridge between the two settings[END_REF]), or even formal argumentation [1]!

Formal Context Decomposition

In FCA, a formal concept [START_REF] Ganter | Formal Concept Analysis[END_REF] is defined as a pair (X, Y ) ∈ O × P such that

R △ (X) = Y and R t△ (Y ) = X, where R t△ (Y ) = {x ∈ O|R(x) ⊇ Y } = y∈Y R t (y)
is the set X of objects having all properties in Y , and in this case Y is also the maximal set of properties shared by all objects in X. A formal concept (X, Y ) is a maximal sub-rectangle in the formal context, i.e. is such that X × Y ⊆ R. It can be checked that R ▽ gives back the same Galois connection as the one defined from R △ , while R N (or R Π ) induces another connection, which is now described.

Consider the connection defined from R N in a similar formal way as when defining formal concepts. It was proposed by Popescu [START_REF] Popescu | A general approach to fuzzy concepts[END_REF] and studied in a general setting of residuated algebras, but not in the usual Boolean setting. Namely, let us consider pairs (X, Y ) s.t. R N (X) = Y and R tN (Y ) = X. As suggested in [START_REF] Dubois | Possibility theory and formal concept analysis in information systems[END_REF], the pairs (X, Y ) s.t. R N (X) = Y and R tN (Y ) = X allow us to characterize independent sub-contexts (i.e. that have no common objects and no common properties). They are thus of interest for the decomposition of a formal context into smaller independent ones. This is expressed through the following property, proved in [START_REF] Djouadi | Possibility theory and formal concept analysis: context decomposition and uncertainty handling[END_REF][START_REF] Dubois | Possibility theory and formal concept analysis: characterizing independent sub-contexts[END_REF]:

1 Proposition 1. The following properties of pairs (X, Y ) are equivalent 1. R N (X) = Y and R tN (Y ) = X 2. R N (X) = Y and R tN (Y ) = X 3. R Π (X) = Y and R tΠ (Y ) = X 4. R ⊆ (X × Y ) ∪ (X × Y )
Proof. Let us first show that Property 1 implies Property 4. First it is clear that:

R N (X) = Y ⇔ x∈X R(x) = Y ⇔ x∈X R(x) = Y . Denoting X + Y = X × Y , it implies R ⊆ X + Y . Likewise due to R tN (Y ) = X, R t ⊆ Y + X holds.
Finally: R ⊆ (X +Y )∩(Y +X), which equivalently writes: R ⊆ (X ×Y )∪(X ×Y ).

Conversely assume Property 4. Then it is clear that R N (X) ⊆ Y and R tN (Y ) ⊆ X hold since there is no property possessed by any object in X outside Y , and no object outside X that possesses a property outside Y . Suppose R N (X) ⊂ Y , i.e. ∃y * ∈ Y such that property y * is possessed by objects outside X. But then R(x, y * ) = 1 for some x ∈ X, y ∈ Y . So Property 4 does not hold. Contradiction.

The invariance of Property 4 with respect to complementation proves that the choice of (X, Y ) versus (X, Y ) in Property 1 is immaterial. Hence the equivalence with Property 2. For Property 3, note that R

N (X) = Y is equivalent to R Π (X) = Y . ⊓ ⊔
Thus, (X, Y ) and (X, Y ) are two independent sub-contexts in R, in the sense that there is no object / property pair (x, y) of the context R either in X × Y or in X × Y . The above proposition does not involve any minimality in the inclusion Property 4 of the above proposition. In particular, the pair (O, P) trivially satisfies it. However, this result leads to a decomposition of R into a disjoint union of minimal independent sub-contexts. Indeed, suppose two pairs (X 1 , Y 1 ), (X 2 , Y 2 ) satisfy Proposition 1. It implies that for instance, the pair

(X 1 ∩ X 2 , Y 1 ∩ Y 2 ) satisfies it (it can be checked that R N (X 1 ∩ X 2 ) = Y 1 ∩ Y 2 ),
and likewise with any element of the partition refining both partitions (X 1 , X 1 ) and (X 2 , X 2 ). Due to point 4 of Proposition 1, it yields

R ⊆ ((X 1 × Y 1 ) ∪ (X 1 × Y 1 )) ∩ ((X 2 × Y 2 ) ∪ (X 2 × Y 2 )),
where the intersection on the right-hand side comes down to the union of subcontexts (

X 1 ∩ X 2 ) × (Y 1 ∩ Y 2 ), (X 1 ∩ X 2 ) × (Y 1 ∩ Y 2 ), (X 1 ∩ X 2 ) × (Y 1 ∩ Y 2 ), (X 1 ∩ X 2 ) × (Y 1 ∩ Y 2 )
. The decomposition of R into minimal subcontexts is achieved by taking the following intersection

(X,Y ):R N (X)=Y,R tN (Y )=X (X × Y ) ∪ (X × Y ).
1 We again provide the proof for the sake of self-containedness. Thus, through the notions of formal sub-contexts and of formal concepts, one sees two key aspects of granulation at work. Namely, on the one hand independent sub-contexts are separated, while inside each sub-context, formal concepts (X, Y ) are identified where each object in X is associated with each property in Y . However, objects in the extension of a formal concept may not be fully similar since they may also possess properties outside the intension of the concept. They are only similar with respect to the properties associated to the formal concept. Thus, the classical Galois connection founding formal concept analysis (associated with the actual possibility operator), and the other connection induced by the actual necessity operator, respectively embed two basic ideas associated with the idea of a cluster (see, e.g., [START_REF] Höppner | Fuzzy Cluster Analysis: Methods for Classification. Data Analysis and Image Recognition[END_REF]), namely 1. any pair of elements in a cluster should be closely related in some sense, and 2. any element of a cluster should be sufficiently separated from any element outside it.

Moreover, formal concept analysis is also useful for conceptual clustering, where clusters should be associated with labels, obtained in this case as a conjunction of the properties shared by the objects in the cluster [START_REF] Carpineto | GALOIS: an order-theoretic approach to conceptual clustering[END_REF]. Such an idea can be also stated in terms of graph clustering, taking advantage of an exact parallel between formal concept analysis and bipartite graph analysis [START_REF] Gaume | A parallel between extended formal concept analysis and bipartite graphs analysis[END_REF], as viewing an (ideal) cluster as a group of vertices 1. either with no missing link inside the group, 2. or with no link with vertices outside the group.

These two complementary views are also clearly at the basis of cluster analysis for unipartite graphs [START_REF] Schaeffer | Graph clustering[END_REF].

In practice, it is important to introduce some tolerance in the evaluation of the similarity between the members of a cluster and in the separatedness of the clusters, leading to a more permissive and approximate view of granules or clusters; see, e.g., [START_REF] Gaume | Clustering bipartite graphs in terms of approximate formal concepts and sub-contexts[END_REF].

The other (mixed) connections R Θ (S) = T and R tΛ (T ) = S where Θ, Λ ∈ {Π, N, ∆, ∇} with Θ = Λ are also worth studying. They have still to be better understood and to be investigated systematically. See [START_REF] Djouadi | Possibility-theoretic extension of derivation operators in formal concept analysis over fuzzy lattices[END_REF] for a preliminary discussion, and [START_REF] Bělohlávek | Optimal triangular decompositions of matrices with entries from residuated lattices[END_REF][START_REF] Bělohlávek | Discovery of optimal factors in binary data via a novel method of matrix decomposition[END_REF] for results in the graded case.

Graded Links Between a Property and an Object

Fuzzy extensions of FCA where R is a fuzzy relation in L O×P with L often taken as the unit interval have been proposed early [START_REF] Bělohlávek | Fuzzy Galois connections[END_REF][START_REF] Burusco | The study of the L-fuzzy concept lattice[END_REF]. However, the development of a fuzzy formal concept analysis theory requires an appropriate algebra of fuzzy sets [START_REF] Bělohlávek | Fuzzy Galois connections[END_REF][START_REF] Bělohlávek | What is a fuzzy concept lattice[END_REF]. While many theoretical studies have been developed, the different gradual interpretations of a fuzzy formal context have not been much discussed. Following [START_REF] Djouadi | Graduality, uncertainty and typicality in formal concept analysis[END_REF], this section highlights some basic issues regarding the fact that a "fuzzy" or graded extension of binary formal contexts may convey different semantics: graded satisfaction of properties vs. uncertainty.

Gradual Properties: Unipolar Vs. Bipolar Scale Interpretation

In this first interpretation, the values in the table (which are scalars in L) may be understood as providing a refinement of the cross marks. Namely, they represent to what extent an object has a property, while in the classical model, this relationship was not a matter of degree. It is important to remark that in this view, we do not refine the absence of a property for an object (the blank mark is always replaced by the bottom element 0 of L). This view will be referred to as the positive unipolar interpretation. In this interpretation, R t (y) (resp. R(x)) is considered as the support of the fuzzy set of objects (resp. properties) satisfying the property y (resp. the object x). One could also consider the opposite convention namely the negative unipolar interpretation where degrees would represent to which extent an object does not have a property and equivalently provide a refinement of the blank marks.

The most commonly used interpretations, through existing FCA proposals, are implicitly based on the positive unipolar interpretation that allows to map a formal context with quantitative attributes into a fuzzy formal context. In this spirit, conceptual scale theory [START_REF] Wolff | Concepts in fuzzy scaling theory: order and granularity[END_REF] may be used to achieve a suitable (Boolean) representation by successive subsumptions.

Example 2. For instance, the formal context illustrated in Table 2 is obtained from Table 1 by a conceptual scaling of both many-valued attributes "Age"and "Salary". As can be seen, we have two sets of properties with obvious subsumption relations between them. Pairs ({P eter, Sophie, M ike, Joe}, {age ≥ 20, salary ≥ 1000}), ({Sophie, M ike}, {age ≥ 20, age ≥ 25, salary ≥ 1000, salary ≥ 1200}), or ({M ike}, {age ≥ 20, age ≥ 25, age ≥ 30, salary ≥ 1000, salary ≥ 1200, salary ≥ 1400}) are formal concepts. Knowing the ages and the salaries, the formal context R 2 can be re-encoded in a more compact way, using two fuzzy sets 'young' and 'small' with decreasing membership functions, as illustrated in Table 3.

Observe also that R 3 offers a more precise representation of initial data than Table 2. The context in Table 3, event though more compact than Table 2 highlights the fact that Mike, and to a lesser extent Sophie are not very young and have a salary that is not really low. It constitutes in some sense the negative of the picture shown on Table 1. Note that the type of representation on Table 3 can be obtained even without providing interpretable fuzzy sets and thus, by normalizing in L the domain of attribute values. This approach is used in [START_REF] Messai | Many-valued concept lattices for conceptual clustering and information retrieval[END_REF].

Another interpretation of the degrees, maybe more in the standard spirit of fuzzy logic would be to replace both the cross marks and the blank marks by values in the scale L (L = [0, 1]). Then L possesses a mid-point acting as a pivoting value between the situations where the object possesses the property to some extent and the converse situation where the object possesses the opposite property to some extent. Under this view, a fuzzy formal concept should be learnt together with its negation. This view corresponds to a bipolar scale interpretation.

Uncertainty

Neither the standard FCA approach nor its fuzzy extension are equipped for representing situations of partial or complete ignorance. To this end, in the Boolean case, we need to introduce a proper representation of partial uncertainty including ignorance in the relational table of the formal context. One may think of introducing gradations of uncertainty by changing crosses and blanks in the table into probability degrees, or by possibility or necessity degrees. In the probabilistic case, one number shall assess the probability that a considered property holds for a given object (its complement to 1 corresponding to the probability it does not hold). However, this is assuming precise knowledge on the probability values, which is not really appropriate if we have to model the state of complete ignorance. It is why we investigate the use of the possibilistic setting in the following.

In the possibilistic setting, crosses may be replaced by positive degrees of necessity for expressing some certainty that an object satisfies a property. The blanks could be refined by possibility degrees less than 1, expressing that it is little possible that an object satisfies a property. However, this convention using a single number in the unit interval for each entry in the context may be misleading as when the number replaces a blank or a cross, the meaning of the number is not the same.

In the possibilistic setting, possibility and necessity functions are related by the duality relation N (A) = 1 -Π(A), that holds for any event A, where A denotes the opposite event [START_REF] Dubois | Possibility Theory: An Approach to Computerized Processing of Uncertainty[END_REF]. Then, for entries (x, y) in the table, we use a representation as a pair of necessity degrees (α, 1β) where α = N ((x, y) ∈ R) (resp. 1-β = N ((x, y) ∈ R)) corresponds to the necessity (certainty) that object x has (resp. does not have) property y. Moreover, we should respect the property min(α, 1β) = 0, since min(N (A), N (A)) = 0 in agreement with complete ignorance, in which case nothing (i.e., neither A nor A) is even somewhat certain. Pairs (1,0) and (0,1) correspond to completely informed situations where it is known that object x has, respectively does not have, property y. The pair (0,0) reflects total ignorance, whereas pairs (α, 1β) s.t. 1 > max(α, 1β) > 0 correspond to partial ignorance.

An uncertain formal context is thus represented by

R U = {(α(x, y), 1 -β(x, y)) | x ∈ O, y ∈ P} where α(x, y) ∈ [0, 1] , β(x, y) ∈ [0, 1].
A relational database with fuzzily-known attribute values is theoretically equivalent to the fuzzy set of all ordinary databases corresponding to the different possible ways of completing the information consistently with the fuzzy restrictions on the attribute values. So, an uncertain formal context may be viewed as a weighted family of all standard formal contexts obtained by changing uncertain entries into sure ones. More precisely, one may consider all the completions of an uncertain formal context. This is done by substituting entries (x, y) that are uncertain, i.e., such that 1 > max(α(x, y), 1β(x, y)) by a pair (1,0), or a pair (0,1). Replacing (α(x, y), 1β(x, y)) by (1, 0) is possible at degree β(x, y), the possibility that x has property y. Similarly, replacing (α(x, y), 1β(x, y)) by (0, 1) is possible at degree 1α(x, y), the possibility that x does not have the property y. In this way, one may determine to what extent a particular completion (a context C) is possible, by aggregating the possibility degrees associated with each completed entry (using min operator). Formally, one can write π(C) = min(min (x,y):(x,y)∈C β(x, y), min (x,y):(x,y) ∈C 1α(x, y)).

Likewise the degree of possibility that (

X, Y ) is a formal concept of R U is π(X, Y ) = sup{π(C) : C such that (X, Y ) is a formal concept of C}.
Useful completions are those where partial certainty becomes full certainty. Indeed, given an uncertain formal context and a threshold pair (u, v), let us replace all entries of the form (α, 0) such that α u with (1, 0) and entries of the form (0, 1β) such that 1β v with (0, 1). All such replacements have possibility 1 according to the above formula. Remaining entries, which are more uncertain, can be systematically substituted either by (1,0), or by (0,1). Considering, the two extreme cases where all such entries are changed into (1,0) and the case when where all such entries are changed into (0,1) gives birth to upper and lower completions, respectively. In this way, two classical (Boolean) formal contexts, denoted R * (u,v) and R * (u,v) are obtained as respective results of the two completions. They allow to determine, for a given threshold (u, v), maximal extensions (resp. minimal intensions) and minimal extensions (resp. maximal intensions) of uncertain formal concepts. It is clear that R * (u,v) ⊆ R * (u,v) . Let us illustrate the idea with an example. Example 3. Table 4 exhibits a formal context where some entries are pervaded with uncertainty. Let us examine the situation regarding formal concepts. Take u = 0.7, v = 0.5 for instance. In context R * (0.7,0.5) , examples of formal concepts are pairs ({6, 7, 8}, {c, d, e}), or ({5, 6, 7, 8}, {d, e}), or ({2, 3, 4}, {g, h}), although with u = 0.9, the last formal concept would reduce to ({2, 3}, {g, h}), i.e., the extent of the concept is smaller. Now consider R * (0.7,0.5) , where the entries with low certainty levels (either in favor or against the existence of the link between x and y) are turned into positive links. Then, ({2, 3, 4}, {g, h}) remains unchanged as a formal concept, while a larger concept now emerges, namely ({5, 6, 7, 8}, {c, d, e}). However, one may prefer to consider the results obtained from R * (0.7,0.5) , where only the almost certain information is changed into positive links. In the example, if we move down u to 0.5, and use R * (0.5,0.5) we still validate the larger former concept ({5, 6, 7, 8}, {c, d, e}). This illustrates the fact that becoming less and less demanding on the level of certainty, may enable the fusion of close concepts (here ({6, 7, 8}, {c, d, e}), and ({5, 6, 7, 8}, {d, e}), providing a more synthetic view of the formal context. This small example is intended to illustrate several points. First of all, it should be clear that being uncertain about the existence of a link between an object and a property is not the same as being certain about a gradual link. Second, under uncertainty, there are formal concepts whose boundaries are not affected by uncertainty, while others are. Lastly, regarding certain enough pieces of information as fully certain may help simplifying the analysis of the formal context. Besides, the proposed setting may also handle inconsistent information by relaxing the constraint min(α, 1β) = 0. This would amount to introducing paraconsistent links between objects and properties. -The parallel of FCA with PoTh leading to the introduction of new operators extends to conceptual pattern structures [START_REF] Ferré | A logical generalization of formal concept analysis[END_REF][START_REF] Ganter | Pattern structures and their projections[END_REF], where the description ∂(x) of an object x, may, e.g., be a possibilistic knowledge base [START_REF] Assaghir | A possibility theory-oriented discussion of conceptual pattern structures[END_REF]; -Applications of FCA to the fusion of conflicting pieces of information issued from multiple sources using pattern structures for labeling sets of possible values in terms of sources supporting them [START_REF] Assaghir | Numerical information fusion: lattice of answers with supporting arguments[END_REF]; -The clustering of sets of objects on the basis of approximate concepts [START_REF] Dubois | Possibility theory and formal concept analysis: characterizing independent sub-contexts[END_REF][START_REF] Gaume | Clustering bipartite graphs in terms of approximate formal concepts and sub-contexts[END_REF],

with labeling of the clusters [START_REF] Navarro | Clustering sets of objects using conceptsobjects bipartite graphs[END_REF]; -The building of conceptual analogical proportions [START_REF] Miclet | From analogical proportions in lattices to proportional analogies in formal concepts[END_REF] on the basis of the formal definition of analogical proportions in non-distributive lattices [START_REF] Miclet | Looking for analogical properties in a formal concept analysis setting[END_REF], conceptualization and analogical reasoning being two basic cognitive activities [START_REF] Barbot | Analogical Proportions and the Factorization of Information in Distributive Lattices[END_REF].

Fig. 1 .

 1 Fig. 1. Square of opposition

Fig. 2 .

 2 Fig. 2. Square of oppositions induced by a relation R and a subset S

Fig. 3 .

 3 Fig. 3. Cube of opposition in formal concept analysis

5

  More Lines for Further Research Let us briefly conclude this survey of works in FCA inspired by PoTh by mentioning other examples of lines of interest for further research:

  Example 1. The table below presents a formal context. Pairs ({6, 7, 8}, {c, d, e}), or ({5, 6, 7, 8}, {d, e}), or ({2, 3, 4}, {g, h}) are examples of formal concepts, while ({5, 6, 7, 8}, {a, b, c, d, e}), ({2, 3, 4}, ({f, g, h}), ({1}, {i}) are minimal subcontexts.

Table 1 .

 1 Many-valued relation

	R 1	Pierre Sophie Mike Nahla
	Age	22	28	30	22
	Salary 1100 1300	1500 1500

Table 2 .

 2 Context subsumption

	R 2	Pierre Sophie Mike Nahla
	age ≥ 20	×	×	×	×
	age ≥ 25		×	×
	age ≥ 30			×
	salaire ≥ 1000 ×	×	×	×
	salaire ≥ 1200		×	×	×
	salaire ≥ 1400			×	×

Table 3 .

 3 Context summarization

	R 3	Pierre Sophie Mike Nahla
	age 'young' 1	0.7	0.6	1
	salary 'low' 1	0.8	0.6	0.6

Table 4 .

 4 Uncertain formal concepts

	1 2	3 4	5	6 7 8
	a				×
	b			×	×
	c			(0.5,0) × × ×
	d			×	× × ×
	e			×	× × ×
	f	(0, 0.8)	×	(0, 0.3)
	g	×	× (0.8, 0)
	h	×	× (0.8, 0)
	i ×


I.-C Morȃrescu

W Michiels
email: wim.michiels@cs.kuleuven.be

M Jungers
email: marc.jungers@univ-lorraine.fr

Effect of a distributed delay on relative stability of diffusely coupled systems, with application to synchronized equilibria

Keywords: Nonlinear oscillators, gamma distributed delays, interconnected systems

We perform a linear stability analysis of synchronized equilibria in networks of diffusely coupled oscillators, affected by distributed delays in the coupling and we characterize the structure of the emanating solutions in the bifurcations, under the assumption that the delay kernels are equal to each other. The motivation comes from the fact that valuable quantitative and qualitative information about the occurrence and type of synchronous or partially synchronous solutions can be obtained from this inear analysis. We analyze stability of synchronized equilibria as a function of the parameters of a so-called shifted gamma-distributed delay, which allows to represent or approximate a large class of distributed delay kernels. We also present an asymptotic analysis method, which is particularly suitable for studying the effect of a distribution of delays on stability. We apply the methods to networks of coupled Lorenz systems, where we highlight that stability is favored by a distributed delay compared with a discrete delay with the same average value. Among others, we show that if the coupling is diffusive, the synchronized equilibria become asymptotically stable for large values of the coupling gain, i.e. the systems locally asymptotically synchronize at an equilibrium, independently of the network topology.

INTRODUCTION

Systems of coupled nonlinear oscillators have been used to describe a variety of phenomena in different fields of science, ranging from electric network modeling to brain activity behavior. The emergence of synchronization in networks of oscillators with couplings given by various types of topologies has received much attention [START_REF] Kuramoto | Chemical Oscillations, Waves and Turbulence[END_REF][START_REF] Pecora | Master stability functions for synchronized coupled systems[END_REF][START_REF] Strogatz | Sync: the emerging science of spontaneous order[END_REF][START_REF] Gfeller | Spectral coarse graining and synchronization in oscillator networks[END_REF]. The influence of discrete fixed delays on the consensus and synchronization of coupled subsystems is considered in [START_REF] Oguchi | Synchronization in networks of chaotic systems with time-delay coupling[END_REF][START_REF] Michiels | Synchronization of delay-coupled nonlinear oscillators: An approach based on the stability analysis of synchronized equilibria[END_REF][START_REF] Morȃrescu | Consensus with constrained convergence rate and time-delays[END_REF][START_REF] Panchuk | Synchronization of coupled neural oscillators with heterogeneous delays[END_REF]. However, as pointed out in several studies [START_REF] Macdonald | Biological delay systems: linear stability theory[END_REF][START_REF] Insperger | Semi-discretization of time-dely systems[END_REF][START_REF] Michiels | Consensus problems with distributed delays, with application to traffic flow models[END_REF][START_REF] Kyrychko | Amplitude death in systems of coupled oscillators with distributeddelay coupling[END_REF], the behavior of interconnected systems is often more accurately reflected by considering distributed delays affecting their interactions. Note in this context that fast varying delays in communication channels can be well approximated by distributed delays, whose kernel is fitted to the probability density function of the delay [START_REF] Michiels | Stabilization of time-delay systems with a controlled, timevarying delay and applications[END_REF].

This paper considers the problem of relative stability in a network of nonlinear oscillators, namely the synchronization problem, where the coupling is diffusive and affected by distributed delays. Precisely we are performing a local stability analysis concerning the linearization of the system around synchronized equilibria. The aim is to characterize the occurrence of synchronous behavior as a function of the parameters of the distributed delay. Preliminary results have been reported in [START_REF] Morȃrescu | Synchronization of coupled nonlinearoscillators with shifted gamma-distributed delays[END_REF]. Here, we are going further by emphasizing the delay distribution effect on stability of the linearized system as well as the information we can get for the nonlinear one. We follow the approach of [START_REF] Michiels | Synchronization of delay-coupled nonlinear oscillators: An approach based on the stability analysis of synchronized equilibria[END_REF][START_REF] Unal | Prediction of partial synchronization in delay-coupled nonlinear oscillators, with application to hindmarsh-rose neurons[END_REF], which is based on analyzing stability and bifurcations of synchronized equilibria. Even though the analysis is local it has been demonstrated in [START_REF] Unal | Prediction of partial synchronization in delay-coupled nonlinear oscillators, with application to hindmarsh-rose neurons[END_REF] that by combining the stability charts and the information on emanating solutions in bifurcations with a characterization of invariant manifolds, valuable qualitative and quantitative information about the occurrence and type of synchronous or partially synchronous motion can be obtained, even if the attractor under consideration is more complicated than an equilbrium (e.g., a chaotic attractor).

The contributions of the paper are two-fold. First, we will provide systematic tools for analyzing the stability and bifurcations of synchronized equilibria when the delays are modeled by a gammadistribution kernel with a gap. This tool is based on the linearized system around the corresponding synchronized equilibria. By an appropriate choice of parameters such a kernel allows to model a broad class of distributed delays [START_REF] Morȃrescu | Stability crossing curves of shifted gamma-distributed delay systems[END_REF]. The applicability will be demonstrated with a case-study on coupled Lorenz systems. Second, we will pay particular attention to analyzing the effect of a delay distribution on synchronized equilibria, in comparison to the case of a discrete ( point-wise) delay, where we use the computed stability crossing curves as well as the asymptotic methods of [START_REF] Michiels | Stabilization of time-delay systems with a controlled, timevarying delay and applications[END_REF].

Intuitively the stabilizing effect of a distributed delay might be understood from the fact that, in contrast to a discrete delay, it acts as a low-pass filter, and, hence, it suppresses the propagation of high frequency signals which might otherwise be destabilizing. In [START_REF] Michiels | Stabilization of time-delay systems with a controlled, timevarying delay and applications[END_REF] a systematic study of the stabilizing effect of fast varying delays on stability is investigated, based on a distributed delay comparison system obtained by averaging. In [START_REF] Sipahi | Stability of traffic flow behavior with distributed delays modeling the memory effects of the drivers[END_REF] distributed delays are used in modeling the driver's reaction in a car following system, and compared with existing models with discrete delay. In [START_REF] Atay | Distributed delays facilitate amplitude death of coupled oscillators[END_REF] it is illustrated that distributed delay facilitates amplitude death of coupled oscillators, meaning that stability islands in parameter spaces of equilibria merge in such a way that the oscillations disappear. At the end of this paper we will analyze a form of amplitude death in networks of coupled Lorenz systems, which exhibit chaotic behavior without coupling, but, for sufficiently large coupling gain, synchronize at an equilibrium. Finally, the potentially destabilizing effect of approximating distributed delays by discrete delays is important in developing robust prediction based controllers for systems with input and output delays, as it has been shown that approximating the distributed delays in the control law, which originate from the inclusion of finite impulse response (FIR) filters, by discrete delays using a numerical quadrature rule, might destabilize the () control loop [START_REF] Van Assche | Implementation of a distributed control law for a class of systems with delay[END_REF]. The instability mechanism, as well as remedies, are explained in, e.g., [START_REF] Mondié | Finite spectrum assignment of unstable time-delay systems with a safe implementation[END_REF] and the references therein.

The structure of the paper is as follows. In Section 2 we describe the interconnected system and outline the main assumptions. In Section 3 we perform a linear stability analysis of synchronized equilibria as a function of the parameters of a shifted γ-distributed delay. In Section 4 we analyze the beneficial effect of distributed delays using an asymptotic method. Finally, in Section 5 the application to coupled Lorenz systems is considered.

PROBLEM FORMULATION AND PRELIMINARIES

We consider a set of p identical oscillators coupled via a fixed, directed network topology. Each oscillator has a behavior described by the following nonlinear dynamics

ẋi (t) = f (x i (t)) + Bu i (t), y i (t) = Cx i (t), i = 1, . . . , p (1)
where x i ∈ R nx , i = 1, . . . , p, B, C ∈ R nx×m , and f : R n → R n is twice continuously differentiable.

Assumption 1

We assume that for u i = 0 the system has at least one unstable equilibrium of focus type, i.e., the Jacobian has two complex conjugate eigenvalues in the open right half plane, which we denote by

x * in what follows.

The network topology is described by the directed graph

G(V, E, G) (2)
characterized by the node set V = {1, . . . , p}, a set of edges E and a weighted adjacency matrix G defined by the entries α k,l . The edge (k, l) ∈ E exists if and only if α k,l = 0. Moreover, we suppose that self-loops are not present in the graph, i.e. α l,l = 0, ∀l ∈ V . Each vertex corresponds to an oscillator and the edges describe the coupling strength between them. We assume a diffusive coupling between the oscillators, where the information sent from other systems is affected by a distributed delay. More precisely, we couple the systems (1) by means of the control law

u i (t) = k p l=1 α i,l ∞ 0 g(θ)(y l (t -θ) -y i (t))dθ , i = 1, . . . , p, (3)
where k > 0 represents the gain parameter and g(•) stands for the γ-distributed delay with a gap kernel given by

g(θ) =    0, θ < τ, (θ-τ ) n-1 e -θ-τ T T n (n-1)! , ξ ≥ τ, (4) ()
where n ∈ N, T > 0 and τ ≥ 0. Note that g(θ) ≥ 0 for all θ ≥ 0 and ∞ 0 g(θ)dθ = 1. The average delay satisfies

τ av := ∞ 0 g(θ)θdθ = τ + nT.
A motivation for the consideration of this type of delay distribution can be found in [START_REF] Morȃrescu | Stability crossing curves of shifted gamma-distributed delay systems[END_REF].

The coupled system (1)-( 3) has an equilibrium

(x * , . . . , x * ) ∈ R n×p , (5)
which we denote the synchronized equilibrium.

Assumption 2

The graph G is strongly connected.

Assumption 3

The adjacency matrix G is stochastic i.e.

p l=1 α i,l = 1, i = 1, . . . , p.
The first assumption is natural in the context of synchronization while the second assumption will be motivated in what follows. The second assumption ensures that synchronous behavior is preserved if it characterizes the past, i.e., it guarantees the presence of a synchronization manifold.

We recall that Perron-Frobenius theorem for stochastic matrices assures that 1. G has a simple eigenvalue equal to 1 and [1, . . . , 1] is the corresponding right-eigenvector.

2. All eigenvalues of G have modulus smaller than or equal to 1.

In what follows we denote the eigenvalues of G as λ i (G), i = 1, . . . , p, where we take the following convention that λ 1 (G) = 1.

Following the coordinate transformation in [START_REF] Michiels | Synchronization of delay-coupled nonlinear oscillators: An approach based on the stability analysis of synchronized equilibria[END_REF], we can bring the system (1) and (3) in the form

ẋ1 (t) = f (x 1 (t)) + BkC p l=1 α 1,l ∞ 0 g(θ)(x 1 (t -θ) -x 1 (t))dθ + BkC p l=1 α 1,l ∞ 0 g(θ)e l (t -θ)dθ, (6)
    ė2 (t)
. . .

ėp(t)     =     f (x 1 (t) + e 2 (t)) -f (x 1 (t))
. . .

f (x 1 (t) + ep(t)) -f (x 1 (t))     -k                   p l=1 α 2,l . . . p l=1 α p,l          ⊗ BC              e 2 (t)
. . .

ep(t)     + k G ⊗ BC     ∞ 0 g(θ)e 2 (t -θ)dθ . . . ∞ 0 g(θ)ep(t -θ)dθ     -          k          p l=1 α 1,l - p l=1 α 2,l . . . p l=1 α 1,l - p l=1 α p,l          ⊗ BC          ∞ 0 g(θ)(x 1 (t -θ) -x 1 (t))dθ, (7) ()
where e i = x i -x 1 , i = 2, . . . , p and σ( G) = σ(G) \ {1}. Assumption 3 guarantees the preservation of synchronized solutions (that is, e 2 (s) = . . . = e p (s) = 0, ∀s ≤ t implies ė2 (t) = . . . ėp (t) = 0). The forward invariant set, described by

x 1 = • • • = x p is called the synchronization manifold.
Note from (6) that the coupling is invasive, in the sense that the dynamics on the synchronization manifold are affected by the coupling. This property is inherent to diffusive coupling affected by transmission delays.

SHIFTED γ-DISTRIBUTED DELAYS: STABILITY ANALYSIS OF SYNCHRONIZED EQUILIBRIA

We perform a stability analysis of the equilibrium (5) in the parameter space (k, T, τ ) and characterize the emanating solutions in the bifurcations.

Linearizing the system (1) and (3) around the synchronized equilibrium (x * , . . . , x * ), we obtain

    v1 (t) . . . vp (t)     = I ⊗ (A -BkC)     v 1 (t) . . . v p (t)     + kG ⊗ BC     ∞ 0 g(θ)v 1 (t -θ)dθ . . . ∞ 0 g(θ)v p (t -θ)dθ     , (8)
where A = ∂f ∂x (x * ).

Characteristic equation

The characteristic function of ( 8) is given by

∆(s, k, T, τ ) = det F (s, k, T, τ ), (9)
where the characteristic matrix F is defined by

F (s, k, T, τ ) = I ⊗ (sI -A + kBC) -G ⊗ kBC e -sτ (1 + sT ) n . ( 10
)
Using the factorization G = P ΛP -1 , where Λ ∈ C p×p is triangular and P ∈ C p×p , the characteristic function may be decomposed as

∆(s, k, T, τ ) = p i=1 ∆ i (s, k, T, τ ), (11)
where

∆ i (s, k, T, τ ) = det F i (s, k, T, τ ), F i (s, k, T, τ ) = sI -A + kBC -kBCλ i (G) e -sτ (1 + sT ) n , i = 1, . . . , p. (12)
()

Eigenspaces and behavior on the onset of instability

We investigate the eigenspace of the characteristic matrix [START_REF] Michiels | Consensus problems with distributed delays, with application to traffic flow models[END_REF], corresponding to a characteristic root. The presence of eigenvalues of G of multiplicity larger than one, leads to discussions related to the dimension of the eigenspaces and significantly complexifies the characterization of solutions of the linearized system. For the sake of clarity and in order to keep the reader focused on the main objective of the section, we simplify the analysis by restricting ourselves to the generic case where all the eigenvalues of the adjacency matrix G are simple. Let E i be the eigenvector of G corresponding to the eigenvalue λ i (G), i = 1, . . . , p.

If for some l ∈ {1, . . . , p}, the equation ∆ l (s, k, T, τ ) = 0 has a simple root at s = ŝ such that

F l (ŝ, k, T, τ ) V = 0, V ∈ C n×1 , (13)
then it can be verified that

F (ŝ, k, T, τ ) (E l ⊗ V ) = 0. ( 14
)
This implies that the linearized system ( 8) has an exponential solution

    ν 1 (t) . . . ν p (t)     = c(E l ⊗ V )e λt = c     E l,1 V . . . E l,p V     e λt , (15)
with the constant c depending on the initial conditions. In words, in an exponential solution of ( 8), corresponding to a zero of ∆ l (s, k, T, τ ), the relation between the state variables of an individual subsystem is determined by the null vector of F l , while the relation between the corresponding state variables of the different subsystems is solely determined by the eigenvector E l , corresponding to the l-th eigenvalue of the adjacency matrix G. This implies that all modes can be classified in at most p types, based on the relations between the behavior of the different subsystems. Note for example that bifurcations of the synchronized equilibrium due to imaginary axis zeros of f 1 are bifurcations where the synchronized behavior is preserved in the emanating solution, as

E 1 = [1 1 • • • 1] T .
In [START_REF] Unal | Prediction of partial synchronization in delay-coupled nonlinear oscillators, with application to hindmarsh-rose neurons[END_REF], for the case of discrete delay, the above information is used to predict the occurrence of synchronous or partially synchronous regimes as a function of parameters of the coupling, even if the attractor is different from an equilibrium. Thereby it complements the conventional, Lyapunov based analysis. To explain the main idea, let us assume that we have p = 4 agents and that all unstable characteristic roots are due to zeros of f 1 and

f k , k ∈ {2, 3, 4} such that E k = [1 1 -1 -1] T .
In addition, assume that the set M := x ∈ C 4nx : x 1 = x 2 ; x 3 = x 4 is forward invariant with respect to the dynamics of the coupled system, i.e., it is a so-called partial synchronization manifold. From (15), we then conclude that close to the synchronized equilibrium the solutions are repelled but in such a way that the synchrony between agents 1 and 2, and the synchrony between agents 3 and 4 are maintained. This preservation of partial synchrony might be an indication (but no proof) that the manifold M is stable, and the corresponding partially synchronous motion occurs.. See further Figure 2 and Figure 6. The special structure of the unstable eigenvectors and in particular the "Hopf 90" bifurcations in Figure 2 are an indication (but not more than that) of stable () motion characterized by a 90 degree phase shift. Such a behavior is effectively observed in Figure 6.

Computation of stability crossing curves in parameter spaces

Due to the specific form of the characteristic matrix [START_REF] Michiels | Consensus problems with distributed delays, with application to traffic flow models[END_REF], which can be pre-multiplied with (1 + sT ) n , standard methods for the stability and bifurcation analysis of systems with discrete delays can be used. In Section 5 we will use the package DDE-BIFTOOL [START_REF] Engelborghs | Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL[END_REF] in order to compute stability determining characteristic roots corresponding to the synchronized equilibria. Stability regions in the (T, τ ) and (k, τ ) parameter spaces will be computed by the package as Hopf bifurcation curves of the nonlinear coupled system. Here the computational cost can be significantly reduced by taking into account the factorization [START_REF] Michiels | Stability and Stabilization of Time-Delay Systems: An Eigenvalue Based Approach[END_REF] of the characteristic matrix, as well as the following property.

Proposition 1

The following invariance properties of imaginary axis characteristic roots hold:

1. For i ∈ {1, . . . , p}, the property ∆ i (jω, k, T, τ ) = 0 for ω > 0 implies that ∆ i jω, k, T, τ + 2πl ω = 0, where l = 1, 2, . . ..

For

i, l ∈ {1, . . . , p}, such that |λ i (G)| = |λ l (G)|, and ω > 0, the equality ∆ i (jω, k, T, τ ) = 0 implies ∆ l (jω, k, T, τ ), where τ = τ + ∠(λ l (G))-∠(λi(G)) ω .
For computing stability crossing curves specifically in the (T, τ ) parameter space, a semianalytical method can also be used, which extends the approach for single input-single output systems of [START_REF] Morȃrescu | Stability crossing curves of shifted gamma-distributed delay systems[END_REF] . More precisely, defining the frequency crossing set Ω as the set collection of all ω > 0 such that there exists a parameter pair (T, τ ) satisfying ∆(jω, k, T, τ ) = 0, we have the following result.

Proposition 2

For a given gain k * and any ω ∈ Ω, all corresponding (T, τ ) values satisfying ∆ i (jω, k * , T, τ ) = 0 can be computed by

T = T l = 1 ω | z l (ω) | -2/n -1 1/2 (16)
τ = τ l,m = 1 ω (-∠z l (ω) -n arctan(ωT ) + 2mπ), m = 0, ±1, ±2, .... ( 17
)
where z l (ω), l = 1, . . . , n ω , satisfy | z l (ω) |< 1 and are eigenvalues of the pencil

Λ(z) = (-jωI + A -k * BC) + zλ i (G)(k * BC).

Proof

Let us consider the following variable transformation

z = e -jωτ (1 + jωT ) n . () Then, the equation ∆ i (jω, k * , T, τ ) = 0 rewrites as det (jωI -A + kBC -kBCλ i (G)z) = 0.
From the definition of z it is clear that the admissible solution satisfy |z| < 1. For such a solution z l , the corresponding pair (T, τ ) can be simply computed by [START_REF] Morȃrescu | Consensus with constrained convergence rate and time-delays[END_REF] and [START_REF] Morȃrescu | Stability crossing curves of shifted gamma-distributed delay systems[END_REF].

Proposition 2 gives rise to a frequency sweeping procedure for computing stability crossing curves in the (T, τ ) parameter space, where as a function of ω the eigenvalues inside the unit disk of pencil Λ(z) need to be computed. Note that, as the factor λ i (G) can be absorbed into z, only one eigenvalue computation per frequency is necessary. For a characterization of the crossing directions of imaginary axis roots, we refer to [START_REF] Morȃrescu | Stability crossing curves of shifted gamma-distributed delay systems[END_REF] as the results directly carry over.

The following result characterizes stability for large values of the gain parameter.

Theorem 1

Consider the coupled systems (1) and ( 3). Let T > 0 be given. Assume that for τ = 0 the synchronized equilibrium is asymptotically stable for sufficiently large k. Then there is a constant k such that for every k > k, the synchronized equilibrium is asymptotically stable independently of τ ≥ 0.

Proof

We show that for large values of k eigenvalues cannot cross the imaginary axis when varying parameter τ . The theorem then follows from the assumed stability for τ = 0.

The equation

∆ i (jω, k, T, τ ) = 0 is equivalent to det I -(jωI -A + kBC) -1 kBCλ i (G) e -jωτ (1 + jωT ) n = 0.
A necessary solvability condition is given by

ρ (jωI -A + kBC) -1 kBC λ i (G) (1 + jωT ) n = 1,
with ρ(•) the spectral radius. This condition is always violated for large k. Indeed, in the complex plane the nonzero eigenvalues of the matrix

j ω k I - A k + BC -1 BCλ i (G), (18)
converge to the curve

Ω ≥ 0 → 1 jΩ + 1 λ i (G)
as k → ∞, uniformly in the parameter ω ≥ 0. Furthermore, we have:

1 1 + jΩ λ i (G) ≤ |λ i (G)| ≤ 1, ∀Ω ≥ 0 () and 1 |(1 + jωT ) n | < 1, ∀ω > 0.
It follows that zeros on the imaginary axis are not possible for large values of k.

ASYMPTOTIC ANALYSIS OF THE DELAY DISTRIBUTION EFFECT ON STABILITY

In Section 3 we outlined how to compute stability information as a function of the parameters n, T of the shifted γ-distributed delay. In this section we analyze the effect on stability of a transition from a discrete delay to a distributed delay. Here we use an asymptotic method originally propesed in [START_REF] Michiels | Stabilization of time-delay systems with a controlled, timevarying delay and applications[END_REF] and later rediscovered in [START_REF] Atay | Delayed feedback control near hopf bifurcation[END_REF]. The obtained results are asymptotic, yet independent of the shape of the distributed delay.

The kernel of the distributed delay (4) can be expressed as

g(θ) = w(θ -τ av ) (19)
where the function w : R → R satisfies the following properties,

w(θ) ≥ 0, ∀θ ∈ R, w(θ) = 0, ∀θ ≤ -τ av , ∞ -∞ w(θ)dθ = 1, ∞ -∞ θw(θ)dθ = 0. (20)
The characteristic equation corresponding to the synchronized equilibrium can now be written as

Π p i=1 det sI -A + kBC -kBCλ i (G)e -sτav η(s) = 0, where η(s) = ∞ -τav
w(θ)e -sθ dθ.

In this way we have, in the frequency domain description, separated the contribution of the delay in a factor exp(-sτ av ), related to the average delay, and a factor η(s), related to distribution of the delay around this average value. Note that, specifically for the kernel (4), we have

w(θ) =    0, θ < -nT, (θ+nT ) n-1 e -θ+nT T T n (n-1)! , θ ≥ -nT, η(s) = e nT s (1+sT ) n ⇒ e -sτav η(s) = e -sτ (1+sT ) n ,
which is in agreement with the results in the previous section. Note further that the function s → η(s) is entire and satisfies

η(0) = 1, η (0) = 0, η (0) > 0. (21)
Let us now, instead of ( 4) and [START_REF] Panchuk | Synchronization of coupled neural oscillators with heterogeneous delays[END_REF], consider the scaled kernel

g δ (θ) = 1 δ w θ -τ av δ (22)
() where δ ∈ (0, 1]. In words, δ is a scaling factor describing the distribution of the delay around the average value. The smaller δ, the more the delay is concentrated around the average value, and in the limit δ → 0+, function g δ becomes a Dirac impulse, corresponding to a discrete delay τ av . This is clarified in Figure 1.

0
With the scaled kernel, the characteristic equation becomes

h(s, η(sδ)) = 0, (23) where h(s, η(sδ)
) := Π p i=1 det sI -A + kBC -kBCλ i (G)e -sτav η(sδ) . (24)
Consistent with the above observation on the limit δ → 0+, the contribution of the delay in [START_REF] Unal | Prediction of partial synchronization in delay-coupled nonlinear oscillators, with application to hindmarsh-rose neurons[END_REF] reduces to e -sτav for δ = 0.

Consider an isolated characteristic root s 0 for δ = 0, i.e., for the case of discrete delay. From (24) we can compute the sensitivity of the root with respect to parameter δ at δ = 0. This sensitivity describes the initial effect on the root of distributing the delay around its original version. More precisely, by the implicit function theorem, Equation (24) describes a function δ ∈ [0, δ m ] → r(δ) for some δ m > 0 satisfying r(0) = s 0 and h(r(δ), η(r(δ)δ)) = 0 for all δ ∈ [0, δ m ). An application of the chain rule, taking into account [START_REF] Sipahi | Stability of traffic flow behavior with distributed delays modeling the memory effects of the drivers[END_REF], yields

   r (0) = 0, r (0) = - ∂h ∂η (s0,1) ∂h ∂s (s0,1) η (0)s 2 0 .
A corollary related to the stabilizing effect of a delay-distribution is the following, where we used g (0) > 0.

Proposition 3

Assume that for δ = 0, the rightmost characteristic roots of the synchronized equilibrium of the coupled system are simple and on the imaginary axis. Denote them by jω i , i = 1, . . . , m. If ∂h ∂η (jω i , 1) ∂h ∂s (jω i , 1)

< 0 i = 1, . . . , m, (25) ()
then the synchronized equilibrium is asymptotically stable for small values of δ. If for some i ∈ {1, . . . , m} the inequality is reversed, the synchronized equilibrium is unstable for small values of δ.

It is important to note that in the derivation of the above result only the properties (20) of w are used, hence, it is largely independent of the way in which the delay is distributed. Finally, considering the behavior of a characteristic root s 0 as a function of δ as well as the nominal delay τ av around (δ, τ av ) = (0, τ av,0 ) leads to a root function (δ, τ av ) → r(δ, τ av ). A straightforward computation leads to

∂ 2 r(0, τ av ) ∂δ 2 < 0 (> 0) ⇔ ∂r(0, τ av ) ∂τ av < 0 (> 0), (26)
hence, the sensitivity of a root with respect to a distribution of the delay can be obtained as a byproduct from its sensitivity with respect to a change of size of the (average) delay.

In [START_REF] Michiels | Stabilization of time-delay systems with a controlled, timevarying delay and applications[END_REF] the effect of fast time-varying discrete delays on stability of (stand-alone) control system was analyzed. In the first step the effect of a fast varying discrete delay and a distributed delay were correlated via an averaging based argument. In the second step the distributed delay was analyzed using the above asymptotic approach. The three examples in [START_REF] Michiels | Stabilization of time-delay systems with a controlled, timevarying delay and applications[END_REF], including a model for a high speed cutting process, demonstrate the stabilizing effect of a fast variation and distribution of the delay.

The qualitative properties of the spectrum of time-delay systems, described in Chapter 1 of [START_REF] Michiels | Stability and Stabilization of Time-Delay Systems: An Eigenvalue Based Approach[END_REF], give an intuitive argument: in our intensive experiments with characteristic roots computations of delay equations it was often observed (but not always) that roots in the upper half plane follow the exponential envelope curves containing the spectrum and described by Proposition 1.10 of [START_REF] Michiels | Stability and Stabilization of Time-Delay Systems: An Eigenvalue Based Approach[END_REF], that is, they cross the imaginary axis from left to right and with decreasing imaginary part. An example where this statement can actually be proven is the first order system ẋ(t) = -kx(t -τ ), k > 0, based on the analytical computation of the crossing direction (using, e.g., the forumulas of Chapter 4 of [START_REF] Michiels | Stability and Stabilization of Time-Delay Systems: An Eigenvalue Based Approach[END_REF]). Example 2 in the article [START_REF] Lichtner | The spectrum of delay differential equations with large delay[END_REF] on the asymptotic behavior of roots for large delay illustrates a much more complex scenario where multiple critical frequencies and crossings in both directions occur (see Figure 3.1 in the reference), characterized by decreasing imaginary part at the crossings as well. If this crossing behavior happens, a delay distribution has a stabilizing effect, by relation (26). In the next section we illustrate the stabilizing effect for networks of coupled Lorenz systems, where we obtain results independent of the network topology.

CASE STUDY: COUPLED LORENZ SYSTEMS

In this section the nonlinear oscillators (1) are specified as Lorenz systems:

     ẋi,1 (t) = σ(x i,2 (t) -x i,1 (t)), ẋi,2 (t) = rx i,1 (t) -x i,2 (t) -x i,1 (t)x i,3 (t) + u i,1 (t), ẋi,3 (t) = -bx i,3 (t) + x i,1 (t)x i,2 (t) + u i,2 (t), y i,1 (t) = x i,2 (t), y i,2 (t) = x i,3 (t) -r, i = 1, . . . p, ( 27
) ()
where (28)

u i = [u i,1 u i,2 ] T , y i = [y i,1 y i,2 ] T .
Note that for u i ≡ 0 each Lorenz system has three equilibria given by

(0, 0, 0), ± b(r -1), ± b(r -1), r -1 , (29)
the latter two corresponding to unstable foci. Thus, Assumption 1 holds. Furthermore, with the parameter values (28) it exhibits a chaotic attractor [START_REF] Strogatz | Nonlinear dynamics and chaos: with applications to physics, biology, chemistry and engineering[END_REF].

If we linearize the coupled system ( 27) and ( 3) around the synchronized equilibrium

(x * , . . . , x * ), x * = (± b(r -1), ± b(r -1), r -1), ( 30
)
then we obtain the linear system [START_REF] Lichtner | The spectrum of delay differential equations with large delay[END_REF], where the matrices are specified as

A =    -σ σ 0 1 1 ∓ b(r -1) ± b(r -1) ± b(r -1) -b    , B = C =    0 0 1 0 0 1    .
It is easy to show that the stability of the linearized system does not depend on which equilibrium x * in (30) is considered, Therefore, we will restrict ourselves to the one in the positive octant.

In what follows we analyze the stability properties of the synchronized equilibria (30) in the (k, T, τ ) parameter space. First we study the asymptotic behavior for large values of the gain parameters in §5.1. For the standard parameters (28) this will allow us to make assertions about stability regions, stability switches and emerging behavior, which do not depend on the network topology. Next we present a numerical example in §5.2.

Behavior for large gain values

For T = 0, i.e., the case of discrete delay, the behavior of the stability crossing curves for large gain values is described by the following result from [START_REF] Michiels | Synchronization of delay-coupled nonlinear oscillators: An approach based on the stability analysis of synchronized equilibria[END_REF].

Theorem 2

Consider a network of coupled Lorenz systems (27) with parameters (28) and coupling [START_REF] Engelborghs | Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL[END_REF]. Assume that T = 0. Then there exists a number k > 0 and a function

τ * : [ k, ∞] → R + , k → τ * (k), (31)
satisfying the following properties:

1. there is a constant k > k such that for every k > k, the synchronized equilibrium has two characteristic roots in the open right half plane for all τ ∈ [0, τ * ], while it is asymptotically stable for τ ∈ (τ * , τ * + ), with sufficiently small; 2. at τ = τ * Hopf bifurcation of the synchronized equilibrium of ( 27) and (3) occurs that is synchrony preserving (for decreasing delay);

() 3. for all k ∈ [ k, ∞] we can factor τ * (k) = υ(k) k where lim k→∞ υ(k) = 0.586004. ( 32
)
Furthermore, the number k and the function (31) are independent of the number of subsystems and of the network topology.

In addition, in [START_REF] Michiels | Synchronization of delay-coupled nonlinear oscillators: An approach based on the stability analysis of synchronized equilibria[END_REF] asymptotic synchronization has been proven that for k sufficiently large and kτ sufficiently small: boundedness of solutions proven using the Lyapunov functional

V = p i=1 γ i {y i (t) T y i ) + k/2 p i=1 y i (θ) T y(θ)dθ},
with [γ 1 . . . γ p ] the left eigenvector of G corresponding to eigenvalue 1, while stability of the synchronization manifold is inferred from a small gain argument. In fact, for k large and kτ sufficiently small synchronized chaotic behavior is observed.

The following theorem shows that the qualitative behavior for large values of k is significantly different if T = 0.

Theorem 3

Consider a network of coupled Lorenz systems (27) with parameters (28) and coupling [START_REF] Engelborghs | Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL[END_REF]. Assume that T > 0. Then there is a constant k such that for every k > k, the synchronized equilibrium is asymptotically stable independently of τ ≥ 0.

Proof

Given Theorem 1, it suffices to show that for τ = 0 all characteristic roots are in left half plane for sufficiently large k.

Let i ∈ {1, . . . , p}. As k → ∞ the function ∆ i (s, k, T, 0) k 2 = det N i (s, k), with N i (s, k) =      s + σ -σ 0 -1 k s+1 k -λi(G) (1+st) n -1 √ b(r-1) k - √ b(r-1) k - √ b(r-1) k s+b k -λi(G) (1+sT ) n -1     
, uniformly converges on compact subsets of C to the function

1 -λi(G) (1+sT ) n 2 (s + σ). (33)
From Rouché's theorem it follows that, as k → ∞, 2n + 1 zeros of ∆ i (s, k, T, 0) converge pairwise to the 2n + 1 zeros of (33) (taking multiplicity into account). Similarly, letting s = s k we get ()

∆ i (sk, k, T, 0) k 3 = det(sI -P i (s, k)),
where

P i (s, k) =      -σ k σ k 0 1 k -1 k + λi(G) (1+skT ) n -1 - √ b(r-1) k √ b(r-1) k √ b(r-1) k -b k + λi(G) (1+skT ) n -1      .
Thus, ∆i(sk,k,T,0)

k 3
uniformly converges on compact sets in the s-plane to s(s + 1) 2 . Rouché's theorem then implies that for sufficiently large k the function ∆ i (s, k, T, 0) has two zeros equal to k.s 1,2 , where

lim k→∞ s1,2 = -1,
hence they move off the infinity without leaving the open left half plane.

The function ∆ i (s, k, T, 0) has exactly 2n + 3 zeros, induced by rank(BC) = 2. We conclude from the previous observations that all these zeros are in the open left half plane for large k if the zeros of (33) are in the open left half plane. This is the case for i ∈ {2, . . . , p}. The case i = 1 deserves special attention, since for λ 1 (G) = 1, (33) has two (rightmost) zeros equal to zero. The latter means that two zeros of f 1 (s, k, T, 0) approach the origin as k → ∞. It remains to show that these approach the imaginary axis from the left.

With ŝ = ks we define g(ŝ, k) := f 1 ŝ k , k, T, 0 . By the expansion

1 (1 + ŝT /k) n -1 = - ŝT n k + O ŝT n k 2
and, again, an application of Rouché's theorem, it follows that, as k → ∞, two zeros of g(•, k) converge to the finite eigenvalues of the pencil -A + ŝ(T n)BC. With the parameters (28) these eigenvalues are given by ŝ1,2 = 1 T n (-1.333 ± 1.926j). We conclude that the two zeros of f 1 (s, k, T, τ ) that converge to the origin as k → ∞ satisfy

s 1,2 = 1 k ŝ1,2 + O 1 k 2 ,
hence, (s 1,2 ) < 0 for large k.

Numerical examples

We consider a ring topology with unidirectional coupling, described by the adjacency matrix

G =       0 • • • 0 1 1 0 . . . . . . 1 0       ∈ R p×p , (34) ()
which has the following properties:

λ l (G) = e j 2π(l-1) p , E l =        1 e -j 2π(l-1) p . . . e -j 2π(p-1)(l-1) p       
for l = 1, . . . , p. If ( 13) is satisfied for ŝ = jω, ω > 0, then the emanating solution [START_REF] Morȃrescu | Synchronization of coupled nonlinearoscillators with shifted gamma-distributed delays[END_REF] becomes

    ν 1 (t) . . . ν p (t)     = c          V e jωt V e jωt-2π(l-1) p V e jωt-4π(l-1) p . . . V e jωt-2(p-1)π(l-1) p          . ( 35
)
It can be interpreted as a traveling wave solution, where the agents follow each other with a phase shift of 360(l -1)/p degrees. Therefore, if the characteristic root ŝ on the imaginary axis corresponds to a Hopf bifurcation of the original nonlinear system (27) and ( 3) for a critical value of some free parameter, we refer to this bifurcation as a "Hopf 360(l -1)/p" bifurcation. For p = 4 subsystems the stability regions in the (k, τ ) parameter space are displayed in Figure 2, for T = 0 (discrete delay). These results are in accordance with Theorem 2. The blue curve segments correspond to the case where condition ( 25) is satisfied, the red segments to the case condition is not satisfied. We conclude that for k > 0.1 the (initial) effect of a delay distribution is stabilizing.

In Figure 3 we depict the corresponding stability crossing curves for a shifted γ-distributed delay with parameters T = 0.01 and n = 1. The stability independent of τ for large k is expected from Theorem 3. The similarity between the different curves is a consequence of the invariance properties in Proposition 1.

To relate the asymptotic results on the stabilizing effect of a delay distribution, indicated in Figure 2 by the colors, and the stability regions for the distributed delay with parameters T = 0.01 and n = 1, we superpose the stability regions shown in Figures 2 and3 on Figure 4. It should be noted that in the horizontal axis we now display the average delay τ av = τ + nT .

In order to validate the above results and illustrate the information that can be obtained from the (local) analysis of synchronized equilibria (see also the discussion at the end of Section 3.2), we performed simulations of the nonlinear coupled system. First, we considered the case of discrete delays, T = 0. We take p = 4, k = 25 and perform a simulation for the values of τ indicated in Figure 2 (left) with crosses, where the solution is initialized with a random perturbation of a synchronized equilibrium (the perturbation on each state variable is obtained from a normal distribution with zero mean and standard deviation one). For τ = 0.008 synchronized chaotic behavior is observed, see the top-left pane in Figure 5. This is consistent with the Hopf 0 bifurcation (synchrony preserved in the unstable manifold) and the asymptotic synchronization proven in [START_REF] Michiels | Synchronization of delay-coupled nonlinear oscillators: An approach based on the stability analysis of synchronized equilibria[END_REF]. For τ = 0.05 the stability of the synchronized equilbrium is confirmed (top-right pane of Figure 5). For τ = 0.15, depicted in the bottom pane of Figure 5, the behavior becomes quite interesting. The system exhibits a chaotic like attractor where the solutions of each subsystem are irregularly 25) is satisfied. Whenever the blue curve bounds a stability region, the latter means that these values of k and τ correspond to an asymptotically stable synchronized equilibrium for small δ > 0 (replacing the discrete delay by a distributed delay), i.e., the stability region grows locally. swung from the neighborhood of one unstable equilibrium to the neighborhood of another unstable equilibrium, as for an uncoupled system. However, when zooming in on the solutions, it turns out that a phase shift of π/2 degrees between the four systems is observed, see Figure 6. This is in agreement with Figure 2 (left): the delay value τ lies in a region where the synchronized equilibria is unstable. The two surrounding Hopf 90 bifurcations curves imply that solutions are repelled from the synchronized equilibria, but with the property that the solutions in the unstable manifold are characterized by a π/2 phase shift. Roughly speaking, solutions are the kicked away, but the kicks on the subsystems occur in a periodic, wavy manner. This illustrates that looking at the structure of eigenvectors at stability switches of synchronized equilibria may provide valuable information about the dynamics also beyond equilibria. Next, we considered the case of distributed delay with n = 1, T = 0.01 and k = 25. Our simulations for different values of τ confirmed the delay-independent stability. See Figure 7 for the simulation for τ = 0.15 (indicated with a cross in the left pane of Figure 3). Finally, Figure 8 shows the stability regions in the (T, τ ) parameter space for k = 12 and k = 25.

() In order to illustrate the applicability to other network configurations we consider a network of 8 Lorenz systems, where the communication graph is described by

G =               
               , (36)
with eigenvalues {1, 0.2624 ± 0.5678j, -0.1676 ± 0.4310j, -0.4698 ± 0.1681j, -0.2500}.

(37)

The coupling delay is described by a shifted γ-distributed delay. In Figure 9 we plot the stability regions of the synchronized equilibria as a function of the parameters k and τ for n = 1 and T = 0.01. Notice, once more, the delay-independent stability for large values of k, a property independently of the network. Compared to Figure 3 the local maxima of the stability crossing curves are not equal to each other, which can be explained by the fact that not all eigenvalues (37) have equal magnitude (see also Proposition 1).

CONCLUSION

We gained qualitative and quantitative insight in the effect of a delay distribution on relative stability in delay-coupled nonlinear oscillations, thereby modeling latency in the coupling by shifted γdistributed delay. This was done by studying the stability regions of the synchronized equilibrium in the (coupling gain-delay) parameter space and by characterizing the structure of the emanating solutions in the bifurcations. In addition, we studied in a semi-analytical way, the effect of a transition from a discrete to a distributed delay, focusing on the effect on (relative) stability.

Applying the results to the particular case of coupled Lorenz systems we have shown that, independently of the network topology, for sufficiently large coupling gains the distribution of the delay has a stabilizing effect on the stability of the synchronized equilibrium. The fact that the synchronization mechanism is a high-gain mechanism, based on dominating the destabilizing terms in the synchronization error dynamics, suggests that a similar conclusion can be drawn for other types of nonlinear oscillators.

The quantitative tools for stability analysis can also be applied to the case where the coupling is described instead by

however, the qualitative results are expected to be different. A major difference between (3) and ( 38) is that coupling (38) is non-invasive, in the sense that it vanishes on the synchronization manifold. As a consequence, the dynamics on this manifold are described by the dynamics of one uncoupled subsystem, and if the latter only has unstable equilibria, the synchronized equilibria of the network can never be stable.

Figure 1 .

 1 Figure 1. The scaled kernel g δ for different values of δ and parameters T = τ = 1, n = 2. For δ = 1 we have g δ = g. If δ is decreased the scaled kernel becomes more concentrated around τav = τ + nT . The limit δ → 0 corresponds to a Dirac impulse, associated with a discrete delay.

  The parameter values are given by σ = 10, r = 28, b = 8/3.

Figure 2 .

 2 Figure 2. Stability regions of synchronized equilibria of 4 coupled Lorenz systems in a ring configuration, for T = 0. The numbers refer to the number of characteristic roots in the right half plane. The Hopf bifurcation curves are classified based on type of emanating solution (35). The blue color indicates that condition (25) is satisfied. Whenever the blue curve bounds a stability region, the latter means that these values of k and τ correspond to an asymptotically stable synchronized equilibrium for small δ > 0 (replacing the discrete delay by a distributed delay), i.e., the stability region grows locally.

Figure 3 .

 3 Figure 3. Stability regions of the synchronized equilibria of 4 coupled Lorenz systems in a ring configuration.Distributed delay with parameters T = 0.01 and n = 1.

Figure 4 .Figure 5 .

 45 Figure 4. Stability regions of synchronized equilibria as a function of k and τav, for T = 0 (dashed curves) and T = 0.01, n = 1 (full curves).

1 Figure 6 .

 16 Figure 6. Simulations for 4 ring coupled systems. Zoom of Figure 5 (bottom) for the time-windows [58, 60] (left) and [98 100] (right). Notice the phase shifts of π/2.

Figure 7 .Figure 8 .Figure 9 .

 789 Figure 7. Simulation for 4 ring coupled systems, corresponding to parameter values indicated with a cross in Figure 3 (left). The delay values are τ = 0.15, n = 1, T = 0.01.

  )(y l (t -θ) -y i (t -θ))dθ , i = 1, . . . , p,

The work of I.-C. Morȃrescu and M. Jungers has been supported by the ANR project COMPACS-"Computation Aware Control Systems", ANR-13-BS03-004-02. The work of W. Michiels has been supported by the Programme of Interuniversity Attraction Poles of the Belgian Federal Science Policy Office (IAP P6-DYSCO), by OPTEC, the Optimization in Engineering Center of the KU Leuven, and the project G.0712.11N and G.0717.11N of the Research Foundation -Flanders (FWO).