Fernando Peruani

Monojit Choudhury

Animesh Mukherjee

Niloy Ganguly

Emergence of a non-scaling degree distribution in bipartite networks: a numerical and analytical study

Keywords: 89, 75, -k -Complex systems 89, 75, Fb -Structures and organization in complex systems

We study the growth of bipartite networks in which the number of nodes in one of the partitions is kept fixed while the other partition is allowed to grow. We study random and preferential attachment as well as combination of both. We derive the exact analytical expression for the degree-distribution of all these different types of attachments while assuming that edges are incorporated sequentially, i.e., a single edge is added to the growing network in a time step. We also provide an approximate expression for the case when more than one edge are added in a time step. We show that depending on the relative weight between random and preferential attachment, the degree-distribution of this type of network falls into one of four possible regimes which range from a binomial distribution for pure random attachment to an u-shaped distribution for dominant preferential attachment.

A bipartite network is a graph which connects two distinct sets (or partitions) of nodes, which we will refer to as the top and the bottom set. An edge in the network runs between a pair of a top and a bottom node but never between a pair of top or a pair of bottom nodes (see Fig. 1). Typical examples of this type of networks include collaboration networks such as the movie-actor [1][2][3][4][5][START_REF] Amaral | Proc. Natl. Acad. Sci[END_REF], articleauthor [START_REF] Newman | [END_REF][8][9][10][11], and board-director [12,13] network. In the movie-actor network, for instance, the movies and actors are the elements of the top and the bottom set respectively, and an edge between an actor a and a movie m indicates that a has acted in m. The actors a and a ′ are collaborators if both have participated in the same movie, i.e., if both are connected to the same node m ′ . The concept of collaboration can be extended to include so diverse phenomena represented by bipartite networks as the citypeople network [14], in which an edge between a person and a city indicates that the person has visited that particular city, the word-sentence [15,16], bank-company [17] or donor-acceptor network, which accounts for injection and merging of magnetic field lines [18].

Several models have been proposed to synthesize the structure of bipartite networks when both partitions grow unboundedly over time [1][2][3][4]16]. It has been found that for such growth models when each incoming top node con-nects through preferential attachment to bottom nodes the emergent degree distribution of bottom nodes follows a power-law [1]. Another important property of bipartite networks is that the clustering coefficient cannot be measured in the standard way [2], and has to be measured as a cycle of four connections [19].

On the other hand, bipartite networks, where one of the partitions remains fixed over time (i.e., the number of bottom nodes are constant), have received comparatively much less attention. Recently it was shown through numerical simulations that restrictions in the growth rate of the partitions can lead to non-scaling degree distribution highly sensitive to the parameters of the growth model [20]. However, there is still no systematic and analytical study of this kind of networks. Realizations of this type of bipartite networks include numerous relevant systems such as the interaction between the codons and genes as well as amino acids and proteins in biology and elements and compounds in chemistry. We can also include in this group those networks in which one partition can be considered to be in a pseudo-steady state while the other one keeps on growing at a much faster rate. For instance, it is reasonable to assume that for the city-people network [14], the city growth rate is zero compared with the population growth rate. Other examples of this type of networks could be the phoneme-language [START_REF] Choudhury | Proceedings of COLING-ACL[END_REF][START_REF] Mukherjee | [END_REF], in linguistics, or train-station [START_REF] Sen | [END_REF]24] in logistics.

In this work, we study the growth of bipartite networks in which the number of nodes in one of the partitions is kept fixed. We explore random and preferential attachment as well as a combination of both. We obtain an exact analytical expression for the degree distribution of the bottom nodes assuming that the attachment is sequential, i.e., a single edge is added to the network in one time step. We also present an approximate solution for the case of parallel attachment, i.e., when more than one edge are incorporated into the network in a given time step. We show that, depending on the relative weight of random to preferential attachment, the degree-distribution of this type of network falls into one of four possible regimes which range from a binomial distribution for pure random attachment to an u-shaped distribution for dominant preferential attachment. For combinations of random and preferential attachment the degree-distribution asymptotically tends to beta-distribution with time.

The growth model. -We consider the case in which the top partition grows with time while the number of nodes in the bottom partition N is kept constant. We grow the network in the following way. At each time step a new node is incorporated to the top set. Then, µ edges are connected from the new node to the nodes in the bottom set (see Fig. 1). The probability of attaching a new edge to the bottom node i is A(k t i ), where k t i refers to the degree of the bottom node i at time t. We refer to A(k t i ) as the attachment kernel and define it as:

A(k t i ) = γk t i + 1 N j=1 (γk t j + 1) (1)
where the sum in the denominator runs over all bottom nodes, and γ is a model parameter which controls the relative weight of random to preferential attachment. γ can be thought of as γ = 1/α, where α is a positive constant known in previous models as initial attractiveness [START_REF] Dorogovtsev | Evolution of Networks: From Biological Nets to the Internet and WWW[END_REF].

There is a subtlety related to the attachment kernel and µ that is worth to mention. The stochastic process can be performed in such a way that the attachment of the µ incoming nodes is done sequentially, i.e., one edge is attached per time step. This implies that the denominator of A(k t i ) has to be updated for each incoming node (and hence an edge), and that t refers to the event of incorporating a new edge to the bottom set. Alternatively, the attachment of the µ new edges can be done in parallel. This implies that the new µ edges have all the same probability of attaching to bottom node i. In this case, t refers to the event of incorporating a new node to the top set.

There are two significant limits to consider: γ = 0 and γ -→ ∞. For γ = 0, Eq. ( 1) reduces to A(k t i ) = 1/N , which implies that all bottom nodes have the same probability of being selected by an incoming edge. This limit corresponds to pure random attachment. For γ -→ ∞ Eq. ( 1) reduces to A(

k t i ) = k t i / N j=1 (k t j )
, which means that higher degree bottom nodes have higher probability of being selected. This case corresponds to pure preferential attachment.

Stochastic simulations have been performed with the initial condition where all bottom nodes at time t = 0 have zero degree, i.e., initially no edges are connected to the bottom nodes.

Evolution equation for sequential attachment. -Now we aim to derive an evolution equation for the degree distribution of the bottom nodes. We focus on sequential attachment. Let p k,t be the probability of finding a randomly chosen bottom node with degree k at time t. We recall that t refers to the t-edge attachment event. p k,t can be defined as p k,t = n k,t /N , where n k,t refers to the number of nodes in the bottom set with degree k at time t, and ... denotes ensemble average, i.e. average over realizations of the stochastic process described above. We express the evolution of p k,t in the following way:

p k,t+1 = (1 -A(k, t))p k,t + A(k -1, t)p k-1,t (2)
where A(k, t) refers to the probability that the incoming edge lands on a bottom node of degree k. A(k, t) can be easily derived from Eq. ( 1) and takes the form:

A(k, t) = γk + 1 γt + N (3)
The reasoning behind Eq. ( 2) is the following. The probability of finding a bottom node with degree k at time t + 1 decreases due to the number of nodes, which have a degree k at time t and receive an edge at time t + 1 therefore acquiring degree k + 1, i.e., A(k, t)p k,t . Similarly, this probability increases due to the number of nodes that at time t have degree k -1 and receives an edge at time t + 1 to have a degree k, i.e., A(k -1, t)p k-1,t . Hence the net increase in the probability can be expressed as in Eq. ( 2).

According to what was done in the stochastic simulations, we assume that at time t = 0 all bottom nodes have zero degree, which implies the initial condition p k,t=0 = δ k,0 , where δ k,0 is the Kronecker delta function.

Exact solution for sequential attachment. -We look for the exact analytical solution of sequential attachment. For this purpose we express Eq. (2) as:

p t+1 = M t p t = t τ =0 M τ p 0 ( 4
)
where p t denotes the degree distribution at time t and is defined as p t = [p 0,t p 1,t p 2,t . . .] T , p 0 is the initial condition expressed as p 0 = [1 0 0 . . .] T , and M τ is the evolution matrix at time τ which is defined as:

M τ =      1 -A(0, τ ) 0 0 0 . . . A(0, τ ) 1 -A(1, τ ) 0 0 . . . 0 A(1, τ ) 1 -A(2, τ ) 0 . . . . . . . . . . . . . . . . . .     
(5) Since our initial condition is a vector with zeros in all position except in the first one, all the relevant information, i.e., the degree distribution of the bottom nodes, is in the first column of t τ =0 M τ . A close inspection to the evolution of this column, explicitly using Eq. ( 3), reveals that the k-th element of it, which corresponds to p k,t , can be expressed as:

p k,t = t k k-1 i=0 (γi + 1) t-1-k j=0 (N -1 + γj) t-1 m=0 (γm + N ) (6)
for k ≤ t, and p k,t = 0 for k > t, and where we have defined -1 i=0 (...) = 1, and

t k refers to the combina- torial number t!/[(t -k)!k!].
Eq. ( 6) is the exact solution of Eq. ( 2) using as initial condition p k,t=0 = δ k,0 , i.e., the analytical expression of the degree distribution of the bottom nodes when sequential attachment is applied.

In the limit of γ = 0 Eq. ( 6) reduces to:

p k,t = t k 1 N k 1 - 1 N t-k (7)
for k ≤ t, and p k,t = 0 for k > t. In other words, Eq. ( 7) is the solution of the sequential problem when pure random attachment is applied.

