Chris Muelder
email: muelder@cs.ucdavis.edu

Arnaud Sallaberry
email: arnaud.sallaberry@lirmm.fr

Kwan-Liu Ma
email: ma@cs.ucdavis.edu

Improved Cluster Tracking for Visualization of Large Dynamic Graphs

Analysis and visualization of dynamic graphs is a challenging problem. Clustering can be applied to dynamic graphs in order to generate interactive visualizations with both high stability and good layout quality. However, the existing implementation is naïve and unoptimized. Here we present new algorithms to improve both the temporal clustering results and the efficiency of the cluster tracking calculation, and evaluate the results and performance.

Introduction

In recent years, the domain of network visualization has yielded many techniques for exploring large graphs. Most of these deal with static networks and are based on graph drawing algorithms (see for example [START_REF] Hachul | An experimental comparison of fast algorithms for drawing general large graphs[END_REF] or [START_REF] Muelder | Rapid graph layout using space filling curves[END_REF]), clustering techniques (see [START_REF] Schaeffer | Graph clustering[END_REF] for an introduction), or exploratory methods (see for example [START_REF] Van Ham | Interactive visualization of small world graphs[END_REF], [START_REF] Abello | ASK-GraphView: A large scale graph visualization system[END_REF] or [START_REF] Archambault | GrouseFlocks: Steerable exploration of graph hierarchy space[END_REF]). In contrast, fewer works have been devoted to the exploration of dynamic graphs. A dynamic graph is an evolving graph where vertices and edges are added and removed over time. Examples of such graphs include social networks, dependency graphs in software engineering, website hyperlinks, router networks, collaboration networks, etc.

When creating a node-link diagram for a dynamic graph, not only does the layout need to consider graph topology, but also the stability between time-steps. This generally forces a trade-off between layout quality and stability, as a perfectly stable layout would sacrifice layout quality, and naively calculating ideal layouts would not offer stability. While there are a number of existing methods for creating these layouts, they have not been shown to scale well to large dynamic graphs.

A dynamic clustering based approach for visualization of dynamic graphs can provide an overview of the entire dynamic graph over time, yield high quality layouts for every timestep, minimize node motion between time steps to provide stability and preserve the user's mental map, and allow for interactive exploration even under random access patterns. One existing approach consists of first clustering each time-step independently to guarantee good locality for every time-step, then tracking the clusters between time-steps, and finally arranging the clusters and nodes such that nodes that are consistant are stationary and transitional node motion is minimized [START_REF] Sallaberry | Clustering, visualizing, and navigating for large dynamic graphs[END_REF]. This produces a temporal arrangement which we directly visualize as a timeline, and which is used to define layouts for a node link diagram for each time-step which both meets general layout criteria (namely cluster co-location and short average edge length) and where node motion is minimized between time-steps.

However, there are some limitations to this approach. First, by only considering pairwise timesteps, it is impossible to track clusters that disappear for some amount of time before reforming, which results in the formation of extraneous dynamic clusters. Second, the association between timesteps was previously cacluated very naïvely, with pairwise comparison and Jaccard calculation for each possible association.

Here, we describe improvements to this approach that resolve both of these issues. The first issue is addressed by an algorithm for preserving unused dynamic clusters from old time steps and allowing them to be reincorporated into new time steps. And the second issue is addressed through a heavy optimization of the association algorithm that not only reduces the number of clusters that have to be compared but which also enables direct calculation of the Jaccard index. The result of the combination of these improvements is a more compact, robust, dynamic clustering that can be computed much more efficiently.

Related Works

A common method for visualizing dynamic graphs is to animate the transitions between time-steps [START_REF] North | Incremental layout in DynaDAG[END_REF][START_REF] Diehl | Graphs, they are changing[END_REF][START_REF] Erten | GraphAEL: Graph animations with evolving layouts[END_REF][START_REF] Garey | Dynamic graph drawing of sequences of orthogonal and hierarchical graphs[END_REF][START_REF] Boitmanis | Visualizing internet evolution on the autonomous systems level[END_REF][START_REF] Frishman | Online dynamic graph drawing[END_REF]. This approach yields dynamic visualization with nodes appearing, disappearing and moving to produce readable layout for each time-step. Alternatively, multiple time-steps can be staticly placed next to each other using "Small Multiples" [START_REF] Tufte | Envisionning Information[END_REF]. This eases the comparison of distant time-steps but the area devoted for each time-step is small and this reduces the readability of each graph. An empirical study to compare the advantages and drawbacks of these approaches ("Animation" vs. "Small Multiples") has been performed by [START_REF] Archambault | Animation, small multiples, and the effect of mental map preservation in dynamic graphs[END_REF]. A major issue for both methods is to ensure the stability of the layout [START_REF] Kumar | Visual exploration of complex time-varying graphs[END_REF][START_REF] Frishman | Online dynamic graph drawing[END_REF][START_REF] Brandes | A quantitative comparison of stress-minimization approaches for offline dynamic graph drawing[END_REF][START_REF] Hu | Embedding, clustering and coloring for dynamic maps[END_REF]. A stable layout helps preserve the user's mental map as there is less movement between time-steps, but sacrifices quality in terms of readability for later time-steps as their layout depends on previous time-steps. Many experiments has been proposed to examine the effect of preserving the mental map in dynamic graphs visualization [START_REF] Purchase | How important is the "mental map"? -an empirical investigation of a dynamic graph layout algorithm[END_REF][START_REF] Saffrey | The "mental map" versus "static aesthetic" compromise in dynamic graphs: A user study[END_REF][START_REF] Purchase | Extremes are better: Investigating mental map preservation in dynamic graphs[END_REF]. The results of [START_REF] Purchase | Extremes are better: Investigating mental map preservation in dynamic graphs[END_REF] were quite surprising because the most effective visualizations were the extreme ones, i.e. the ones with very low or high mental map preservation: visualizations with medium preservation performed less well. The approach described in this paper aims to achieve high mental map preservation.

An interesting visualization dealing with dynamic large directed graphs has been proposed by [START_REF] Burch | Parallel edge splatting for scalable dynamic graph visualization[END_REF]. Vertices are ordered and positioned on several vertical parallel lines, and directed edges connect these vertices from left to right. Each time-step's graph is thus displayed between two consecutive vertical axes. [START_REF] Hu | Embedding, clustering and coloring for dynamic maps[END_REF] proposed a method based on a geographical metaphor to visualize clustered dynamic graphs. However, their approach requires one global clustering over time, while ours allows nodes to be transferred in order to create better local clusterings and to capture the evolutions of the communities.

Finding a partition of the nodes of a static graph according to its structure is a well studied problem; [START_REF] Schaeffer | Graph clustering[END_REF] has published a good overview of graph clustering methods. But clustering a dynamic graph is a less studied problem.

Clustering

A dynamic graph can be defined formally as an agglomerate graph G = (V, E) and an ordered sequence of subgraphs

S = {G 1 = (V 1 , E 1 ), G 2 = (V 2 , E 2 ), ..., G k = (V k , E k )}
where each G t is the subgraph of G at time t. V, V 1 , V 2 , ..., V k are finite and non-disjointed sets of nodes, E, E 1 , E 2 , ..., E k are finite and non-disjointed sets of edges such that