Parallel attachment. -We focus on parallel attachment, i.e., when more than one edge are added per time step. We do not aim to derive an exact analytical expression for the degree distribution of this problem but provide a reasonable approximation. We recall that for parallel attachment t refers to the event of incorporating a new top node. We assume that µ ≪ N and expect Eq. ( 2) to be a good approximation of the process after replacing A(k, t) with A p (k, t). We define A p (k, t) as The term µ in the denominator appears since in this case the total degree of the bottom nodes at any point in time is µt rather than t as in Eq. ( 3). The numerator contains a µ since at each time step there are µ edges that are being incorporated into the network rather than a single edge.

A p (k, t) = (γk + 1) µ γµt + N . ( 8
)
It is important to mention here that Eq. ( 2) cannot exactly represent the stochastic parallel attachment because it explicitly assumes that in one time step a node of degree k can only get converted to a node of degree k +1. Clearly, the incorporation of µ edges in parallel allows the possibility for a node of degree k to get converted to a node of degree k + µ. The correct expression for the evolution of p k,t reads:

p k,t+1 = (1 - µ i=1 A(k, i, t))p k,t + µ i=1 A(k -i, i, t)p k-i,t (9)
where A(k, i, t) represents the probability at time t of a node of degree k of receiving i new edges in the next time step. We expect Eq. ( 2) to be a good approximation of Eq. ( 9) when A(k, 1, t) ≫ A(k, i, t) where i > 1.

The solution of Eq. ( 2) with the attachment kernel given by Eq. ( 8) reads:

p k,t = t k k-1 i=0 (γi + 1) t-1-k j=0 N µ -1 + γj t-1 m=0 γm + N µ ( 10
)
We expect Eq. ( 10) to approximate the degree distribu-tion of the stochastic process with parallel attachment for µ ≪ N . This means that we cannot expect the approximation to hold for large values of γ or µ/N . In the limit of random attachment, i.e., γ = 0, Eq. ( 10)

becomes p k,t = t k µ N k 1 -µ N t-k .
Figs. 2(a)-(c) and 3(a)-(b) show a comparison between stochastic simulations and Eq.( 10) and prove that Eq. ( 10) is a good approximation for µ ≪ N and relative low values of γ. For large values of γ, as said above, the approximation fails. However, for µ = 1 Eq. ( 10) reduces to Eq. ( 6), which in this case is the exact solution, and then the theory works for all values of γ (see Figs. 2(d) and 3(c)).

From random to preferential attachment. -Fig. 2 shows that there is a clear transition from random to preferential attachment. At γ = 0 (see Fig. 2(a)) we observe that p k,t is centered around the maximum (mode of the distribution) which shifts with time at a speed of µ/N per time step, while the width of the distribution also spreads with time. This behavior corresponds to a situation in which all bottom nodes receive roughly the same amount of edges with time. The well defined maximum tells us about the average number of edges each bottom node has, while the variance of the distribution indicates the presence of fluctuations around that mean value which increases with time.

For 0 < γ < 1 the distribution is not longer symmetric (see Fig. 2(b)). Bottom nodes having small degree rarely receive an edge, and so p k,t decays slowly for small value of k. However the distribution still exhibits a maximum, mode of the distribution, which shifts with time (see Fig. 3(a)).

For 1 ≤ γ ≤ (N/µ)-1 the distribution looses the (local) maximum and becomes monotonically decreasing (see Fig. 2(c)). We can always find a bottom node with small degree because small degree nodes hardly get an edge. On the other hand, there are very few nodes with high degree, and these ones receive almost all incoming edges. The temporal evolution of the distribution for this range of γ is shown in Fig. 3(b).

For γ > (N/µ)-1 the distribution described by Eq. ( 10) exhibits an u-shape. As said above, we cannot expect Eq. (10) to approximate the stochastic process for such large values of γ. Stochastic simulations performed in this range of γ for µ > 1 are very noisy and the u-shape cannot be obtained by averaging over few simulations. However, for µ = 1 we can illustrate the u-shaped distribution in a clear way, see Fig. 2(d). As in the previous case, we still can find a bottom node with small degree because small degree nodes hardly get an edge (see maximum at k = 0). But on the other hand, we can be sure that there is at least one node with very large degree, because the node with the largest degree at time t -1 very likely is going to get an edge at time t, in an effect like "winner takes all" (see peak at k = t). The node with largest degree keeps on increasing its degree with time, and so p k,t has two peaks, one located at k = 0 and the other one at k = t (see Fig.

3(c)).

Beta Distribution. -In the following we offer a quantitative analysis of the transition by showing that p k,t behaves asymptotically with time as a beta-distribution for γ > 0.

For t ≫ η, where η = N/(γµ), we can approximate the products in Eq. ( 10) by gamma-functions and apply Stirling's approximation. After some algebra we obtain:

p k,t ≃ C -1 (k/t) γ -1 -1 (1 -k/t) η-γ -1 -1 (11)
where C is a normalization constant defined by 3 shows a comparison between stochastic simulations (circles), the theoretical solution given by Eq. ( 10) (black dashed curve), and the approximation given by Eq. ( 11) (red solid curve) for two different values of γ at various times. In Fig. 3 it can be observed that p k,t approaches asymptotically to Eq. ( 11) (compare the black dashed and the red solid curves). Notice that Eq. ( 11) does not have any fitting parameter and represents a betadistribution f (x; α, β) of the variable x = k/t and fixed parameters α = γ -1 and β = η -γ -1 .

C = t 0 (k ′ /t) γ -1 -1 (1 -k ′ /t) η-γ -1 -1 dk ′ . Fig.
For 0 < γ < 1, α > 1 and β > 1 the mode of the distribution is given by (α -1)/(α + β -2) = ((γ -1 ) -1)/(η -2). This can be easily verified by taking the first derivative of Eq. ( 11) equal to zero. From this we learn that in this range of γ the maximum of the distribution k max is located at k max = t((γ -1 ) -1)/(η -2) (see Fig. 3(a)). In the limit of γ -→ 0 we retrieve the behavior of k max observed for pure random attachment, i.e., k max = t(µ/N ).

At γ = 1, α = 1 and β > 0, the moving peak is no longer found, i.e., the mode of the distribution is located for all times at k max = 0 (see Fig. 3(b)). This condition also holds for 1 < γ ≤ (N/µ) -1.

For γ > (N/µ)-1, there is another regime for the degree distribution. For α < 1 and β < 1, p k,t becomes u-shaped with a peak fixed at k = 0 and the other one shifting with t. For µ = 1 the additional peak is located at k = t (see Fig. 3(c)).

Concluding remarks. -We have studied the growth of bipartite networks in which the number of nodes in the bottom set is kept fixed. We consider random and preferential node attachment as well as the combination of both. We have derived the degree distribution evolution equation for sequential and parallel attachment of nodes. For sequential attachment we have provided the exact analytical solution of the problem. For parallel attachment we have obtained an approximate expression for the degree distribution. Through simulations we have provided numerical evidence which shows that the approximation for parallel attachment is reasonable when µ ≪ N and γ is small. Finally, we have shown that for both, sequential and parallel attachment, the degree-distribution falls into one of four possible regimes: a) γ = 0, a binomial distribution whose mode shifts with time, b) 0 < γ < 1, a skewed distribution which exhibits a mode that shifts with time, c) 1 ≤ γ ≤ (N/µ) -1, a monotonically decreasing distribution with the mode frozen at k = 0, and d) γ > (N/µ) -1, an u-shaped distribution with peaks at k = 0 and k = t.

Our results might be useful to explain the dynamical growth of various systems like the speech sound inventories of the world's languages, which can be rendered a bipartite structure as explained through the phoneme-language network in [START_REF] Choudhury | Proceedings of COLING-ACL[END_REF][START_REF] Mukherjee | [END_REF]. A detailed study of the parameter γ leading to the degree distribution of the network can then shed some light on the amount of randomness/preference that has gone into shaping the evolution of such a system.

Fig. 1 :

 1 Fig. 1: Scheme of a growing bipartite network. In the example, the number of fixed bottom nodes is N = 4. Each top nodes arrives with µ = 2 edges. The top node D represents a new incoming node.

Fig. 2 :

 2 Fig. 2: The four possible degree distributions depending on γ. Symbols represent average over 5000, in (a)-(c), and 50000, in (d), stochastic simulations. The dashed curve is the theory given by Eq. (10). From (a) to (c), t0 = 1000, N = 1000 and µ = 20. (a) at γ = 0, p(k, t) becomes a binomial distribution. (b) γ = 0.5, the distribution exhibits a maximum which shifts with time for 0 ≤ γ < 1. (c) γ = 1, p(k, t) does not longer exhibit a shifting maximum and the distribution is a monotonically decreasing function of k for 1 ≤ γ ≤ (N/µ) -1. (d) γ = 2500, t0 = 100, N = 1000 and µ = 1. p(k, t) becomes an u-shaped curve for γ > (N/µ) -1.

Fig. 3 :

 3 Fig. 3: Temporal evolution of p(k, t) for various values of γ. Symbols represent average over 5000, in (a)-(b), and 50000, in (c), stochastic simulations. The black dashed curve is the theory, represented through Eq. (10). The red solid curves correspond to the approximation given by the beta-distribution, Eq. (11). (a) γ = 0.5, and (b) γ = 1. N = 1000 and µ = 20. Compare with Fig. 2(b) and (c). (c) γ = 2500, N = 1000 and µ = 1, see Fig. 2(d).

This work was supported by the Indo-German (DST-BMBT) project grant. The authors would like to extend their gratitude to Dr. Andreas Deutsch, Dr. Lutz Brusch and Dr. Luis G. Morelli for their valuable comments and suggestions. MC and AM would like to thank Media Lab Asia and Microsoft Research India respectively for financial assistance. They would also like to extend their gratitude to Prof. Anupam Basu for providing them with the laboratory infrastructure. FP would also like to acknowledge the hospitality of IIT-Kharagpur and financial support through Grant No. 11111.


Fernando Peruani

Tobias Klauss

Andreas Deutsch

Anja Voss-Boehme

Traffic jams, gliders, and bands in the quest for collective motion

Keywords: numbers: 87.18.Gh, 87.10.Hk, 05.65.+b, 05.70.Ln

We study a simple swarming model on a two-dimensional lattice where the self-propelled particles exhibit a tendency to align ferromagnetically. Volume exclusion effects are present: particles can only hop to a neighboring node if the node is empty. Here we show that such effects lead to a surprisingly rich variety of self-organized spatial patterns. As particles exhibit an increasingly higher tendency to align to neighbors, they first self-segregate into disordered particle aggregates. Aggregates turn into traffic jams. Traffic jams evolve toward gliders, triangular high density regions that migrate in a well-defined direction. Maximum order is achieved by the formation of elongated high density regionsbands -that transverse the entire system. Numerical evidence suggests that below the percolation density the phase transition associated to orientational order is of first-order, while at full occupancy it is of second-order.

Self-propelled particle (SPP) systems are found at all scales in nature. Examples in biology range from human crowds [1] and animal groups [START_REF] Cavagna | Proc. Natl. Acad. Sci[END_REF][START_REF] Bhattacharya | [END_REF], down to insects [4,5], bacteria [START_REF] Zhang | Proc. Natl. Acad. Sci[END_REF], and even to the microcellular scale with, e.g., the collective motion of microtubules driven by molecular motors [START_REF] Schaller | [END_REF]. SPP systems are not restricted to living systems. There are examples in non-living matter, as for instance in driven granular media [8][9][10]. Interestingly, the statistical properties of the large-scale self-organized patterns emerging in SPP systems depend only on a few microscopic details: the symmetry associated to the selfpropulsion mechanism of the particles, which can be either polar [11][12][13] or apolar [14], the symmetry of the velocity alignment mechanism, which can be either ferromagnetic [11,12] or nematic [13,14], and very importantly, the presence or absence of volume exclusion effects, as we are going to discuss here. In addition, the nature of the supporting space where the particles move plays also a crucial role: the dimension of the space, and whether this space is continuous [11][12][13][14] or discrete [15][16][17][18].

In this letter, we focus on polar SPPs moving on a two-dimensional lattice that align their orientation, respectively their moving direction, via a local ferromagnetic alignment mechanism. We explicitly model volume exclusion effects: nodes can be occupied at most by one particle. We show that such effects introduce a coupling between particle speed, local density and local alignment that lead to a surprisingly rich variety of selforganized spatial patterns unseen in previous swarming models. As particles exhibit an increasingly higher tendency to align to neighbors, the system passes through three distinct phases. For weak alignment strength, the system exhibits orientational disorder, while particles self-segregate, Fig. 1(a). Within this initial phase, there is a transition from a spatially homogeneous to an ag-gregate phase. The onset of orientational (polar) order marks the beginning of the second phase which is characterized by the emergence of locally ordered, high density regions: traffic jams, Fig. 1(b). As the tendency to align is enhanced, traffic jams evolve toward triangular high density aggregates that migrate in a well-defined direction. We refer to these dynamical traffic jams as gliders, Fig. 1(c). The third phase emerges when the particles self-organize into highly ordered, elongated, high density regions: bands, Fig. 1(d). In contrast to the traveling bands observed in off-lattice SPP models with ferromagnetic alignment [12], these bands are formed by particles aligned to the long axis of the band and are rather static.

We find evidence that the phase transition to orientational order is discontinuous below the percolation threshold. When the lattice is fully occupied, the system reduces to the classical planar Potts model and the phase transition to orientational order is undoubtedly of second-order [24,25]. Previous lattice swarming models were found to exhibit a continuous phase transition from a homogeneous to a condensed phase in 1D [17,18], while in 2D, both, first and second-order transitions to orientational order have been claimed. Bussemaker et al. reported a second-order transition in a cellular automaton model with 4 moving directions [15], while Csahók and Vicsek found, for a lattice-gas model with 6 moving directions, a weakly-first-order transition to collective migration [16]. Here we show that all these phenomena occur in a minimal 2D lattice swarming model, but where in contrast to previous models the system dynamics is dominated by volume exclusion effects that lead to a completely novel spatial self-organization of particles.

Model. Particles move on a two-dimensional lattice with periodic boundary conditions, and have four possible orientations: up, down, left, and right, whose associated vectors are v 1 , v 2 , v 3 , and v 4 , respectively. The state of a particle is given by its position on the lattice and its orientation. As it will become clear below, the orientation of a particle fully determines its moving direction. The particles are able to perform two actions: i) they can change their orientation, and ii) they can migrate in the direction given by their orientation. These actions have associated transition rates which specify the average number of events per time unit. Let us start out with the reorientation transition rate T R that a particle at x with orientation v turns its orientation into direction w:

T R ((x, v) → (x, w)) = exp(g y ∈A(x) w|V(y) ) , (1)
where the sum runs over the nearest lattice neighbors of x, represented by A(x). The vector V(y) returns the orientation of the particle placed in node y, if there is any, and the null vector otherwise. The symbol .|. indicates the inner product between two vectors, while g is a parameter which controls the alignment sensitivity. For positive g, eq. ( 1) defines a stochastic ferromagnetic alignment mechanism. As result of this alignment, first nearest neighbor particles tend to be aligned.

The migration rate is defined in the following way:

T M ((x, v) → (y = x + v, v)) = v 0 if node y is empty 0 if node y is occupied. (2)
Thus, a particle at position x and pointing in direction v migrates to the neighboring node y = x + v with a transition rate v 0 as long as the node at y is empty. If the y-node is occupied, the particle will not jump. The only action that is allowed to the particle in this situation is to change its orientation. At this point, it is important to understand how this continuous-time process is simulated. Let us assume that at time t 0 the system is in a given state. Then, we compute the time at which the next event will take place in the system, i.e., we calculate t 1 . Now we have to decide which is the event that will take place at time t 1 . We choose at random one out of all the possible events that could take place, but we weight each of these events according to their associated transition rate, eqs. (1)

and (2)

. This procedure is an adaption of the classical Gillespie algorithm [19] to interacting particle systems [20].

Results. We start by fixing the migration rate v 0 and particle density d, and use g as control parameter. The degree of orientational order in the system at time t is characterized by the (global) orientation

m(t) = (1/N )| x V(x)|,
where N is the total number of particles in the system, the sum runs over all lattice sites, and V(x) is defined as above. Fig. 2 shows the behavior of the mean orientation m as function of the alignment sensitivity g, with . . . a temporal average taken once the system reaches the steady-state after a short transient. The system exhibits a phase transition to orientational order above a critical g * for all densities, as long as v 0 > 0. Here we focus on 0 < d < d p , where d p refers to the (site) percolation threshold in a 2D (square) lattice, d p ∼ 0.59. The system exhibits three phases with g, labeled I, II, and III by increasing alignment sensitivity g. Phase I corresponds to g < g * and is characterized by exhibiting no macroscopic orientational order. Fig. 3 shows that within phase I there is a dynamic phase transition from a spatially homogeneous to an aggregate phase as g is increased. The degree of aggregation is characterized by the average cluster size k , which is computed as the temporal average of k = k k p(k, t), with p(k, t) = k n k (t)/N , where n k (t) is the number of clusters of mass k at time t and N = d L 2 is the number of particles in the system. The figure shows that there exists a critical value g a above which a phase transition to aggregation occurs, with g a < g * . As g approaches g * , more than 85% of the particles in the system form a large aggregate. The transition point between phase I and II is at g * , where the curves m (g, L) for different system sizes L meet. Interestingly, g * seems to be independent of the density d, as confirmed with simulations for various system sizes with density d = 0.2, 0.3, and 0.4 (data not shown). Moreover, g * coincides with the critical point for the full occupancy problem, i.e., d = 1. obtained from the m(t) time series for values of g close to g * do not exhibit a Gaussian shape as expected for a second order transition, but rather a bimodal distribution (the arrow indicates the second peak) as expected for first order transitions. The coexistence of several particle configurations (or "phases") exhibiting different degrees of ordering, i.e., values of m, is evident for values of g deep into phase II. Fig. 4(c) shows that m(t) jumps between well-defined values corresponding to different spatial patterns. Values of m(t) ∼ 1/ √ 2 are associated to gliders, while higher values of m(t) correspond to bands. Lower values of m(t) are due to the presence of traffic jams. Gliders are dynamical traffic jams moving backwards with respect to their mean average orientation. Their presence affects the temporal evolution of the center of mass of the system, x cm (t), which exhibits ballistic motion whenever there is a glider, while otherwise is Brownian. As result of this, the average speed of the center of mass V cm peaks at values of g where gliders are more stable [START_REF]Vcm has been computed as Vcm = |xcm(t + τ )xcm(t)|/τ t, with τ = 10 4 and[END_REF], see Fig. 5. Gliders are remarkably different from traffic jams observed in 2D traffic models [START_REF] Biham | [END_REF], arguably due to the presence of the alignment mechanism, Eq. ( 1). How frequently gliders appear and for how long they survive, depends on the value of g and L. For example, for g = 1.4, m(t) displays excursions from low to high values that reflect the fact the system alternates between traffic jams and gliders. For an illustration of this dynamics, see [23]. For g > g b , i.e. in phase III, the only stable configurations is a band. In contrast to g * , g b is highly dependent on L.