V = V 1 [ V 2 [ ... [ V k and E = E 1 [ E 2 [ ... [ E k .
What we need is to create a time-varying clustering, i.e. a set of clusters evolving over time. The clustering method we describe here is a two step algorithm. The first step consists of partitioning the nodes for each time step independently. Then, we associate these clusters through time to derive time-varying clusters.

Time-step Clusterings

To calculate a dynamic clustering, we first find a partition for each time step, i.e. a set of clusterings C = {C 1 , C 2 , ..., C k } where C t = {c t 1 , c t 2 , ..., c t lt } is a partition of the nodes V t of G t . In this paper, we call each C t a "time-step clustering" where c t i is the "time-step cluster" i at time t, and c t i ✓ V t for each i 2 J1, ltK, V t = c t 1 [ c t 2 [ ...c t l and c t i \ c t j = ; for each pair (i, j) 2 J1, ltK 2 . Our algorithm is based on the so-called modularity function of [START_REF] Newman | Graph clustering[END_REF]. It represents the sum of the number of edges linking nodes of the same clusters minus the expected such sum if edges were distributed at random. For a graph G t = (V t , E t ) and a partition C t of its nodes, the modularity Q(C t ) is defined by:

Q(C t ) = 1 2|E t | X u,v2Vt  A uv k u k v 2|E t | (c t (u), c t (v))
where |E t | is the number of edges, A uv is 1 if there is an edge between u and v and 0 otherwise,

k u = P v A uv is the number of edges attached to u, c t (u) is the time-step cluster of C t containing u, (c t (u), c t (v)) is 1 if c t (u) = c t (
v) and 0 otherwise. A partition that maximizes this function helps to discover clusters of densely connected communities. Moreover, as shown by [START_REF] Noack | Modularity clustering is force-directed layout[END_REF], optimizing the modularity is the same as optimizing an energy function in graph layout. This equivalence implies that our layout based on such a clustering algorithm yields a good representation of the graph.

The problem of finding a partition that maximizes the modularity is hard, and the corresponding decision problem is NP-complete [START_REF] Brandes | Maximizing modularity is hard[END_REF]. We use the heuristic proposed by [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF], which works well in terms of both the quality of the results and the computation time. Initially, each node belongs to its own cluster. Then pairs of clusters are recursively merged such that the modularity of the partitioning increases. If two possible merges involve the same cluster, the merge that improves the modularity the most is performed.

Time-varying Clustering

Overview

The previous approach [START_REF] Sallaberry | Clustering, visualizing, and navigating for large dynamic graphs[END_REF] was to compare each time-step cluster in the current time-step pairwise with the time-step clusters of the previous time-step according to the Jaccard index, and then to itratively and greedily associate the time-step clusters that most closely match into the same time-varying cluster, halting when all clusters are assigned or the similarity falls below a user-defined threshold. If there are any remaining new clusters that do not have a match, they are considered new clusters, and so they start new time-varying clusters. And any remaining time-varying clusters present in the previous time-step that were not assigned a cluster in the current time-step were discarded, as there was no match.

Our approach here operates fairly similarly, with two key differences. The first improvement is that rather than discarding time-varying clusters after they fail to find a match, we retain them in the matching algorithm so that they can return if that portion of the network returns to a similar enough state, even if there are many intervening time-steps. This reduces the overall number of time-varying clusters, and can aid in identifying certain patterns (such as periodicity, split/merges, clustering instability, etc...). The second improvement is to reduce the number of pairwise clusters to evaluate, by only comparing clusters that share at least one vertex. This optimization method has the added bonus that it can be used to calculate the Jaccard index very efficiently.

Details

We define a time-varying clustering of a dynamic graph G as a set of time-varying clusters

V C = {V C 1 , V C 2 , ..., V C l }. Each of these time-varying clusters is an ordered sequence V C i = {vc 1 i , vc 2 i , ..., vc k i }
where k is the number of time steps and each vc t i is a subset of the vertices V t at time t. That is, each time-varying cluster V C i is a cluster whose membership can evolve over time, where vc t i represents the set of nodes in the cluster i at time t. As the number of time-step clusters can change between timesteps, not every time-varying cluster is populated at every timestep, and the total number of time-varying clusters l can be larger than the number of time-step clusters at any time step.

We start from an empty set V C of time-varying clusters and we create a time-varying cluster V C i for each time-step cluster c 1 i of the first time-step clustering C 1 . The set of nodes of these time-varying clusters V C i at time 1 are initialized with the time-step clusters c 1 i :

vc 1 i c 1 i .
Then, for each subsequent timestep t, we want to compute similarities between each time-step cluster c t i 2 C t and potential time-varying clusters V C i . In our implementation, we use the Jaccard index to compute the similarities. For two clusters c t 1 i and c t j , this is defined by the equation

|c t 1 i \ c t j |/|c t 1 i [ c t j |.
There are two main advantages in using this metric.

First it takes into account the number of shared nodes as well as the total number of nodes, which guarantees homogeneity between consecutive steps of a time-varying cluster. Secondly it returns a value normalized between 0 and 1 which is helpful for empirically defining a threshold.

In the original algorithm, the association step was performed by comparing each time-step cluster c t i 2 C t to every time-step cluster c t 1 j 2 C t 1 to create a similarity matrix, which costs O(|C t | ⇤ |C t 1 |) times the cost of the similarity calculation. Associations are then performed greedily, starting with the largest matrix value, and stopping when either one set of clusters is exhausted or the remaining matrix values are less than the threshold. Then, any remaining clusters in C t were assigned new time-varying clusters in V C. And any remaining clusters in C t 1 were discarded, and their corresponding time-varying cluster in V C was terminated for the remainder of the execution.

However, we found that this process led to the creation of many new time-varying clusters, because sometimes the network might revert to a clustering previously encountered where the corresponding time-varying cluster was already terminated, so the method would create a new time-varying cluster. To resolve this, we preserve the most recent time-step cluster c u k for each time-varing clusters in V C that would have been terminated, and compare against those as well, where 0 < u < t 1. If an older c u k is more similar than any time-step cluster c t 1 j 2 C t 1 , then the system will select and revive the time-varying cluster of c u k instead of that of one of active time-varying clusters V C i . This adds additional complexity to the algorithm, but produces more robust and compact results, as time-varying clusters are allowed to reform instead being terminated.

Another issue is that many of the time-step clusters are disjoint sets of nodes, and thus we do not need to compute their similarity. In our improved version, rather than computing the entire matrix, we consider only the clusters that share at least one node. We do this by computing a list of candidate clusters CC t i for each time-step cluster c t i , which we define as

CC t i = {(vc ⇤ a , |c t i \ vc ⇤ a |), (vc ⇤ b , |c t i \ vc ⇤ b |), ...}
, where each vc ⇤ j is the most recent timestep of V C j . This can be computed relatively efficiently by iterating over each node n 2 c t i : for each n we take any vc ⇤ j that contains n, then either add vc ⇤ j to CC t i with a paired value of 1, or if it is already in the list, we simply increment its paired value. To make this process even more efficient, we use a lookup table to map each node n to its existing vc ⇤ j clusters. In the original process, this would be trivial to do, as each node n could only exist in one previous cluster vc ⇤ j . However, since we are now preserving terminated clusters, it is possible for n to have multiple previous clusters. So we build a hash to map each node n to its prior clusters by iterating over each vc ⇤ j and inserting a pointer to vc ⇤ j for each node m 2 vc ⇤ j . This elimates the need for any searching. While this process appears to add some computational overhead as it iterates over every node instead of working with the clusters, this is entirely offset in the calculation of the Jaccard index, as we have already computed the size of each intersection |c t i \ vc ⇤ j |. From this, we can calculate the Jaccard index in O(1) time as J(c t i , vc ⇤ j ) =

|c t i \vc ⇤ j | |c t i |+|vc ⇤ j | |c t i \vc ⇤ j |
. So while this optimization does add additional memory overhead, the computation is much more efficient.

Our visualization is based not just on a clustering, but on an ordering of the time-varying clusters (in the next sections, we use the word "cluster" instead of "time-varying cluster" to simplify the notation). Then, nodes are also ordered within each cluster. A node that moves from a cluster V C a to a cluster V C b is involved in the node ordering of both V C a and V C b , e.g. it can be the 6th node of V C a and the 3rd node of V C b .

Ordering clusters

The stability of the layout is one of the main goals of our method: we want to easily see the evolution of the clusters and also be able to follow nodes that move between clusters. As the layout depends on the ordering, clusters need to be ordered in such a way that two clusters exchanging many nodes are close to each other.

We do this by first creating a weighted quotient graph QG = (V QG , E QG , !) defined by the relationships between the time-varying clusters V C of G. Each node of V QG represents a cluster of V C, i.e. V QG V C. There is an edge in E QG between V C i and V C j if and only if there is at least one node in the sets of V C i that is also in a set of V C j . The weight function ! is a function ! : E QG ! N defined for each edge e = (V C i , V C j ) as the number of transferred of nodes between sets of V C i and sets of V C j .

Next we need to find an ordering of these clusters, i.e. a permutation : V QG ! {1, 2, ..., |V QG |} that minimizes the function:

LA (QG) = X uv2EQG u,v2VQG !(uv) • | (u) (v)|
This function is called the Linear arrangement function (LA) and finding an ordering that minimizes it is known as the Minimum Linear Arrangement Problem, MinLA [START_REF] Petit | Experiments on the minimum linear arrangement problem[END_REF].

MinLA is NP-hard and the corresponding decision problem is NP-complete [START_REF] Garey | Dynamic graph drawing of sequences of orthogonal and hierarchical graphs[END_REF]. Many heuristics have been proposed to find a satisfying solution. A list of these methods and an experiment has been proposed by [START_REF] Petit | Experiments on the minimum linear arrangement problem[END_REF]. More recently, [START_REF] Koren | A multi-scale algorithm for the linear arrangement problem[END_REF] have proposed a new heuristic that is a good compromise between computation time and quality of the results.

Ordering nodes

The second ordering step consists in finding a permutation of the nodes within each cluster V C i of V C. Since we want to maximize stability, we calculate this permutation over all time, so that nodes will not move within a cluster, even if this leaves gaps at some time-steps. As with the clusters, this is another MinLA problem. Let vc i be the set of nodes of V C i : vc i = S 1tk vc t i . Then the permutation is defined as ' i : vc i ! {1, 2, ..., |vc i |}. This ordering needs to take into account the ordering of the clusters computed previously: for example, if a node v moves only once from a cluster V C a to a cluster V C b and if (V C a ) < (V C b ), then v should lie at the upper extremity of V C a (high ' a (v)) and at the lower extremity of V C b (low ' b (v)). To find the permutation ' i , we first compute for each node v of vc i the median of the clusters V C i v belongs to:

median i (v) = X V Cj ;v2vcj (V C j ) |{V C j ; v 2 vc j }|
Then, the permutation ' i is the ordering obtained by sorting the nodes of vc i according to their median value:

median i (v) < median i (u) , ' i (v) < ' i (u).

Visualization

The visualization methods we employ focus on representing the evolving clusters in dynamic graphs. We employ two views: a time-line inspired by [START_REF] Ogawa | Software evolution storylines[END_REF]; [START_REF] Tanahashi | Design considerations for optimizing storyline visualizations[END_REF] that provides an overview of the entire dynamic graph, and a more traditional node-link view for individual time-steps. Both of these views are derived from the clustering and ordering methods described earlier. Moreover, since the clustering and ordering are computed as a preprocessing step, the computation times of the visualizations are linear, which makes it possible to obtain real-time, interactive navigation of the dynamic graph.

Time-line view

The time-line view depicts an overview of the nodes' arrangement into clusters and of the nodes motion between clusters. Each node is represented as a line where the x-position is time and the y-positions corresponds the cluster the nodes belong to at each given time and its position within the cluster. Figure 1(a) shows an example. For reference purposes, in this diagram the time-steps are represented with vertical grey lines (from t=1 to 5) positioned along the x-axis that represents time. There are 8 plotted lines (including black, blue and green ones), which correspond to 8 nodes. There are four clusters on the y-axis, and a horizontal line is in front of one of them when the corresponding node belongs to it. Clusters are positioned according to the ordering computed by the pre-processing algorithm, from bottom to top (e.g. the cluster labelled 4 is the cluster V C a with a such that 4 = (V C a )).

As an example of reading this plot, consider the blue line. The corresponding node v belongs to the cluster 4 at times 1 and 2, and it belongs to the cluster 3 at times 3 and 4, since the blue line moves from cluster 4 to cluster 3 at time 3. Also, the blue node is no longer in the graph at the time-step 5 and that the green node appears in the graph at the time-step 2.

Lines in the clusters are positioned according to the orderings ' i computed during the preprocessing step. In this way, a node v that moves from a cluster V C a to a cluster V C b with

(V C a ) < (V C b
) is likely to be positioned at the upper extremity of V C a (high ' a (v)) and at the lower extremity of V C b (low ' b (v)). This technique reduces edge crossings and improves the readability of the view.

Clusters are separated by a constant gap to clarify their distinctions. The height of a cluster V C i corresponds to the size the set vc i of all the nodes that belong to it at least for one timestep. As an example, |vc b | = 4 (see the red circle 1) and |vc c | = 3 with c such that 2 = (V C c ) (see the red circles 2 and 3). Thanks to this, a node has always the same position in the same cluster, so there will be no bends when a node remains in the same cluster and the area devoted to a cluster remains the same. FIG. 1: (a) The time-line gives an overview of the clusters and of the nodes moving from clusters to clusters. Each horizontal or bent line is a node. Vertical grey lines represent timesteps, from 1 to 5. Y-axis represents clusters, e.g. the blue line near the cluster 4 at time-steps 1 and 2 stands for a node v that belongs to V C a with a such that 4 = (V C a ), at the time 1 and at the time 2 (it belongs to vc 1 a and vc 2 a but not to vc 3 a , vc 4 a and vc 5 a ). (b) Time-step view of the graph used in the example of Figure 1(a). It shows the graph at the 3rd time-step. We don't display the edges here. Nodes represented as disks are positioned along a Hilbert's curve, represented by the black bended line.

Time-step view

The second view is a node-link diagram that shows the graph at any selected time-step. The layout is based on the technique of [START_REF] Muelder | Rapid graph layout using space filling curves[END_REF], which maps a 1-D ordering of nodes to a space-filling curve to define the layout.

Since we have already computed a stable ordering of nodes, it is sensible to map this same ordering onto the space-filling curve. In the timeline, the height of each line at any time step corresponds to the pre-computed ordering. So, we can reuse these y-positions as a 1-D layout for that time-step, then map the nodes directly to a space-filling curve by placing the nodes at the corresponding distance along the curve. This is done by normalizing both the 1-D layout and the length of the curve, then calculating the position of each node by recursively mapping it to the curve in constant time, as in the original paper [START_REF] Muelder | Rapid graph layout using space filling curves[END_REF]. Figure 1(b) shows an example of this node positioning on the same example as the one presented in Figure 1(a) for time-step 3. In this diagram, we use a Hilbert curve, but we also use Peano curve, a Gosper curve, and an H-curve, and the user can switch between these curves as desired (see [START_REF] Haverkort | Locality and bounding-box quality of twodimensional space-filling curves[END_REF] for a summary of well-known space-filling curves).

One interesting property of a space-filling curve is known as the Worst-Case Locality [START_REF] Haverkort | Locality and bounding-box quality of twodimensional space-filling curves[END_REF]. This property guarantees that the euclidean distances between nodes in the layout are bounded by the distances of the same nodes in the onedimensional layout. So, the proximities of elements (nodes/clusters) depend directly on the ordering. As the ordering is based primarily on the connectivity of the networks, this guarantees layout quality metrics, such as tightly connected groups of nodes being placed close together with a good aspect ratio, and short average edge lengths.

Since a node has always the same position in the time-line when it is in the same cluster and the area devoted to a cluster remains the same, its placement in the layout will also be constant. This ensures the stability of the layout. Even the distance that nodes move is minimized, as the ordering is such that clusters that exchange many nodes are placed closer together.

As the layout itself runs in linear time, the visualization can be updated interactively by the user and we can even easily play the sequence of graphs and animate the transitions with graphs of tens of thousands nodes/edges (see our previous paper [START_REF] Sallaberry | Clustering, visualizing, and navigating for large dynamic graphs[END_REF] for more details).

As we use a clustering hierarchy, we can also employ the hierarchical edge bundling technique of [START_REF] Holten | Hierarchical edge bundles: Visualization of adjacency relations in hierarchical data[END_REF] which improves the readability of the graph. Control points of the spline linking a node v and a node u are defined by the path through the clustering hierarchy, and placed according to the clusters' centroids.

Interaction and navigation

One of the most useful features of our approach is that any time-step can be laid out quickly and directly, without needing to iterate over the other time-steps. The benefit of this is that it enables random access. That is, users can find interesting time-steps in the time-line and skip between them directly. We enable this form of interaction by letting the user simply click in the time-line on the time-step that they want to load. We also include the more traditional approach of simply animating over the entire dynamic graph. In either case, the positions of nodes that move are interpolated between time-steps so that the user can follow their motion. Within the node-link diagram itself, we can also allow for traditional graph interaction, such as selection, or focus+context zooming.

Discussion

We have described two algorithmic improvements over the previous approach. First, our new approach is tuned to reduce the total number of persistant time-varying clusters, which will improve the space utilization of the resulting visualization. Second, we have substantially optimized the cluster association step, which greatly reduces the time it takes to process the data. Here, we present some case studies and quantitatively evaluate both improvements.

To test our improvements, we ran both improvements as well as the original approach on several datasets. First, we ran them on a social network dataset collected from the Rimzu social networking site, as used by [START_REF] Frishman | Online dynamic graph drawing[END_REF]. While quite small and straightforward, this dataset makes a good baseline for comparison against existing works. Next, we evaluated the MIT Reality dataset [START_REF] Eagle | CRAWDAD data set mit/reality (v. 2005-07-01[END_REF]. In this dataset, there is a small, strong core of about 80 people, but hundreds of periferal nodes that do persist through the data and which generally only have one connection. Due to their transitivity, they do not contriubte much to the overall structure of the network, so we evaluate both the entire network and just the core network without the external nodes. Finally, we evaluate the improvement on a dataset of the autonomous systems of the internet, derived from data collected by the Oregon Route Views project and as used in several existing works [START_REF] Muelder | Advanced Visualization Techniques for Abstract Graphs and Computer Networks[END_REF][START_REF] Sallaberry | Clustering, visualizing, and navigating for large dynamic graphs[END_REF].

In each case, the results followed our expectation, as is shown in Table 1. The inclusion of additional comparisons against dead clusters succeeds in reducing the total number of dynamic time-varying clusters, but requires more comparisons, and hence takes more computation time. However, this is entirely offset by the optimization, which reduces the computation time by around a factor of 20 or more. TAB. 1: Results of our improvements. We found that our methods reduce the number of timevarying clusters by up to half, and reduce computation time by a factor of about 20 or more.

Original

Conclusion and Future Work

This paper presents an incremental improvement to a previous work. While the visual representation remains the same, the underlying clustering algorithm has been greatly improved and optimized. The evaluation follows accordingly, and demonstrates our improvements quantitatively.

While the results described in this paper are promising, they are still further ways they can be improved. Most importantly, the cluster and node ordering steps are still extremely slow and we are investigating ways to improve upon this. One such method is to assign more localized orderings, where we find the actual time range used by each dynamic cluster, and more compactly arrange the clusters accordingly, similar to the work of [START_REF] Tanahashi | Design considerations for optimizing storyline visualizations[END_REF]. However, this algorithm is also very time consuming, and we are investigating heuristic improvements.


Ioan Doré Landau
email: landau@lag.ensieg.inpg.fr

Aurelian Constantinescu

Adaptive rejection of unknown finite band disturbances -Theory and applications

It is a great pleasure for us to contribute to this symposium dedicated to Masayoshi Tomizuka on the occasion of his sixtieth birthday. The important contributions of Masayoshi Tomizuka to advanced control of mechanical systems have been largely recognized by the control community. In particular he has dedicated a lot of research efforts to "adaptive control" and "repetitive control" (for rejection of periodic disturbances). This paper also discusses an adaptive approach to disturbance rejection but for a different class of disturbances. Both Tomizuka's approach to disturbance rejection and ours share some common grounds. We do expect to explore jointly in the future the connections between the two approaches.

Introduction

One of the basic problems in control is the attenuation (rejection) of unknown disturbances without measuring them. The common framework is the assumption that the disturbance is the result of a white noise or a Dirac impulse passed through the "model of the disturbance". While in general one can assume a certain structure for such "model of disturbance", its parameters are unknown and may be time varying. This will require to use an adaptive approach. To be more specific, the disturbances considered can be defined as "finite band disturbances". This includes single or multiple narrow band disturbances or sinusoidal disturbances. Furthermore for robustness reasons the disturbances should be located in the frequency domain within the regions where the plant has enough gain (see explanation in section 3).

Solutions for this problem, provided that an "image" of the disturbance can be obtained by using an additional transducer, have been proposed by the signal processing community and a number of applications are reported ([?, ?, ?, ?]). However, these solutions (inspired by Widrow's technique for adaptive noise cancellation ([?])) ignore the possibilities offered by feedback control systems and require an additional transducer. The principle of this signal processing solution for adaptive rejection of unknown disturbances is illustrated in figure 1, where E{y 2 } represents the variance of the output y. The basic idea is that a "well located" transducer can provide a measurement, highly correlated with the unknown disturbance. This information is applied to the control input of the plant through an adaptive filter (in general a Finite Impulse Response -FIR) whose parameters are adapted such that the effect of the disturbance upon the output is minimized. The disadvantages of this approach are:

• It requires the use of an additional transducer.

• Difficult choice for the location of this transducer (it is probably the main disadvantage). • It requires the adaptation of many parameters. To achieve the rejection of the disturbance (at least asymptotically) without measuring it, a feedback solution can be considered. As mentioned earlier, the common framework is the assumption that the disturbance is the result of a white noise or a Dirac impulse passed through the "model of the dis-turbance"3 . Several problems have been considered within this framework leading to adaptive feedback control solutions:

1. Unknown plant and disturbance models ([?]). 2. Unknown plant model and known disturbance ([?, ?]). 3. Known plant and unknown disturbance model ([?, ?, ?, ?, ?, ?, ?, ?, ?]).

The present paper will focus on the last case, since this is the situation encountered in many applications. Among the various approaches considered for solving this problem, the following ones may be mentioned:

1. Use of the internal model principle ([?, ?, ?, ?, ?, ?, ?, ?, ?, ?]). 2. Use of an observer for the disturbance ([?, ?]). 3. Use of the "phase-locked" loop structure considered in communication systems ([?, ?]).

Of course, since the parameters of the disturbance model are unknown, all these approaches lead to an adaptive implementation which can be of direct or indirect type.

From the user point of view and taking into account the type of operation of existing adaptive disturbance compensation systems one has to consider two modes of operation of the adaptive schemes:

• Self-tuning operation (the adaptation procedure starts either on demand or when the performance is unsatisfactory and the current controller is frozen during the estimation/computation of the new controller parameters). • Adaptive operation (the adaptation is performed continuously and the controller is updated at each sampling).

This paper explores the use of the internal model principle for the rejection of unknown time-varying finite band disturbances. The other related approaches will be briefly reviewed in Section 6.

Using the internal model principle, the controller should incorporate the model of the disturbance ([?, ?, ?, ?]). Therefore the rejection of unknown disturbances raises the problem of adapting the internal model of the controller and its re-design in real-time.

One way for solving this problem is to try to estimate in real time the model of the disturbance and re-compute the controller, which will incorporate the estimated model of the disturbance (as a pre-specified element of the controller). While the disturbance is unknown and its model needs to be estimated, one assumes that the model of the plant is known (obtained for example by identification). The estimation of the disturbance model can be done by using standard parameter estimation algorithms (see for example [?, ?]). This will lead to an indirect adaptive control scheme. The principle of such a scheme is illustrated in figure 2. The time consuming part of this approach is the redesign of the controller at each sampling time. The reason is that in many applications the plant model can be of very high dimension and despite that this model is constant, one has to re-compute the controller because a new internal model should be considered. This approach has been investigated in [?, ?, ?]. However, by considering the Youla-Kucera parametrization of the controller (known also as the Q-parametrization), it is possible to insert and adjust the internal model in the controller by adjusting the parameters of the Q polynomial (see figure 3). It comes out that in the presence of unknown disturbances it is possible to build a direct adaptive control scheme where the parameters of the Q polynomial are directly adapted in order to have the desired internal model without recomputing the controller (polynomials R 0 and S 0 in figure 3 remain unchanged). The number of the controller parameters to be directly adapted is roughly equal to the number of parameters of the denominator of the disturbance model. In other words, the size of the adaptation algorithm will depend upon the complexity of the disturbance model.

This paper focuses on the direct feedback adaptive control for the case of unknown and time-varying frequency narrow band disturbances. The direct adaptive control scheme to be presented([?]) takes advantage of the Youla-Kucera parametrization for the computation of the controller. This algorithm takes its roots from an idea of Tsypkin ([?]) [START_REF] Anderson | From Youla-Kucera to identification, adaptive and nonlinear control[END_REF] . A similar approach has been considered in [?] for an application to a chemical reactor but a theoretical analysis of the scheme is not provided. A related paper is [?, ?] where the application field is the active noise control in an acoustic duct [START_REF] Bengtsson | Output regulation and internal models -a frequency domain approach[END_REF] . For evaluation purposes (complexity and performance) an indirect adaptive control scheme based on the Internal Model Principle will be also presented.

The paper is organized as follows. Section 2 is dedicated to a brief review of the plant, disturbance and controller representation as well as of the Internal Model Principle. Some robustness issues are addressed in section 3. The direct and the indirect adaptive control schemes for disturbance rejection are presented in sections 4 and 5, respectively. Other approaches proposed in the literature are discussed in section 6. The active suspension on which the proposed algorithms have been tested, as well as the real-time results obtained are presented in section 7. Some concluding remarks are given in section 8.

Plant representation and controller structure

The structure of a linear time invariant discrete time model of the plant (used for controller design) is:

G(z -1 ) = z -d B(z -1 ) A(z -1 ) = z -d-1 B * (z -1 ) A(z -1 ) , (1)
with: d = the plant pure time delay in number of sampling periods ;

A = 1 + a 1 z -1 + . . . + a nA z -nA ; B = b 1 z -1 + . . . + b nB z -nB = q -1 B * ; B * = b 1 + . . . + b nB z -nB +1 ,
where A(z -1 ), B(z -1 ), B * (z -1 ) are polynomials in the complex variable z -1 and n A , n B and n B -1 represent their orders [START_REF] Beranek | Noise and Vibration Control Engineering: Principles and Applications[END_REF] . The model of the plant may be obtained by system identification. Details on system identification of the models considered in this paper can be found in [?, ?, ?, ?, ?, ?].

Since in this paper we are focused on regulation, the controller to be designed is a RS-type polynomial controller ([?, ?]) -see also figure 8.

The output of the plant y(t) and the input u(t) may be written as:

y(t) = q -d B(q -1 ) A(q -1 ) • u(t) + p 1 (t) ; (2)
S(q -1 ) • u(t) = -R(q -1 ) • y(t) , (3)
where q -1 is the delay (shift) operator (x(t) = q -1 x(t + 1)) and p 1 (t) is the resulting additive disturbance on the output of the system. R(z -1 ) and S(z -1 ) are polynomials in z -1 having the orders n R and n S , respectively, with the following expressions:

R(z -1 ) = r 0 + r 1 z -1 + . . . + r nR z -nR = R ′ (z -1 ) • H R (z -1 ) ; (4)
S(z -1 ) = 1 + s 1 z -1 + . . . + s nS z -nS = S ′ (z -1 ) • H S (z -1 ) , (5)
where H R and H S are pre-specified parts of the controller (used for example to incorporate the internal model of a disturbance or to open the loop at certain frequencies).

We define the following sensitivity functions:

• Output sensitivity function (the transfer function between the disturbance p 1 (t) and the output of the system y(t)):

S yp (z -1 ) = A(z -1 )S(z -1 ) P (z -1 ) ;

• Input sensitivity function (the transfer function between the disturbance p 1 (t) and the input of the system u(t)):

S up (z -1 ) = - A(z -1 )R(z -1 ) P (z -1 ) , (7)
where

P (z -1 ) = A(z -1 )S(z -1 ) + z -d B(z -1 )R(z -1 ) = A(z -1 )S ′ (z -1 ) • H S (z -1 ) + z -d B(z -1 )R ′ (z -1 ) • H R (z -1 ) (8)
defines the poles of the closed loop. In pole placement design, P (z -1 ) is the polynomial specifying the desired closed loop poles and the controller polynomials R(z -1 ) and S(z -1 ) are minimal degree solutions of (8) where the degrees of P , R and S are given by: n

P ≤ n A + n B + d -1, n S = n B + d -1 and n R = n A -1.
Using the equations ( 2) and (3), one can write the output of the system as:

y(t) = A(q -1 )S(q -1 ) P (q -1 ) • p 1 (t) = S yp (q -1 ) • p 1 (t) . (9)
For more details on RS-type controllers and sensitivity functions see [?]. Suppose that p 1 (t) is a deterministic disturbance, so it can be written as

p 1 (t) = N p (q -1 ) D p (q -1 ) • δ(t) , (10)
where δ(t) is a Dirac impulse and N p (z -1 ), D p (z -1 ) are coprime polynomials in z -1 , of degrees n Np and n Dp , respectively. In the case of stationary disturbances the roots of D p (z -1 ) are on the unit circle. The energy of the disturbance is essentially represented by D p . The contribution of the terms of N p is weak compared to the effect of D p , so one can neglect the effect of N p . Internal Model Principle: The effect of the disturbance given in (10) upon the output:

y(t) = A(q -1 )S(q -1 ) P (q -1 ) • N p (q -1 ) D p (q -1 ) • δ(t) , (11)
where D p (z -1 ) is a polynomial with roots on the unit circle and P (z -1 ) is an asymptotically stable polynomial, converges asymptotically towards zero if and only if the polynomial S(z -1 ) in the RS controller has the form:

S(z -1 ) = D p (z -1 )S ′ (z -1 ) . (12)
In other terms, the pre-specified part of S(z -1 ) should be chosen as H S (z -1 ) = D p (z -1 ) and the controller is computed using (8), where P , D p , A, B, H R and d are given [START_REF] Bodson | Rejection of periodic distrubances of unknown and time-varying frequency[END_REF] .

Using the Youla-Kucera parametrization (Q-parametrization) of all stable controllers ([?, ?]), the controller polynomials R(z -1 ) and S(z -1 ) get the form:

R(z -1 ) = R 0 (z -1 ) + A(z -1 )Q(z -1 ) ; (13)
S(z -1 ) = S 0 (z -1 ) -z -d B(z -1 )Q(z -1 ) . (14)
The (central) controller (R 0 , S 0 ) can be computed by poles placement (but any other design technique can be used). Given the plant model (A, B, d) and the desired closed-loop poles P one has to solve:

P (z -1 ) = A(z -1 )S 0 (z -1 ) + z -d B(z -1 )R 0 (z -1 ) . (15)
Equations ( 13) and ( 14) characterize the set of all stabilizable controllers assigning the closed loop poles as defined by P (z -1 ) (it can be verified by simple calculations that the poles of the closed loop remain unchanged). For the purpose of this paper Q(z -1 ) is considered to be a polynomial of the form:

Q(z -1 ) = q 0 + q 1 z -1 + . . . + q nQ z -nQ . ( 16
)
To compute Q(z -1 ) in order that the controller incorporates the internal model of the disturbance one has to solve the diophantine equation:

S ′ (z -1 )D p (z -1 ) + z -d B(z -1 )Q(z -1 ) = S 0 (z -1 ) , (17)
where D p (z -1 ), d, B(z -1 ) and S 0 (z -1 ) are known and S ′ (z -1 ) and Q(z -1 ) are unknown. Equation ( 17) has a unique solution for S ′ (z -1 ) et Q(z -1 ) with:

n S0 ≤ n Dp + n B + d -1, n S ′ = n B + d -1, n Q = n Dp -1 .
One sees that the order n Q of the polynomial Q depends upon the structure of the disturbance model.

Robustness considerations

As it is well known, the introduction of the internal model for the perfect rejection of the disturbance (asymptotically) will have as effect to raise the maximum value of the modulus of the output sensitivity function S yp . This may lead to unacceptable values for the modulus and the delay margins if the controller design is not appropriately done (see [?]). As a consequence, a robust control design should be considered assuming that the model of the disturbance is known, in order to be sure that for all situations an acceptable modulus margin and delay margin are obtained.

On the other hand at the frequencies where perfect rejection of the disturbance is achieved one has S yp (e -jω ) = 0 and

S up (e -jω ) = A(e -jω ) B(e -jω ) . (18)
Equation ( 18) corresponds to the inverse of the gain of the system to be controlled. The implication of equation ( 18) is that cancellation (or in general an important attenuation) of disturbances on the output should be done only in frequency regions where the system gain is large enough. If the gain of the controlled system is too low, |S up | will be large at these frequencies. Therefore, the robustness vs additive plant model uncertainties will be reduced and the stress on the actuator will become important. Equation ( 18) also implies that serious problems will occur if B(z -1 ) has complex zeros close to the unit circle (stable or unstable zeros) at frequencies where an important attenuation of disturbances is required. It is mandatory to avoid attenuation of disturbances at these frequencies.

Since on one hand we would not like to react to very high frequency disturbances and on the other hand we would like to have a good robustness it is often wise to open the loop at 0.5f s (f s is the sampling frequency) by introducing a fixed part in the controller H R (q -1 ) = 1 + q -1 (for details see [?] and section 2).

Direct adaptive control for disturbance attenuation

The objective is to find an estimation algorithm which will directly estimate the parameters of the internal model in the controller in the presence of an unknown disturbance (but of known structure) without modifying the closed loop poles. Clearly, the Q-parametrization is a potential option since modifications of the Q polynomial will not affect the closed loop poles. In order to build an estimation algorithm it is necessary to define an error equation which will reflect the difference between the optimal Q polynomial and its current value.

In [?], such an error equation is provided and it can be used for developing a direct adaptive control scheme. This idea has been used in [?, ?, ?, ?]. Using the Q-parametrization, the output of the system in the presence of a disturbance can be expressed as:

y(t) = A(q -1 )[S 0 (q -1 ) -q -d B(q -1 )Q(q -1 )] P (q -1 ) • N p (q -1 ) D p (q -1 ) • δ(t) = S 0 (q -1 ) -q -d B(q -1 )Q(q -1 ) P (q -1 ) • w(t) , (19)
where w(t) is given by (see also figure 3):

w(t) = A(q -1 )N p (q -1 ) D p (q -1 ) • δ(t) = A(q -1 ) • y(t) -q -d • B(q -1 ) • u(t) . ( 20
)
In the time domain, the internal model principle can be interpreted as finding Q such that asymptotically y(t) becomes zero. Assume that one has an estimation of Q(q -1 ) at instant t, denoted Q(t, q -1 ). Define ε 0 (t + 1) as the value of y(t + 1) obtained with Q(t, q -1 ). Using (19) one gets:

ε 0 (t + 1) = S 0 (q -1 ) P (q -1 ) • w(t + 1) - q -d B * (q -1 ) P (q -1 ) Q(t, q -1 ) • w(t) . (21)
One can define now the a posteriori error (using Q(t + 1, q -1 )) as:

ε(t + 1) = S 0 (q -1 ) P (q -1 ) • w(t + 1) - q -d B * (q -1 ) P (q -1 ) Q(t + 1, q -1 ) • w(t) . (22)
Replacing S 0 (q -1 ) from the last equation by ( 17) one obtains

ε(t + 1) = [Q(q -1 ) -Q(t + 1, q -1 )] • q -d B * (q -1 ) P (q -1 ) • w(t) + v(t + 1) , (23)
where v(t) = S ′ (q -1 )D p (q -1 ) P (q -1 )

• w(t) = S ′ (q -1 )A(q -1 )N p (q -1 )

P (q -1 ) • δ(t)
is a signal which tends asymptotically towards zero. Define the estimated polynomial Q(t, q -1 ) as: Q(t, q -1 ) = q0 (t)+q 1 (t)q -1 + . . . + qnQ (t)q -nQ and the associated estimated parameter vector : θ(t) = [q 0 (t) q1 (t) . . . qnQ (t)] T . Define the fixed parameter vector corresponding to the optimal value of the polynomial Q as: θ = [q 0 q 1 . . . q n Q ] T . Denote:

w 2 (t) = q -d B * (q -1 ) P (q -1 ) • w(t) (24)
and define the following observation vector:

φ T (t) = [w 2 (t) w 2 (t -1) . . . w 2 (t -n Q )] . (25)
Equation ( 23) becomes

ε(t + 1) = [θ T -θT (t + 1)] • φ(t) + v(t + 1) . (26)
One can remark that ε(t) corresponds to an adaptation error ([?]). From equation ( 21) one obtains the a priori adaptation error:

ε 0 (t + 1) = w 1 (t + 1) -θT (t)φ(t) ,
with w 1 (t + 1) = S 0 (q -1 ) P (q -1 )

• w(t + 1) ;

w 2 (t) = q -d B * (q -1 ) P (q -1 ) • w(t) ; (27)
w(t + 1) = A(q -1 ) • y(t + 1) -q -d B * (q -1 ) • u(t) , (28)
where B(q -1 )u(t + 1) = B * (q -1 )u(t). The a posteriori adaptation error is obtained from (22):

ε(t + 1) = w 1 (t + 1) -θT (t + 1)φ(t) .
For the estimation of the parameters of Q(t, q -1 ) the following parameter adaptation algorithm is used ([?]):

θ(t + 1) = θ(t) + F (t)φ(t)ε(t + 1) ; (30) ε(t + 1) = ε 0 (t + 1) 1 + φ T (t)F (t)φ(t) ; (31)
ε 0 (t + 1) = w 1 (t + 1) -θT (t)φ(t) ; (32)
F (t + 1) = 1 λ 1 (t)   F (t) - F (t)φ(t)φ T (t)F (t) λ1(t) C) + φ T (t)F (t)φ(t)   . ( 33
) 1 ≥ λ 1 (t) > 0; 0 ≤ λ 2 (t) < 2 (34)
where λ 1 (t), λ 2 (t) allow to obtain various profiles for the evolution of the adaption gain F (t) (for details see [?, ?] and section 7).

In order to implement this methodology for disturbance rejection (see fig-

ure 3), it is supposed that the plant model

z -d B(z -1 ) A(z -1 ) is known (identified)
and that it exists a controller [R 0 (z -1 ), S 0 (z -1 )] which verifies the desired specifications in the absence of the disturbance. One also supposes that the degree n Q of the polynomial Q(z -1 ) is fixed, n Q = n Dp -1, i.e. the structure of the disturbance is known.

The following procedure is applied at each sampling time for adaptive operation:

1. Get the measured output y(t + 1) and the applied control u(t) to compute w(t + 1) using ( 29). 2. Compute w 1 (t + 1) and w 2 (t) using ( 27) and ( 28) with P given by (15). 3. Estimate the Q-polynomial using the parametric adaptation algorithm ( 30) -( 33). 4. Compute and apply the control (see figure 3):

S 0 (q -1 ) • u(t + 1) = -R 0 (q -1 ) • y(t + 1) -Q(t, q -1 ) • w(t + 1) . (35)
For the self tuning operation of the adaptive scheme, the estimation of the Q-polynomial starts once the level of the output is over a defined threshold. A parameter adaptation algorithm (30)- (33) with decreasing adaption gain is used and the estimation is stopped when the adaption gain is below a prespecified level [START_REF] Bodson | Adaptive algorithms for the rejection of sinusosidal disturbances with unknown frequency[END_REF] . During estimation of the new parameters, the controller is kept constant. The controller is updated once the estimation phase is finished. For a stability analysis of this scheme see [?].

Indirect adaptive control for disturbance attenuation

The methodology proposed in this section concerns the indirect adaptive control for the attenuation of unknown disturbances and consists in two steps: (1) Identification of the disturbance model; (2) Computation of a digital controller using the identified disturbance model.

The disturbance is considered as a stationary signal having a rational spectrum. As such it may be considered as the output of a filter with the transfer function N p (z -1 )/D p (z -1 ) and a white noise as input:

D p (q -1 ) • y(t) = N p (q -1 ) • e(t) or y(t) = N p (q -1 ) D p (q -1 ) • e(t) , (36)
where e(t) represents a gaussian white noise, and

N p (z -1 ) = 1 + n p1 z -1 + . . . + n pn Np z -nN p = 1 + z -1 N * p (z -1 ) ; D p (z -1 ) = 1 + d p1 z -1 + . . . + d pn Dp z -nD p = 1 + z -1 D * p (z -1 )
.

Therefore the disturbance model can be represented by an ARMA model.

As we deal with narrow band disturbances, the filtering effect of the primary path in cascade with the output sensitivity function (when operating in closed loop) around the central frequency of the disturbance can be approximated by a gain and a phase lag which will be captured by the

N p (z -1 ) D p (z -1 ) model.
From equation (36) one obtains:

y(t + 1) = - nD p i=1 d pi y(t -i + 1) + nN p i=1
n pi e(t -i + 1) + e(t + 1) .

The problem is, in fact, an on-line adaptive estimation of parameters in presence of noise ([?, ?]). Equation ( 37) is a particular case of identification of an ARMAX model. One can use for example the Recursive Extended Least Squares method ([?, ?]), which is dedicated to the identification of this type of model. The parameter adaptation algorithm given in ( 30) -( 33) is used.

The controller parameters are frozen while the disturbance model is identified. Once the disturbance model is identified, the controller containing the disturbance dynamics is computed by solving the diophantine equation (8) and using (5) with H S (z -1 ) = Dp (z -1 ) (the identified model of the disturbance). In order to apply this methodology we suppose that the plant model is known (can be obtained by identification). We also suppose that the degrees n Np and n Dp of N p (z -1 ) respectively D p (z -1 ) are fixed.

For the self tuning operation a parameter adaptation algorithm with decreasing adaptation gain is used and estimation is stopped when the adaptation gain becomes smaller than a pre-specified level. In adaptive operation, a parameter adaption algorithm with constant trace or forgetting factor is used. The parameters of the controller have to be re-computed at each sampling instant based on the current estimation of the disturbance model (non vanishing adaptation gain). Unfortunately, when the estimated frequency approaches the true one, the output sensitivity function will have a pair of complex zeros on the unit circle leading to a very strong attenuation of the measured effect of the disturbance. This will make the estimation of the exact frequency almost impossible. As a consequence, in adaptive operation there will be a "bias" on the estimated frequency caused by the need to have a certain level of the measured output to carry on the estimation. This phenomenon has been clearly observed both in simulation and on the real system. [START_REF] Constantinescu | Commande robuste et adaptative d'une suspension active[END_REF] For more details see([?, ?, ?]).

To compute the controller one uses equation 8. If a Youla-Kucera parametrization of the controller is used one will have to solve equation 17. The advantage with respect to the first approach is a reduction of the size of the equation to be solved in real time.

Review of other approaches

Several approaches based on adaptive feedback control have been considered in the literature for the rejection of unknown disturbances when the plant model is known, . Some of them are based on the internal model principle (as the approach presented in this paper). See for example, [?, ?, ?]. Other approaches are based on the use of adaptive observers for the disturbances, such as the ones presented in [?, ?]. Another approach uses the phase-locked loop structures, like in [?, ?]. In what follows we will briefly review these results.

Use of the Internal Model Principle

Feedback solutions for the rejection (attenuation) of unknown disturbances, based on the internal model principle, have been also presented in [?, ?, ?] (the model of the plant is considered known). The adaptive approach presented in these papers is similar to the one presented in this paper.

Approach I. In [?] an adaptive controller design approach based on parametrizing the set of stabilizing controllers using the Youla -Kucera parametrization and the internal model principle is implemented to solve a noise cancellation problem in an acoustic duct. The design approach is presented in [?],

The plant considered is presented in figure 4. It consists of a duct, a pair of speakers, a pair of microphones and their amplifiers, and an anti-aliasing filter. The first speaker is used to generate the noise signal, the second one is used to generate the control signal to cancel the noise at a particular location in the duct (disturbance and control speakers). The two microphones are used as measurement and performance evaluation.

The model of the plant has been obtained by system identification. The identification experiment was performed by first exciting the disturbance [START_REF] Constantinescu | Commande robuste et adaptative d'une suspension active[END_REF] One can use other signals for disturbance estimation in order to avoid this problem. One possibility is to replace the signal y(t) by ȳ(t) = y(t) -q -d B A u(t). speaker with white noise and collecting data from the output of the antialiasing filter. The same identification procedure was repeated for the control speaker. A spectrum analyzer provided the frequency response of the plant based on the collected data.

Spectrum Analyzer

Anti

The experiment considered consists of exciting the disturbance speaker with a single tone sinusoidal signal and running the adaptive controller in order to cancel the noise at the performance microphone. The Youla -Kucera filter to be adjusted online is a first order FIR filter.

The experimental results show a mixed performance for the adaptive systems with good performance only in certain disturbance frequency ranges. At some ranges of the disturbance frequency, the adaptive systems achieved 20 to 50 dB reduction in the performance signal PSD (Power Spectral Density) at the disturbance frequency and also good overall reduction in the PSD of the signal. At some other disturbance frequency ranges, the performance was not as good. It has been concluded that the unmodeled plant nonlinearities are most likely to be responsible for the observed deterioration of the adaptive system performance in experiments as compared to simulations.

Approach II. An algorithm very similar to the one presented in this paper and based on the adaptive feedback control for disturbance rejection is presented in [?].

The application considered to illustrate the performances of the proposed algorithm is the Saccharomyces cerevisiae fed-batch fermentation. It is an example of a system subject to an unstable disturbance (but over a finite time).

The objective of the application was to maximize the biomass productivity in the reactor. Hence, the substrate concentration in the reactor needs to be kept at the critical value in order to maintain a constant ethanol concentration. However, the substrate concentration is unknown a priori and may change from experiment to experiment.

So the problem of maximizing the production of biomass is converted into that of regulating the concentration of ethanol. The plant is considered as an integrator that constantly produces ethanol. The substrate consumption by the cells is modeled as an exponential disturbance which the controller has to reject. The methodology has been successfully implemented in an experimental setup in the laboratory. So the control objective is to keep the concentration of this fermentation product constant.

The control strategy is presented in figure 5. The controller C generates the substrate feed rate F that is capable of maintaining exponential cell growth and the desired production of ethanol. The controller will be able to reject the exponential substrate consumption disturbance w, if the rate at which the substrate is supplied is sufficient. Since the cells grow exponentially, it is expected that the substrate feed rate also be exponential. An adaptive disturbance rejection methodology is used to control the fed-batch fermentation near its optimal operating point. The biomass and ethanol productions rates have been decoupled, as it can be seen in figure 5, where G w relates the rate of substrate oxidation reaction to the rate at which the substrate is fed (the discrete transfer function of the disturbance) and G p relates the substrate feed rate to the ethanol production rate. S and P are the concentrations of the substrate and ethanol, respectively, and P sp is the ethanol concentration setpoint. The derived linear models have been used successfully to implement the adaptive disturbance rejection methodology. The ethanol concentration in the bioreactor is the only measurement required. The only parameter that is estimated on-line is the exponential cell growth rate. Changes in the disturbance structure are automatically compensated.

G w

Use of Adaptive Observers

The design of an output feedback adaptive observer-based compensator for the rejection of additive noisy sinusoidal disturbances with unknown bias, magnitude, phase and frequency is presented in [?]. In this approach the system is considered asymptotically stable, linear and observable, its parameters being supposed known. For a system of order n, the order of the compensator generating asymptotically convergent estimates of the biased sinusoidal disturbance is 2n + 6. The estimated parameters of the sinusoidal disturbance are its magnitude, frequency (with known lower bound), phase and bias. It is shown that the closed-loop error system based on adaptive observer is globally asymptotically stable. The adaptive observer-based compensation technique presented in [?] is illustrated on an example taken from [?]. In order to test the robustness of the control strategy, a more general disturbance has been considered, compared to the disturbance considered in [?]. Hence, bias and additional unmodeled noise have been taken into account, neither of them being considered in [?]. The performance in simulation are interesting. While this scheme does not directly make reference to the "internal model principle", one can however remark that through the observer the model of the disturbance is incorporated in the controller.

An extension of the above approach with explicit reference to the internal model principle for the case of nonlinear plants is discussed in [?].

Use of the Phase-Locked Loop Structure

A direct approach for the rejection of sinusoidal disturbances with unknown frequencies, based on the integration of a phase-locked loop for adaptive feedback control with known plant model is presented in [?]. The disturbance frequency estimation and the disturbance cancellation are performed simultaneously using a single error signal. The elimination of the high-frequency components within the system is done by using a low-pass compensator, no additional filtering being necessary. The knowledge of the frequency response of the plant in the frequency range of interest is required. Because of the lockin range of the phase-locked loops, there exist an upper limit on the initial errors of the disturbance frequencies. The system presented in [?] is able to compensate sinusoidal disturbances in a range ±30%. Experimental results are mentionned in [?]. An applicable scheme of this type would need to combine the direct algorithm with an initialization scheme for providing rough initial estimates of the disturbance frequency.

An application -Adaptive rejection of narrow band disturbances on an active suspension

The active suspension

The structure of the system (the active suspension) used in this paper is presented in figure 6. Two photos of the system are presented in figure 7 (Courtesy of Hutchinson Research Center and Laboratoire d'Automatique de Grenoble). It consists of the active suspension, a load, a shaker and the components of the control scheme. The mechanical construction of the load is such that the vibrations produced by the shaker, fixed to the ground, are transmitted to the upper side of the active suspension. The active suspension is based on a hydraulic system allowing to reduce the over-pressure at the frequencies of the vibration modes of the suspension. Its components are: an elastomer cone (1) which marks the main chamber filled up with silicon oil, a secondary chamber (2) marked by a flexible membrane, a piston (3) attached to a motor, an orifice (4) allowing the oil to pass between the two chambers and a force sensor located between the support and the active suspension.

The controller will act upon the piston (through a power amplifier) in order to reduce the residual force. The sampling frequency is 800Hz. The equivalent control scheme is shown in figure 8. The system input, u(t) is the position of the piston (see figures 6, 8), the output y(t) is the residual force measured by a force sensor. The transfer function (q -d1 C D ), between the disturbance force, u p , and the residual force y(t) is called primary path. In our case (for testing purposes), the primary force is generated by a shaker controlled by a signal

given by the computer. The plant transfer function (q -d B A

) between the input of the system, u(t), and the residual force is called secondary path. The input of the system being a position and the output a force, the secondary path transfer function has a double differentiator behavior. The control objective is to reject the effect of unknown narrow band disturbances on the output of the system (residual force), i.e. to attenuate the vibrations transmitted from the machine to the support via the active suspension. The physical parameters of the active suspension system are not provided by the manufacturer. The system has to be considered as a "black box".

Results obtained on the active suspension

The narrow band disturbance rejection procedure using the direct adaptive control scheme proposed in section 4 is illustrated in real time for the case of the control of the active suspension (presented in section 7.1). In our case the disturbance will be a time-varying frequency sinusoid, so we shall consider n Dp = 2 and n Q = n Dp -1 = 1. Furthermore a comparison with the results obtained in real time with the indirect adaptive control scheme proposed in section 5 will be provided. The identification procedure in open and closed-loop operation for the active suspension is discussed in detail in [?, ?, ?, ?, ?]. The frequency characteristic of the identified primary path model (open-loop identification), between the signal u p sent to the shaker in order to generate the disturbance and the residual force y(t), is presented in figure 9. The first vibration mode of the primary path model is near 32Hz. The frequency characteristic of the iden- tified secondary path model (closed-loop identification), is presented also in figure 9. This model has the following complexity: n B = 14, n A = 16, d = 0. The identification has been done using as excitation of the piston a PRBS (Pseudo Random Binary Sequence) with frequency divider p = 4 (for details on the PRBS signals see [?]). There exist several very low damped vibration modes on the secondary path, the first one being at 31.8Hz with a damping factor 0.07. The identified model of the secondary path has been used for the design and implementation of the controller.

The central controller (without the internal model of the disturbance) has been designed using the pole placement method and the secondary path identified model. A pair of dominant poles has been fixed at the frequency of the first vibration mode (31.8Hz), with a damping ξ = 0.8, and the other poles of the model have been considered as auxiliary desired closed loop poles. In addition a pre-specified part H R = 1 + q -1 (R = H R R ′ ) which assures the opening of the loop at 0.5f s has been introduced and 10 auxiliary poles at 0.7 have been added to the desired closed-loop poles. The resulting nominal controller has the following complexity: n R = 14, n S = 16 and it satisfies the imposed robustness constraints on the sensitivity functions(for the design methodology see[?]) [START_REF] Constantinescu | Direct controller order reduction by identification in closed loop applied to a benchmark problem[END_REF] .

In order to evaluate the performances of direct and indirect methods in real time, time-varying frequency sinusoidal disturbances between 25 and 47Hz have been used (the first vibration mode of the primary path is near 32Hz).

For both direct and indirect adaptive control methods, two protocols have been defined: one for a self-tuning operation, the other for an adaptive operation.

• Protocol 1 : Self-tuning operation

The system operates in closed loop with a frozen controller. As soon as a change of the disturbance is detected (by measuring the variance of the residual output), the estimation algorithm is started with the last frozen controller in operation. When the algorithm converges (a criterion has to be defined -see below), a new controller is computed and applied to the system. The adaptation algorithm is stopped and one waits for a change of frequency.

• Protocol 2 : Adaptive operation

The estimation algorithm works permanently (once the loop is closed) and the controller is recomputed at each sampling. The adaptation gain in this case does not tend asymptotically to zero. • Start up: For comparison purpose the system is started in open-loop for both protocols. After 5 seconds (4000 samples) a sinusoidal disturbance of 32Hz is applied on the shaker. The model of the disturbance is estimated and an initial controller is computed (same initial controller for both direct and indirect adaptive control). In the case of the self-tuning operation the adaptation algorithm is stopped while in the case of the adaptive operation the adaptation algorithm continues to be active.

After the start up ends, every 15 seconds (8000 samples) sinusoidal disturbances of different frequency are applied (32Hz, 25Hz,32Hz,47Hz,32Hz).

Protocol 1 : Self-tuning operation. Real time experimental results

The measured residual force obtained in self-tuning operation with the direct adaptation method is presented in figure 10 and with the indirect adaptation method in figure 11 . We note in general a faster convergence speed of the direct adaptive control scheme compared to the indirect one (except for 47Hz). For the self-tuning protocol, the spectral densities of the residual force obtained in open and in closed loop, respectively, using the direct adaptation scheme (after the algorithm converges) are presented in figure 12. The results are given for the three frequencies used: 25, 32 and 47 Hz. We remark that the attenuations are larger than 49 dB for all the frequencies. Similar results are obtained with the indirect adaptation algorithm. For details see [?].

We note the appearance of two harmonics of the first vibration mode of the primary path on the spectral density in open loop when the frequency of the disturbance corresponds with the first resonance mode of the system (32 Hz). They appear in open loop because of the non-linearities of the system at large signals (there is an important amplification of the disturbance at the resonance frequency of the system in open loop). The harmonics do not appear in closed loop operation.

In self-tuning operation, one uses an adaptation gain F (t) with variable forgetting factor, with λ 0 = 0.97 and the initial forgetting factor λ 1 (0) = 0.97 (the forgetting factor is given by λ 1 (t) = λ 0 λ 1 (t-1)+1-λ 0 , with 0 < λ 0 < 1). For the variable forgetting factor the adaptation gain tends asymptotically towards zero. The convergence criterion has been fixed as a threshold on the trace value of the adaptation gain matrix. For details see [?].

The detection of a change of frequency is done using the variance of the measured residual force computed on a sliding window of 50 samples.

Protocol 2 : Adaptive operation. Real time experimental results

The measured residual force obtained in adaptive operation is presented in figure 13 for the direct adaptation method and in figure 14 for the indirect adaptation method. An adaptation gain with variable forgetting factor combined with a constant trace ([?, ?]) has been used in order to be able to track automatically the changes of disturbance characteristics. The low level threshold of the trace has been fixed at 3•10 -9 for the direct algorithm and at 5•10 -7 for the indirect one (note that in the indirect adaptive scheme there are more parameters to estimate than in the direct adaptive scheme). The attenuation obtained with the indirect adaptive scheme in adaptive operation is less good than in the self tuning operation and less good than the one obtained with the direct adaptive scheme. This is certainly caused by the phenomenon discussed in section 5. We note that the direct adaptive control scheme in adaptive operation gives better results than in self tuning operation (compare figures 10 and 13).

The spectral densities of the residual force for the direct adaptive scheme (after the algorithm converges) are similar with those obtained in self-tuning operation (see [?]).

According to the real time results presented above, one can conclude that the direct adaptive control scheme gives better results than the indirect adaptive control scheme, from the point of view of the convergence speed and performance. In addition the direct adaptation scheme is much simpler than the indirect one in terms of number of operations.

Direct adaptive control scheme under the effect of sinusoidal disturbances with continuously time varying frequency

Consider now that the frequency of the sinusoidal disturbance varies continuously and let's use a chirp disturbance signal (linear swept-frequency signal) between 25 and 47Hz. The tests have been done as follows: Start up in closed loop at t = 0 with the central controller. Once the loop is closed, the adaptation algorithm works permanently and the controller is updated (direct approach) at each sampling instant. After 5 seconds a sinusoidal disturbance of 25 Hz (constant frequency) is applied on the shaker. From 10 to 15 seconds a chirp between 25 and 47 Hz is applied. After 15 seconds a 47 Hz (constant frequency) sinusoidal disturbance is applied and the tests are stopped after 18 seconds. The time-domain results obtained in open and in closed-loop (direct adaptive control) are presented in figure 15. We can remark that the performances obtained are very good. Adaptation transients for direct adaptive control

Conclusions

It was shown in this paper that the use of the internal model principle combined with the adaptation of the internal model implemented in a Youla -Kucera parametrized controller allows a very good rejection of the unknown time varying narrow band disturbances without requiring the use of an additional transducer. Two adaptive approaches (direct and indirect adaptation) have been presented and tested comparatively on an active suspension. The results obtained in real time lead us to conclude that the direct adaptive control scheme provides better performance than the indirect adaptive control scheme. Furthermore, from the performances point of view, the adaptive operation is more interesting than the self-tuning one for the direct adaptive control scheme. Moreover, the direct algorithm is much simpler than the indirect one.

A similar approach has been used successfully on a chemical reactor and for noise cancellation in ducts. One can conclude that an efficient technique for rejection of unknown finite band disturbances is available.

For future research it will be interesting to evaluate comparatively this approach with the other approaches indicated in section 6. Establishing connections with "repetitive control" may be also very rewarding.

Fig. 1 .

 1 Fig.1. "Signal processing" approach to rejection of unknown disturbances

Fig. 2 .Fig. 3 .

 23 Fig. 2. Indirect adaptive control scheme for rejection of unknown disturbances

Fig. 4 .

 4 Fig. 4. Block diagram of the noise cancellation system

Fig. 5 .

 5 Fig. 5. Control strategy of the fed-batch fermenter

Fig. 6 .

 6 Fig. 6. Active suspension system (scheme)

Fig. 7 . 1 Fig. 8 .

 718 Fig. 7. Active suspension system (photo)

Fig. 9 .

 9 Fig. 9. Frequency characteristics of the primary and secondary paths

Fig. 10 .

 10 Fig. 10. Time domain results with the direct adaptation method in self-tuning operation

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Time domain results with the indirect adaptation method in self-tuning operation

Fig. 13 .

 13 Fig. 13. Time domain results with the direct adaptation method in the adaptive case (trace = 3 • 10 -9 )

Figure 16 illustrates

 16 Figure16illustrates the adaptation transients on the input and output when a step change of the frequency of the disturbance occurs from 20Hz to 32 Hz respectively. One notes that the convergence of the output requires less than 0.25s This corresponds roughly to 6 periods for 32Hz. Same duration of the adaptation transient are obtained for the other frequencies step changes. These results have to be compared with the transients results given in [?, ?, ?, ?].

Fig. 14 .

 14 Fig. 14. Time domain results with the indirect adaptation method in the adaptive case (trace = 5 • 10 -7 )

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Real-time results obtained with the direct adaptive method and a chirp disturbance: (a) Open loop; (b) Closed loop

Throughout the paper it is assumed that the order of the disturbance model is known but the parameters of the model are unknown (the order can be estimated from data if necessary).

Note that the adaptive rejection of unknown disturbances using the Youla-Kucera parametrization is not considered in the survey [?].

A brief review of the approaches presented in [?] and [?, ?] can be found in section

The complex variable z -1 will be used for characterizing the system's behavior in the frequency domain and the delay operator q -1 will be used for describing the system's behavior in the time domain.

Of course it is assumed that Dp and B do not have common factors.

The magnitude of the adaptation gain gives an indication upon the variance of the parameter estimation error -see for example [?].