Discussion. In the limit of full occupancy, d = 1, particles are frozen in their positions and the only action allowed to them is reorientation. The system defined by Eqs. ( 1) and ( 2) becomes an equilibrium system in this limit, whose order is again characterized by m(t). Since, by definition, there are only four possible orientations, we can safely claim that the model reduces to a 4-state planar Potts model [24]. It has been shown that this model can be reduced to the standard 2-Potts model [25].

In two-dimensions, the standard q-Potts model exhibits a continuous transition for q ≤ 4, and a discontinuous one for q > 4 [24]. Thus, our model exhibits a secondorder transition for d = 1, as confirmed via simulations (data not shown). For v 0 > 0 and d < 1, we are in a pure non-equilibrium scenario. The migration rule, Eq.

(2), breaks detailed balance and prevents us from writing down a free energy. Nevertheless, it is worth to compare our system with its equilibrium counterpart, the diluted Potts model with annealed vacancies. If we represent the absence of particles in a lattice position with an extra vector direction, we end up with a 5-Potts system with full occupancy, instead of a diluted 4-Potts system. In consequence, according to what was said above, the transition would be first-order. The argument applies to the standard Potts model and assumes vacancies are in thermal equilibrium, which is not true for Eq. ( 2). Nevertheless, it helps to realize that a discontinuous dynamic phase transition is quite possible in our non-equilibrium system.

In summary, we have shown through a minimal model that volume exclusion effects, when they are allowed to stop particle motion, can lead to a surprisingly rich variety of self-organized patterns. Such effects introduce a coupling between local density, local orientation and particle speed that strongly affects the large-scale behavior of the system, with the jamming of particles playing a dominant role. This coupling is present in many real systems as in gliding bacteria, animal groups, etc. Certainly, several features of the self-organized patterns described here depend on the discrete nature of the model. Nevertheless, we expect similar phenomena to emerge in off-lattice, continuum symmetry systems. For instance, static traffic jams are probably a robust property of all systems where stagnation can occur. Here we have also learned that the jamming of self-propelled particles can lead to unexpected self-organized structures in two-dimensions like dynamical traffic jams, e.g., gliders. The presence of an alignment mechanism induces (local) orientational order, and provided particles are oriented, density waves of stagnated particles should emerge. The results reported here are a first step toward a deep understanding of the possible phenomena that such a coupling may induce.

We thank A. Greven and C.F. Lee for useful comments.

FIG. 1 :

 1 FIG. 1: (Color online) Example of self-organized spatial patterns: (a) disordered aggregates, (b) traffic jams, (c) gliders, and (d) bands. Particle orientation is indicated by the orientation of the small triangles, and is also color-coded: right (red), left (black), down (green), and up (blue). Parameter values are indicated in Fig. 2.

FIG. 2 :

 2 FIG. 2: (Color online) Mean orientation m vs. the alignment sensitivity g for systems with density d = 0.3 and migration rate v0 = 100. Simulations were carried out for 2 10 7 time steps. The different curves correspond to different system sizes L. The boundary between phase I and II, g * , and the between phase II and III, g b are indicated by the vertical dashed lines. (a) to (d) refer to the simulation snapshots shown in Fig. 1.

FIG. 3 :

 3 FIG.3:(Color online) Phase transition to aggregation (phase I). (a) to (d) correspond to simulation snapshots whose value of g is indicated in (f) (other parameters as in Fig.2). The average cluster size k normalized by the total number of particles, d.L 2 , as function of g is shown in (e) for various system sizes, and in (f) for various densities.

Fig. 4 (

 4 FIG. 4: (Color online) Phase transition to orientational order. Time series m(t), histogram, and typical spatial configurations for phase II, g * < g ≤ g b . Letters in (c) indicate simulation snapshots shown in (d)-(f). Parameters as in Fig. 2.

FIG. 5 :

 5 FIG. 5: (Color online) (a) average speed of the center of mass Vcm as function of g. (b) Trajectories of the center of mass for two different values of g.


R Fónod

D Gontkovič
email: daniel.gontkovic@tuke.sk

Fault Estimation in a Class of First Order Nonlinear Systems

Keywords: Sliding mode state estimator, radial basis function neural network, Lyapunov function, steepest descent rule, linear matrix inequalities

Reformulated principle of fault estimation design for one class of first order continuous-time nonlinear system is treated in this paper, where a neural network is regarded as model-free fault approximator. The problem addressed is presented as approach based on sliding mode methodology with combination of radial basis function neural network to design robust nonlinear fault estimation. The method utilizes Lyapunov function and the steepest descent rule to guarantee the convergence of the estimation error asymptotically. Simulation results show the feasibility of the proposed approach.

I. INTRODUCTION

The theory of sliding mode has emerged as a method capable of use in given robust control systems, as well as state estimation [START_REF] Utkin | Variable structure systems with sliding modes[END_REF], [START_REF] Zak | Systems and Control[END_REF]. The main advantage of sliding mode in state estimation exists in the sliding where the state error trajectories are constrained on the predetermined equivalent system. The sliding mode observer design can be then divided into two phases. First, in the reaching phase, give an error switching law such that the estimated error trajectory is trapped on a switching surface and remain on it thereafter. And second, in the sliding phase, determine the switching surface such that the error dynamics in the sliding mode have good performance.

It should be pointed out that the robustness of a variable structure resides in its sliding phase, but not in its reaching phase. There is no easy way to shape the error dynamics of the reaching phase. In fact, the error dynamic is not completely robust over all time. In addition, if the system suffers from uncertainties or faults, then the error behavior in the sliding mode is not only governed by the switching surface but also determined by faults. In this case, the state estimator stability may not be assured [START_REF] Chen | Design and implementation sliding mode controller based on radial basis function neural network for synchronous reluctance motor[END_REF], [START_REF] Gao | Discrete-time variable structure control systems[END_REF]. However, in many cases sliding mode control is supposed [START_REF] Veselý | Robust model predictive control design with input constraints[END_REF], [START_REF] Yu | Stability analysis of visual servoing with sliding-mode estimation and neural compensation[END_REF], [START_REF] Young | A control engineer's guide to sliding mode control[END_REF], and different observer structures used to estimate sensor, actuator, and system faults [START_REF] Chee | Sliding mode observers for detection and reconstruction of sensor faults[END_REF], [START_REF] Choi | LMI-based sliding-mode sbserver design method[END_REF], [START_REF] Filasová | On sensor faults estimation using sliding mode observers[END_REF].

Recently, neural networks has been successfully used in model-free fault estimation in nonlinear systems. The known supports of these methods give their ability to good performance for approximation of nonlinear functions. Combining sliding mode observer technology with neural networks to design fault estimation structures offer new opportunity in design [START_REF] Chen | Sensor fault diagnosis for a class of time-delay uncertain nonlinear systems using neural network[END_REF].

Radial basis function (RBF) neural networks seems to be one of the neural networks with high approximation and regularization capability [START_REF] Buhmann | Radial Basis Functions: Theory and implementation[END_REF], where the essential phenomenological rationale for the use of RBF rest in the realm of the purpose of feed-forward networks and feature extraction possibility. Gaussian RBF are employed most frequently, since it is bounded, strictly positive and continuous on IR, and the optimization of RBF networks provides better approximation and interpolation capability as compared to the sigmoid functions. However, the performance of RBF neural networks depends on the number of neurons in hidden units, and on methods used for determining the output neuron weights, where training of a RBF neural network is, in general, a challenging nonlinear optimization problem.

The paper is focused on robust system fault estimation and isolation (FDI) approach for one class of nonlinear single input/output (SISO) systems. The FDI scheme is based on sliding mode neural observer [START_REF] Filasová | Some properties of sliding mode residual filter design[END_REF], which is robust against system uncertainty, and fault estimation be realized using the sliding boundary size. When a system fault is occurred the estimate part in the observer for faults is enabled. When a RBF neural network is used to approximate a system fault the fault estimation can be formulated in a simpler way.

The main contribution of the paper is to present a reformulated design method for continuous-time nonlinear SISO systems to solve system fault estimation problem using a sliding mode observer. Used structure is motivated by the need for robustness in model-free fault estimation. Under defined conditions, the stability of the observer is assured, where the state observer error is bounded. In comparison with methods based on analytical models (see e.g. [START_REF] Chen | Robust Model-Based Fault Diagnosis for Dynamic Systems[END_REF], [START_REF] Ding | Model-Based Fault Diagnosis Techniques: Design Schemes, Alghorithms, and Tools[END_REF], [START_REF] Krokavec | Dynamic Systems Diagnosis[END_REF], [START_REF] Nobrega | LMI-based filter design for fault detection and isolation[END_REF]) design conditions for specified system class are here derived using radial basis function neural network.

II. PROBLEM DESCRIPTION

The system under consideration is a first order SISO nonlinear dynamic systems, which is given as q(t)=aq(t)+bu(t)+f (q(t))

(1)

y(t)=q(t) (2)
where q(t) ∈ IR, u(t) ∈ IR,a n dy(t) ∈ IR are the state, input and output variables, respectively, and scalars a ∈ IR, b ∈ IR. A system fault is modelled by an unknown additive function f (q(t)) ∈ IR.

The problem of interest is to design sliding-mode observer based estimation of system faults.

III. OBSERVER STRUCTURE

Since the calculation of f (q(t)) may be not accurate enough, in order to rectify this kind of errors as well as to reduce the number of sensors a sliding mode observer can be used to estimate f (q(t)), i.e. the sliding mode observer of system ( 1), ( 2) can be designed e.g. as qe (t)=aq e (t)+bu(t)+j(y(t)-y e (t)) + s(t) (3)

y e (t)=q e (t) (4)
where q e (t) ∈ IR is the state observer variable, y e (t) ∈ IR is the estimated output, and j is an observer gain. The term s(t) ∈ IR, which represents the unknown f (q(t)), is defined as a switching function of the tracking error of the observer, i.e.

s(t)=̺ sign(e(t)) (5)
where e(t) is the observer estimation error defined as

e(t)=q(t) -q e (t) (6)
and ̺ is a constant factor. The constant ̺ has to be chosen in such way that sliding mode is enforced in the manifold ν e = q(t)-q e (t)=0. Once sliding mode is enforced, the differential equations ( 1) and ( 3) have the same solution and the terms f (q(t)) and s(t) have to be equivalent. During application a low pass filter has to be used to gain the average value s(t) of the discontinuous time function [START_REF] Chen | Robust Model-Based Fault Diagnosis for Dynamic Systems[END_REF]. Model-free sliding mode observer does not require the model of f (q(t)), but it cannot arrive finite time convergence, the sliding mode gain should be bigger than an upper bound, and incomplete information about the rest nominal parameters causes chattering phenomenon. It seems to be natural to construct its estimate f (q(t)|w(t)) depending on a parameter w(t), which can be adjusted online by means of an updating law. In the next, a neural network is used to approximate f (q(t)) and construct a neural observer of the form qe (t)=aq e (t)+bu(t)+ +j(y(t)-y e (t)) + f a (q e (t)|w a (t)) + (1-r(t)) s(t)

)

y e (t)=q e (t) (8)
f a (q e (t)|w a (t)) = w T a (t)ϕ(q e (t)) (9)
Here s(t) ∈ IR is an external feed-forward compensation signal and r(t) is a switch function. The switch function switch between the neural estimator and the sliding mode observer in the dependence on the output error (6), i.e.

r(t)= 1 if e 2 (t) ≥ δ 0 if e 2 (t) <δ ( 10
)
where δ is known upper bound of the neural modelling error.

If r(t)=1 , the observer is pure neural observer, and ( 7), (8) becomes qe (t)=aq e (t)+bu(t)+ +j(y(t)-y e (t)) + f a (q e (t)|w a (t))

)

y e (t)=q e (t) (12)
If after time t 0 e 2 (t) <δ,i . e .r(t)=0then w T a (t) is a constant vector w T a = w T a (t 0 ), and the observer ( 7), (8) becomes pure sliding mode observer qe (t)=aq e (t)+bu(t)+ +j(y(t)-y e (t)) + f a (q e (t)|w a )+s(t)

)

y e (t)=q e (t) (14)
IV. RBF NEURAL NETWORK

The model-free fault estimation employed in the next use a radial basis function neural network (RBFNN) to approximate an unknown fault combining with the sliding mode observer.

Supposing that there are p receptive field units in the neural network hidden layer than the output of the RBFNN is

f a (q e (t)) = p h=1 w h (t)ϕ h (q e (t)) = w T a (t)ϕ(q e (t)) (15)
ϕ h (q e (t)) = exp ⎛ ⎝ - (q e (t) -c h ) 2 σ 2 h ⎞ ⎠ ( 16
)
w T a (t)= w 1 (t) w 2 (t) ••• w p (t) (17)
ϕ T (q e (t)) = ϕ 1 (q e (t)) ϕ 2 (q e (t)) ••• ϕ p (q e (t)) (18) where c h , σ h are the center and width (spread factor) of the neural cell of the h-th neuron in the hidden layer, respectively, and w h (t) ∈ IR is the weight connecting the h-th hidden layer neuron and the network output.

Thus, the optimal weight values of RBF neural network can be considered as follows

w ⋆ =arg min wa∈Ω f sup qe∈Sq |f (q e (t)|w a )-f (q e (t))| (19)
where

Ω f = {w a : w a ≤δ w } (20)
is a valid field of w a (t), δ w is a parameter, and S q ∈ IR is a variable space of the state observer variable. Because radial basis function ϕ(q e (t)) satisfies Lipschitz condition ϕ(q(t))ϕ(q e (t)) 2 ≤ κ ϕ (q(t) -q e (t)) 2 (21)

where κ ϕ can be selected by users, then according to the Stone-Weierstrass theorem, the smooth function f (q(t)) can be written as

f (q(t)) = w ⋆T ϕ(q(t)) + ε(t) (22)
where w ⋆T is a fixed weight matrix of the neural network, and ε(t) is is the smallest approximation (modelling) error satisfying condition

ε 2 (t) ≤ δ (23)
Then, according to (15) is f (q(t)) -f a (q e (t)|w a (t)) = ν(t)= = w ⋆T ϕ(q(t)) + ε(t)w T a (t)ϕ(q e (t)) = = w ⋆T (ϕ(q(t))-ϕ(q e (t))) + ε(t)+ +(w ⋆T -w T a (t))ϕ(q e (t))

where with notations w ⋆T ∆ϕ(t)=w ⋆T (ϕ(q(t))-ϕ(q e (t))) (25)

w T e (t)=w ⋆T -w T a (t) (26)
it yields ν(t)=ε(t)+w ⋆T ∆ϕ(t)+w T e (t)ϕ(q e (t))

V. S LIDING MODE NEURAL OBSERVER STABILITY

The neural sliding mode observer design method given by [START_REF] Filasová | Some properties of sliding mode residual filter design[END_REF] is represented as follows. Assuming that the sliding mode neural observer switches between two models then if r(t)=1then (1), ( 11), [START_REF] Chen | Sensor fault diagnosis for a class of time-delay uncertain nonlinear systems using neural network[END_REF], and (27) implies

ė(t)=aq(t)+bu(t) -aq e (t) -bu(t)- -j(y(t)-y e (t)) + f (q(t))-f a (q e (t)|w a (t)) = =(a -j)e(t)+ν(t) (28)
ė(t)= =(a-j)e(t)+ε(t)+w * T ∆ϕ(t)+w T e (t)ϕ(q e (t))

(29)
respectively, where ∆ϕ(t)=ϕ(q e (t)+e(t))ϕ(q e (t))

Since neural networks have to be discrete updated and e(i) be used in the updating law, then

w e (i +1)=w e (i) -r(i)μ(i)ϕ(q e (i))e(i)= = w e (i) -μ(i)ϕ(q e (i))e(i) (31)
where 0 <μ(i) < 1.

Defining Lyapunov function candidate as

v(w e (i)) = w e (i) 2 = trace w e (i)w T e (i) (32)
then

∆v(w e (i)) = w e (i +1) 2 -w e (i) 2 =
= w e (i) -μ(i)ϕ(q e (i))e(i) 2 w e (i)

2 = = μ 2 (i) ϕ(q e (i)) 2 e(i) 2 - -2μ(i) w e (i) ϕ(q e (i)) e(i) (33)
respectively. Inserting (29) it yields [START_REF] Krokavec | Optimal Stochastic Systems[END_REF] 2μ(i) w e (i) ϕ(q e (i)

) ≥ ≥ 2μ(i) w T e (i)ϕ(q e (i)) = =2μ(i) ė(i)-(a-j)e(i)-ε(i)-w * T ∆ϕ(i) ≥ ≥ μ(i)2 e(i)w * T ∆ϕ(i) - -2e(i)(a-j)e(i) -2 e(i)ė(i) -2 e(i)ε(i) (34)
and it is obvious that (21) implies

2 e(i)w * T ∆ϕ(i) ≤ e(i) 2 + w * T ∆ϕ(i) 2 ≤ ≤ (1 + κ ϕ w * T )e 2 (i) (35)
Since for a smooth stable system and any i>0 yields

ė(i) e(i) ≤ i (36) then -2 e(i)ė(i) ≥-i e(i) 2 = -ie 2 (i) (37)
and -2e(i)(a-j)e(i)=-2(a-j)e 2 (i) ≥ 0 (38)

Because it yields, too

-2 e(i)ε(i) ≥-e(i) 2 -ε(i) 2 = -e(i) 2 -δ (39) then 2μ(i) w T e (i)ϕ(q e (i)) ≥-μ(i)δ+ +μ(i)(1+κ ϕ w * T -i -2(a-j)-1)e 2 (i) (40)
Subsequently, (33) implies

∆v(w e (i)) ≤ μ(i)δ- -μ(i) -μ(i) ϕ(q e (i)) 2 + +κ ϕ w * T -2(a-j) -i e 2 (i) < 0 (41)
Setting

μ(i)= κ ϕ w * T -2(a-j) -i 1+ ϕ(q e (i)) 2 (42)
where μ i > 0 is the learning rate, (33) takes the form ∆v(w e (i)) ≤-μ(i)e 2 (i)+μ(i)δ<0

Since in this case e 2 (i) ≥ δ,t h e n∆v(w e (i)) < 0, v(w e (i)) is bounded, and also w e (i) is bounded.

Here the initial condition w * T can be selected as

κ ϕ w * T < 1+i +2(a -j) (44)
which gives

0 <κ ϕ w * T -i -2(a -j) < 1 (45)
and it is evident that (42) implies 0 <μ(i) < 1.

If at time t = t 0 be r(t)=0the weights become constants, the observer ( 7), ( 8) becomes pure sliding mode taking form [START_REF] Krokavec | Dynamic Systems Diagnosis[END_REF], [START_REF] Nobrega | LMI-based filter design for fault detection and isolation[END_REF]. The error dynamics obtained from (1), ( 2), [START_REF] Krokavec | Dynamic Systems Diagnosis[END_REF], and ( 14) is

ė(t)=aq(t)+bu(t) -aq e (t) -bu(t) -s(t)- -j(y(t)-y e (t)) + f (q(t))-f a (q e (t)|w a )= =(a -j)e(t)+∆f (t) -s(t) (46)
e y (t)=e(t)

where ∆f (t) is neural modelling error when the weight of the neural networks is fixed as w a .

Defining the Lyapunov function of the form v(e(t)) = pe(t)ė(t)=e(t)p(aq(t)+bu(t)+∆f (t)--s(t)-aq e (t)-bu(t)-j(q(t)-q e (t)) (49) v(e(t)) = e(t)p(a e e(t) -s(t)+∆f (t)) < 0

v(e(t)) = 1 2 pe 2 (t) (48)
where a e = a -j

Supposing ∆f (t)=s(t) and setting pj =1(50) implies

p(a -j)=pa -1 < 0 (52)
Thus, for any p satisfying (52) be j = p -1

be introduced conditioned by e(t) =0 ,a n ds(e(t)) = 0 if e(t)=0. This function characterizes the sliding mode in state estimation, and it is evident that with a positive p satisfying (52) the observer error is bounded. The sliding mode neural observer requires two design parameters: switch constant δ and the upper bound ∆ of neural modelling error ∆f (t) when start the sliding mode compensation. Here δ decide when neural network learning is stopped and sliding mode observer is started. The bigger δ is, the shorter training time the neural observer has and since the neural modelling error is bigger, so ∆ should be bigger. Usually ∆ >δ, because δ corresponds to the modelling error with optimal synaptic weight, while ∆ corresponds to the modelling error when e 2 (t) <δ. Note, this parameters are user-defined.

VI. ILLUSTRATIVE EXAMPLE

Considering the synchronous reluctance motor drive system [START_REF] Chen | Design and implementation sliding mode controller based on radial basis function neural network for synchronous reluctance motor[END_REF] described by the equation q(t)=aq(t)+bu(t)+f (q(t))

where q(t)=ω(t) is velocity, the input control u(t) is the electromagnetic torque, J =0 .00076 kgm 2 is the

u(t)=-kq(t)+nω r (t)+ 1 T I t 0 (ω r (τ ) -ω(τ ))dτ
where k =0.0075, n =0.0076, T I =1.5sec.

In this simulation was tested different numbers of hidden nodes with the nest result that after the hidden nodes number is more than 5 the estimation accuracy be not improved a lot.

There is consider no parameter variations except system fault at t =5sec. The simulation results depicted in Figs. 1 and 2 represent the system velocity response and its estimation, and show that the proposed control approach is effective. Figs. 3 and4 provide the fault estimation and synaptic weight adaptation obtained using the sliding mode neural observer.

VII. CONCLUDING REMARKS

This paper explores use of a sliding mode scheme for fault estimation in nonlinear system of the first order. Modified design method is presented where system fault estimation is achieved by the RBF neural network inclusion into the sliding mode observer. The stability analysis is given using the Lyapunov method and LMI observer condition to garanty observer error convergence, where LMI procedure give desirable observer dynamics. It is demonstrated that the sliding mode observer can be used to detection and estimate system faults. The ability to estimate the faults directly is very desirable for fault detection, and preferred to handle problem of fault occurrence and faults behavior.

Fig. 1 .

 1 Fig. 1. Responses of the controlled system

Fig. 2 .

 2 Fig. 2. The corresponding observer output response

Fig. 3 .Fig. 4 .

 34 Fig. 3. System fault and its estimation

ACKNOWLEDGEMENTS

The work presented in this paper was supported by VEGA, Grant Agency of Ministry of Education and Academy of Science of Slovak Republic under Grant No. 1/0256/11. This support is very gratefully acknowledged.


Róbert Fónod
email: robert.fonod@student.tuke.sk

Dušan Krokavec
email: dusan.krokavec@tuke.sk

Actuator Fault Estimation Using Neuro-Sliding Mode Observers

Reformulated principle for designing actuator fault estimation for continuous-time linear MIMO systems, based on neuro-sliding mode observer structure, is presented in this paper. Radial basis function neural network is used as a model-free fault approximator of the unknown additive fault. The method utilizes Lyapunov function and the steepest descent rule to guarantee the convergence of the estimation error asymptotically, where the design parameters can be obtained using LMI techniques. Finally, the proposed fault estimation scheme is applied to a nonlinear water tank system and simulation results illustrate its satisfactory performance.

I. INTRODUCTION

A particular class of observer-based methods uses sliding mode observers (SMO). In SMO a nonlinear discontinuous term is injected into the observer dynamics depending on the output estimation error. These observers are more robust than Luenberger observers, as the discontinuous term enables the observer to reject disturbances, and also a class of mismatch between the system and observer. Supported by existence of linear matrix inequality (LMI) solvers (see e.g. [START_REF] Peaucelle | User's Guide for SeDuMi Interface 1.04[END_REF]), the intention is to reformulate given problem and optimize it over LMI constraints [START_REF] Choi | LMI-based sliding-mode observer design method[END_REF], [START_REF] Filasová | On sensor faults estimation using sliding mode observers[END_REF].

Neural networks can be considered as an alternative modelfree observer because it offers much potential benefit for nonlinear modeling. The dynamic neural networks have been also applied to design a Luenberger type of observer [START_REF] Kim | Nonlinear observer design using dynamic recurrent neural networks[END_REF]. Combining sliding mode observer techniques with neural networks to design fault estimation structures offer new opportunity in design [START_REF] Chen | Sensor fault diagnosis for a class of time delay uncertain nonlinear systems using neural network[END_REF]. The main idea of this combination is to apply them at the same time, where sliding mode term is used to compensate the neural modeling error and the system uncertainties, while the neural network is used to approximate the unknown (in general nonlinear) fault function.

Radial basis function (RBF) neural networks seem to be one of the neural networks with high approximation and regularization capability [START_REF] Buhmann | Radial basis functions: theory and implementations[END_REF], where the essential phenomenological rationale for the use of RBF rest in the realm of the purpose of feed-forward networks and feature extraction possibility. Gaussian RBF are employed the most frequently, since it is bounded, strictly positive and continuous on IR, and the optimization of RBF networks provides better approximation and interpolation capability as compared to the sigmoid functions. However, the performance of RBF neural networks (RBFNN) depends on the number of neurons in hidden layer, and on methods used for determining the network weights, where training of a RBFNN is, in general, a challenging nonlinear optimization problem.

Actuator faults detection and diagnosis (FDD) algorithm for linear systems is presented in this paper. The FDD scheme is based on neuro-sliding mode observer [START_REF] Resendiz | Discrete-time sliding mode neural observer for continuous time mechanical systems[END_REF], which realizes a fault estimation using RBF neural network. When a fault occurs, the neural part in the observer is enabled. The switching function avoids the unnecessary chattering caused by sliding mode phenomena.

The main contribution of the paper is to adapt the work presented in [START_REF] Fónod | Fault estimation in a class of first order nonlinear systems[END_REF], [START_REF] Krokavec | Estimation of nonlinear system parameter faults via sliding mode observer[END_REF] to actuator faults, acting on the system in additive form. Moreover, the idea is extended to a class of continuous-time linear Multi-Input Multi-Output (MIMO) systems, using a bank of observers. Under defined conditions, the stability of the observer is assured, where the state observer error is bounded. In comparison with methods based on analytical models (see e.g. [START_REF] Choi | LMI-based sliding-mode observer design method[END_REF], [START_REF] Krokavec | Discrete-Time Systems[END_REF], [START_REF] Ding | Model-based fault diagnosis techniques: design schemes, algorithms, and tools[END_REF]) design conditions for specified system class are here derived using learning approach and sliding mode technique.

II. PROBLEM DESCRIPTION

Consider a state-space representation of the linear MIMO system in canonical form [START_REF] Filasová | On sensor faults estimation using sliding mode observers[END_REF] with only actuator faults

q(t) = Aq(t) + B u(t) + f a (t) = = Aq(t) + Bu(t) + r i=1 b i f i a (q(t)) (1)
y(t) = Cq(t) (2)
where q(t) ∈ IR n , u(t) ∈ IR r , and y(t) ∈ IR m are the state, input and output variables, respectively, and A ∈ IR n×n , B ∈ IR n×r , C ∈ IR m×n are real matrix parameters of the system. Actuator faults are modelled by an unknown additive vector function f a (q(t)) ∈ IR r , b i ∈ IR n is the i th column vector of the input matrix B, and f i a (q(t)) ∈ IR represents the i th actuator fault.

The task is to design a neuro-sliding mode observer (NSMO) based estimation of the actuator faults using RBF neural networks.

III. NEURO-SLIDING MODE OBSERVER STRUCTURE

The output variable y(t) and the output matrix C can be expressed as

y(t) =      y 1 (t) y 2 (t)
. . .

y m (t)      =      c T 1 c T 2 . . . c T m      q(t) (3)
Thus, the system given by ( 1), (2) can be reduced, for the fault estimation purposes, in the following form

q(t) = Aq(t) + Bu(t) + r i=1 b i f i a (q(t)) (4)
y j (t) = c T j q(t) (5)
where the y j (t) is the j th output variable of the system, and c T j has to satisfy rank c j A T c j . . .

(A T ) n-1 c j T = n (6)
i.e. the pair (A, c T j ) has to be an observable pair (see e.g. [START_REF] Krokavec | Discrete-Time Systems[END_REF]).

In the next, we consider single actuator faults, i.e. only one actuator fault is occurring in the system. The particular NSMO of the simplified system (4), [START_REF] Chen | Sensor fault diagnosis for a class of time delay uncertain nonlinear systems using neural network[END_REF] for the i th actuator fault f i a (q(t)) is expressed in the form qe (t) = Aq e (t) + Bu(t) + l(y j (t)y e (t))

+ +b i f i a (q e (t)|w a (t)) + (1 -r(t))hν(t) (7)
y e (t) = c T j q e (t) (8)
where q e (t) ∈ IR n is the state observer variable, y e (t) ∈ IR is the estimated output, l ∈ IR n is the traditional Luenberger observer gain vector to be designed, ν(t) ∈ IR is an external feed-forward compensation signal, f i a (q e (t)|w a (t)) ∈ IR is the unknown actuator fault approximation using RBFNN, and r(t) is a switch function.

If rank(B) = r, then for each actuator fault f i a (q(t)), it is necessary to construct an own observer of the structure ( 7), [START_REF] Fónod | Fault estimation in a class of first order nonlinear systems[END_REF], i.e. if all actuator faults are monitored, then a bank of r observers (see Fig. 1) have to be designed. In the next, for the sake of simplicity, the index i will be omitted, and the following equivalences are hold

f a (q(t)) ⇔ f i a (q(t)) (9) fa (q(t)|w a (t)) ⇔ f i a (q e (t)|w a (t)) (10)
g ⇔ b i (11)
The state e q (t) and the output e y (t) estimation error are defined as follows e q (t) = q(t)-q e (t), e y (t) = y j (t)-y e (t) = c T j e q (t) [START_REF] Shevitz | Lyapunov stability theory of nonsmooth systems[END_REF] The switch function, which switches between the pure neural and sliding mode observer in the dependence on the output error, i.e. with a = m T m

r(t) = 1, if e y (t) 2 a ≥ τ ; faulty 0, if e y (t) 2 a < τ ; healthy ( 13
)
where τ is a pre-specified threshold of the weighted modelling error, and m will be defined formally later.

IV. NEURAL OBSERVER STRUCTURE AND STABILITY

If r(t) = 1, then ( 7), [START_REF] Fónod | Fault estimation in a class of first order nonlinear systems[END_REF] becomes

qe (t) = Aq e (t) + Bu(t) + l(y j (t) -y e (t))+ +g fa (q e (t)|w a (t)) (14
)

y e (t) = c T j q e (t) (15)
and the observer is pure neural observer.

A. RBF Neural Networks

The model-free fault estimation employed in the next use a radial basis function neural network (RBFNN) to approximate the unknown actuator fault f a (q(t)). Supposing that there are p receptive field units in the neural network hidden layer then the output of the RBFNN is fa (q e (t)|w a (t)) = p h=1 w h (t)ϕ h (q e (t)) = w T a (t)ϕ(q e (t))

(16)

ϕ h (q e (t)) = exp   - q e (t) -c ⋄ h 2 σ 2 h   (17) w T a (t) = w 1 (t) w 2 (t) • • • w p (t) (18)
ϕ T (q e (t)) = ϕ 1 (q e (t)) ϕ 2 (q e (t))

• • • ϕ p (q e (t)) (19)
where c ⋄ h ∈ IR n , 0 < σ h ∈ IR are the center and width (spread factor) of the neural cell of the h th neuron in the hidden layer, respectively, and w h (t) ∈ IR is the weight connecting the h th hidden layer neuron and the network output.

Thus, the optimal weight values of RBF neural network can be expressed as follows

w ⋆ = arg min wa∈Ω f sup q e ∈Sq |f (q e (t)|w a )-f (q e (t))| (20)
where

Ω f = {w a : w a ≤ δ w } (21)
is a valid field of w a (t), δ w is a upper bound parameter, and S q ⊆ IR n is a variable space of q e (t) ∈ S q .

Radial basis function ϕ(q e (t)) satisfies the so-called Lipschitz condition ϕ(q(t))ϕ(q e (t)) 2 ≤ κ ϕ q(t)q e (t) 2 (22

)
where • denotes Euclidian norm of a vector, and κ ϕ > 0 is the Lipschitz constant of (17).

Defining the RBFNN approximation error as

ψ(t) = f a (q(t)) -fa (q e (t)|w a (t)) (23)
then, according to the Stone-Weierstrass theorem, the smooth function f a (q(t)) can be written as

f a (q(t)) = w ⋆T ϕ(q(t)) + ε(t) (24)
where w ⋆T is the fixed weight vector (20), and ε(t) is is the smallest approximation (modelling) error satisfying condition

ε 2 (t) ≤ δ (25)
Then, according to ( 16), ( 23) is

ψ(t) = w ⋆T ϕ(q(t)) + ε(t) -w T a (t)ϕ(q e (t)) = = w ⋆T (ϕ(q(t))-ϕ(q e (t))) + ε(t)+ +(w ⋆T -w T a (t))ϕ(q e (t)) (26)
with notations

∆ϕ(t) = ϕ(q(t)) -ϕ(q e (t)) (27)
w T e (t) = w ⋆T -w T a (t) (28)
it yields

ψ(t) = w ⋆T ∆ϕ(t) + w T e (t)ϕ(q e (t)) + ε(t) (29)

B. Neural Observer Error Dynamics

Combining (4), ( 12), ( 14) and (29) implies the neural observer state error estimation dynamic ėq (t) = Aq(t) + Bu(t) -Aq e (t) -Bu(t)--l(y j (t)-y e (t))+g f a (q(t))-fa (q e (t)|w a (t)) = = (Alc T j )e q (t) + gψ(t)

(30)

ėq (t) = A e +g w * T ∆ϕ(t)+w T e (t)ϕ(q e (t))+ε(t) (31)
respectively, where A e = Alc T j . Since only e y (t) is available, then (12) can be modified as

c j e y (t) = (c j c T j + λI)e q (t) -λIe q (t) (32)
where 0 < λ ∈ IR, and e q (t) = (c j c T j +λI) -1 c j e y (t)+λ(c j c T j +λI) -1 e q (t) (33) e q (t) = me y (t) + λN e q (t) (34)

me y (t) = (I -λN )e q (t) (35)
respectively, where

N = (c j c T j +λI) -1 , m = N c j , Υ = I -λN (36)

C. Neural Network Updating Law

Since NN have to be discrete updated, and only e y (t) can be used in the updating law, the neural network weights can be determined by using the steep descent method w e (i + 1) = w e (i)r(i)µ(i)ϕ(q e (t))g T me y (t)

where µ(i) > 0 is the neural network learning rate.

D. Neural Observer Stability

Defining Lyapunov function candidate as

v(w e (i)) = w T e (i)w e (i) = w e (i) 2 (38)
then with r(i) = 1 yields ∆v(w e (i)) = w e (i + 1)

2 -w e (i) 2 = = w e (i) -µ(i)ϕ(q e (i))g T me y (i) 2 -w e (i) 2 < 0 (39) ∆v(w e (i)) = µ 2 (i) ϕ(q e (i))g T 2 me y (i) 2 - -2µ(i) w T e (i)ϕ(q e (i))g T me y (i) ≤ ≤ α(i) -µ(i)ξβ(i) < 0 (40)
respectively, where

α(i) = µ 2 (i) gϕ T (q e (i)) 2 me y (i) 2 (41)
β(i) = 2 e T q (i)gw T e (i)ϕ(q e (i)) , ξ = Υ (42)
w T e (i)ϕ(q e (i))g T = gw T e (i)ϕ(q e (i))

Since from the Lyapunov theory of stability [START_REF] Shevitz | Lyapunov stability theory of nonsmooth systems[END_REF] yields e T q (i) ėq (i) < 0, then (43) with (31) gives β(i) = = 2 e T q (i) ėq (i)-A e e q (i)-g(ε(i)+w ⋆T ∆ϕ(i)) ≤ ≤ 2 e T q (i)gw ⋆T ∆ϕ(i) -2 e T q (i)A e e q (i) --2 e T q (i)gε(i) -2 e T q (i) ėq (i) (44) For a smooth stable system and any i > 0 yields 2 ėq (i) ≤ i e q (i) (45) then -2 e T q (i) ėq (i) ≥ -i e q (i)

(46)

Since A e < 0 is designed to be stable, then second term in (44) can be rewritten

-2 e T q (i)A e e q (i) ≤ -σ max (A e ) e q (i)

2 ( 47
)
where σ max (A e ) is the largest singular value of the matrix A e . Third term in (44) can be rewritten in the form

-2 e T q (i)gε(i) ≥ -e q (i) 2 -g ε(i) 2 = = -e q (i) 2 -g δ (48)
It is obvious that (22) implies 2 e T q (i)gw ⋆T ∆ϕ(i) ≤ e q (i)

2 + gw ⋆T ∆ϕ(i)) 2 ≤ ≤ 1 + κ ϕ gw ⋆T e q (i) 2 (49)
then with notation σ max (A e ) = σ m µ(i)ξβ(i) ≤ -µ(i)ξ g δ+ +µ(i)ξ κ ϕ gw ⋆Tσ mi e q (i) 2 (50) Thus, combining (41) and (50) gives

∆v(w e (i)) = α(i) -µ(i)ξβ(i) = µ(i)ξ g δ- -µ(i) -µ(i) gϕ T (q e (i)) 2 +1-1 me y (i) 2 + + κ ϕ gw ⋆T -σ m -i ξ e q (i) 2 (51)
Setting

µ(i) = κ ϕ gw ⋆T -σ m -i 1 + gϕ T (q e (i)) 2 (52)
then the Lyapunov function difference (51) takes the form

∆v(w e (i)) = µ(i)ξ g δ -µ 2 (i) me y (i) 2 < 0 (53) ∆v(w e (i)) = -µ(i) µ(i) e y (i) 2 a -ξ g δ < 0 (54)
Because e y (i) 2 a ≥ τ , it is enough to select, such a constant δ, that the sign inside the brackets of ( 54) is positive, then the neural network is stable trained.

Here the initial condition w ⋆T can be selected as

κ ϕ gw ⋆T < 1 + i + σ m (55)
which gives

0 < κ ϕ gw ⋆T -i -σ m < 1 ( 56
)
and it is evident that (52) implies 0 < µ(i) < 1.

V. SLIDING MODE OBSERVER STRUCTURE AND STABILITY

If at time t = t 0 is r(t) = 0, then the weights become constant, i.e. w a = w a (t 0 ), and the observer ( 7), ( 8) becomes pure sliding mode taking the form qe (t) = Aq e (t) + Bu(t) + l(y j (t)y e (t))+ +g fa (q e (t)|w a ) + hν(t) (57)

y e (t) = c T j q e (t) (58)
The discontinuous term ν(t) enters the observer dynamics through the design vector h ∈ IR n , and is defined by

ν(t) = -ρ(t)sign(e y (t)) = -ρ(t) ey(t) ey(t)
if e y = 0, 0 otherwise (59)

where ρ(t) > 0, ρ(t) ∈ IR is the scalar gain function, high enough to enforce the sliding motion.

A. Sliding Mode Observer Stability Conditions

The error dynamics obtained from ( 4), ( 12) and ( 57) is

ėq (t) = A -lc T j e q (t) + g∆f a (t) -hν(t) (60)
where ∆f a (t) = f a (q(t))-fa (q e (t)|w a ) (61)

is the neural modelling error when the weights of the neural network are fixed as w a .

Defining the Lyapunov function of the form v e q (t) = e T q (t)P e q (t) > 0 (62)

where P = P T > 0, P ∈ IR n× , and taking the time derivative of v e q (t) results in v e q (t) = ėT q (t)P e q (t) + e T q (t)P ėq (t) < 0 (63) inserting ( 60) into (63) gives

v e q (t) = e •T q (t)P • e • q (t) < 0 ( 64
)
where e •T q (t) = e T q (t) ∆f a (t) ν(t) (65)

P • =   P (A-lc T j )+(A-lc T j ) T P P g -P h * 0 0 * * 0   < 0 (66)
(Hereafter, * denotes the symmetric item in a symmetric matrix). Defining the congruence transform matrix

T q =   I I I 0   (67)
then pre-multiplying left-hand side of (66) by T q , and righthand side of (66) by T T q gives P • = T q P • T T q = = P (A-lc T j )+(A-lc T j ) T P P (gh) * 0

Since (68) can be rewritten as

P A+A T P P (g -h) * 0 - - P 0 l c T j 0 - c j 0 l T P 0 < 0 (69)
then using the orthogonal complement matrix

c •⊥ j = c j 0 ⊥ = c ⊥ j I (70)
and pre-multiplying left side of (69) by c •⊥ j and right side of (69) by c •⊥T

j gives c ⊥ j P A+A T P c ⊥T j c ⊥ j P (g -h) * 0 < 0 (71) Setting h = g (72)
which gives sufficient condition for the stable sliding mode existence.

VI. DESIGN CONDITIONS

The observer gain vector l can not be directly computed from (73), hence introducing a new LMI variable z ∈ IR n and using assignment

z = P l (74) then (69) implies
P A + A T P -zc T j -c j z T < 0 (75)
The observer gain vector l can be set as

l = P -1 z (76)
Solving the above specified inequalities (73), (75) with respect to the LMI variables P , z the sliding mode observer (57), ( 58) is asymptotically stable. The observer parameters are then given by ( 72) and (76). The sliding mode neural observer [START_REF] Resendiz | Discrete-time sliding mode neural observer for continuous time mechanical systems[END_REF], [START_REF] Fónod | Fault estimation in a class of first order nonlinear systems[END_REF] requires two additional design parameters: switch constant τ and the upper bound ∆ of neural modelling error ∆f a (t) when start the sliding mode compensation. Here, τ decide when the neural network learning is stopped and sliding mode observer is started. The bigger τ , the shorter training time the neural observer has and since the neural modelling error is bigger, so ∆ should be bigger. Usually ∆ > τ , because τ corresponds to the modelling error with optimal synaptic weights, while ∆ corresponds to the modelling error when e y (i) 2 a < τ .

VII. ILLUSTRATIVE EXAMPLE

Considering a nonlinear model of three coupled water tank system [START_REF] Youssef | Piecewise sliding mode decoupling fault tolerant control system[END_REF], linearized around the working points and transformed into the canonical form results in 1 for the first, and c T 2 for the second observer, and solving the inequalities (73) with respect to the LMI variables P 1 and P 2 , the problem was solved as feasible with Self-Dual-Minimization (SeDuMi) package for Matlab [START_REF] Peaucelle | User's Guide for SeDuMi Interface 1.04[END_REF], and giving the solutions

l 1 = P -1 1 z 1 =   0.1141 0.2229 0.0014   , l 2 = P -1 2 z 2 =   0.0096 0.0007 0.1522   It is evident that (72) implies h 1 = b 1 =   0.000 1.234 0.000   , h 2 = b 2 =   0.000 0.000 1.234  
The result of several simulations was, that after the hidden neurons number is more than 7, the estimation accuracy is not improved a lot. The chosen parameters for RBF neural network are µ 1 = µ 2 = 0.2, τ 1 = 10 -3 , τ 2 = 5.10 -4 . The centers c ⋄ h and widths σ h were determined using the Matlab function newrb, and were left out due to the space limitations. The sliding mode gain was chosen as ρ 1 = ρ 2 = 1.

The presented simulation results show, that the method based on neuro-sliding mode observer structure is very well adapted to obtain a faithful estimation of a actuator fault signal, see Fig. 2. The residual signal, Fig. 3, indicates only the fault existence, but says nothing about the shape and amplitude of 4 shows, that the the second RBFNN output is not coupled by the occurrence of the first actuator fault. The results depicted in Fig. 5 and 6 represents the first system output and its estimation, and show the synaptic weights adaptation.

VIII. CONCLUDING REMARKS

In this paper, an approach for actuator fault detection and estimation for a class of linear MIMO systems is investigated. Reformulated design method is presented where the fault estimation is achieved by inclusion of the RBF neural network into the sliding mode observer. The neural networks weights are updated based on a steepest descent rule. Unlike many previous methods in the literature, the proposed fault estimation scheme does not rely on the availability of full state measurements. The stability analysis is given using the Lyapunov method and LMI observer condition to guaranty observer error convergence. Under certain conditions, the estimated signal can approximate the fault signal to any accuracy even in the presence of unmodelled dynamics. Simulation results provided demonstrate and illustrate the effectiveness and capabilities of our proposed actuator fault estimation strategy.

ACKNOWLEDGEMENT

The work presented in this paper was supported by VEGA, Grant Agency of Ministry of Education and Academy of Science of Slovak Republic under Grant No. 1/0256/11. This support is very gratefully acknowledged.

Fig. 1 .

 1 Fig. 1. Bank of observers for all actuator faults estimation

Fig. 2 .

 2 Fig. 2. First actuator fault and its estimation

Fig. 3 .A

 3 Fig. 3. First output estimation error (residual)

2 z 1

 21 second inequality (75) with respect to LMI vector variables z 1 and z = 0.0077 0.1474 0.0014 T z 2 = 0.0160 -0.0011 0.1522 T Thus, l 1 , l 2 can be set according to (76) as

Fig. 4 .Fig. 5 .

 45 Fig. 4. Second actuator fault and its estimation

7 Fig. 6 .Fig. 7 .

 767 Fig. 6. Synaptic weights adaptation of the first observer