B Gatmiri

S Le Pense

P Maghoul

A MULTI-SCALE SEISMIC RESPONSE OF TWO-DIMENSIONAL SEDIMENTARY VALLEYS DUE TO THE COMBINED EFFECTS OF TOPOGRAPHY AND GEOLOGY

It is well-known that the response of a site to seismic excitation depends on the local topographical and geological conditions. The current building codes already take into account unidimensional site effects but ignore complex site effects due to two-dimensional irregular configurations. The aim of this work is to contribute to the establishment of a simple predictive method to estimate site effects. The horizontal ground movements at the surface of sedimentary valleys subjected to SV waves with vertical incidence are calculated by using the HYBRID program, combining finite elements in the near field and boundary elements in the far field (FEM/BEM). A parametric study is conducted to examine the combined effects of topography and geology on the amplification of the response spectrum at various points across the valley. The influence of different parameters is considered, such as filling ratios (from empty to full valleys), impedance contrasts between bedrock and sediments, and dimensions.

Keywords: Seismic site effect; Seismic amplification; Building codes; Seismic response spectrum; Hybrid numerical method; Sedimentary valley; 2D configurations

Introduction

It has been often observed that, after earthquakes, some buildings located on hills or in sedimentary basins suffer greater damage than might have been expected from their distance to the epicentre. The modification of the seismic movement due to local topographical and geotechnical conditions is called site effect.

Even though one-dimensional (1D) site effects have already been incorporated in current building codes, there is still progress to be made in taking into consideration complex site effects due to two-dimensional (2D) irregular configurations.

Much work has been done since the 1980s to improve the understanding of the physical phenomena involved in site effects and their prediction. Analytical solutions have been developed for simple configurations [START_REF] Trifunac | Surface motion of a semi-cylindrical alluvial valley for incident plane SH waves[END_REF][START_REF] Trifunac | D. Trifunac, Scattering of plane SH waves by a semi-elliptical canyon[END_REF]7,[START_REF] Wong | Surface motion of a semi-elliptical alluvial valley for incident plane SH waves[END_REF][START_REF] Sanchez-Sesma | Diffraction of elastic SH waves by wedges[END_REF] . Simultaneously, experimental methods have been developed to make the most of in-situ measurements [START_REF] Nakamura | A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface[END_REF] .

More recently, the development of numerical methods has made it possible to study site effects for more complex configurations. Gatmiri and his coworkers have developed a program combining finite elements in the near field and boundary elements in the far field (FEM/BEM) called HYBRID [START_REF] Gatmiri | Two-dimensional transient wave propagation in anelastic saturated porous media by a Hybrid FE/BE method[END_REF][START_REF] Gatmiri | Applying a new fast numerical method to elasto-dynamic transient kernels in HYBRID wave propagation analysis[END_REF] .

Gatmiri and Arson (2008) 1 proposed a quantitative method to predict horizontal surface displacements. They used independent factors representing the contribution of different geometrical parameters controlling site effects.

Gatmiri et al. (2009) [START_REF] Gatmiri | Site-specific spectral response of seismic movement due to geometrical and geotechnical characteristics of sites[END_REF] showed that in the central zone of a sedimentary valley, geological effects prevail upon topographical effects, whereas from the mid-slope point up to the edge of the valley topographical effects prevail in the site response. They also introduced a criterion for empty valleys, S/A (Surface/Angle), with which the amplification at the top corner increases (Fig. 1).

Based on these previous studies, the aim of this work is to propose a simple criterion combining soil properties and geometrical characteristics, for estimating the amplification of earthquake response spectra in sedimentary valleys.

We computed horizontal displacements at the surface of sedimentary valleys by using the HYBRID code. A parametric study is conducted to evaluate the effects of various parameters on the amplification of the response spectrum. Calculations are made for several points at the surface of valleys (more than 20 points per valley).

The seismic solicitation is a vertically incident SV Ricker wave (Fig. 2).

In the following, soils are assumed to be dry and linearly elastic.

The predominant frequency of the incident signal is fixed and equal to 2 Hz (Fig. 3). The reference is taken as the scattering of the Ricker incident wave on a flat bedrock surface.

Method

A very common representation of seismic movement used by engineers is the response spectrum. This is a plot of the maximum response of a family of damped single-degree-of-freedom oscillators to earthquake ground motion. The pseudoacceleration response spectrum gives an estimation of the maximum acceleration at the base of the building as a function of the natural period of the oscillator T n.

For each case, we calculate the earthquake response spectrum with 5% damping and compare it with the reference spectrum (Fig. 4). The spectral ratio, SR, is defined as the maximum ratio between the pseudo-acceleration SA and the pseudoacceleration of the reference spectrum SA R , for T n between 0.1 to 1 second. The site period, T s , is defined as the natural period for which this ratio SA/SA R reaches its peak.

1D site effects

Soil properties

The mechanical parameters for sediments and bedrock are given in Table 1. The impedance contrast between sediments and bedrock characterises the soil properties "Eq. ( 1)".

β = ρ s C s ρ r C r (1)
Table 1. Mechanical parameters of materials

E (MPa) K (MPa) ν G (MPa) ρ (T/m 3 ) C (m/

Analytical solution

In the first place, we consider the simple case of an horizontal soil layer over a semi-infinite elastic space submitted to an harmonic plane S wave. We know the analytical formula of the ratio between the displacements at the surface point A and at point B located at the interface between the soil layer and bedrock "Eq. ( 2)" (Fig. 5).

u A u B = 1 cos 2 (k s • H) + β 2 • sin 2 (k s • H) (2)
• H is the soil layer height.

• k s = 2πf c C s
is the wavenumber where C s is the shear wave velocity in sediments and f c is the predominant frequency of the input signal. • β is the impedance contrast between sediments and bedrock.

The maximum amplification takes place for specific frequencies. "Eq. (3, 4)"

f n = (2n + 1) C s 4H (3) u A u B max = 1 α (4)

Site period T s

We study the seismic response of horizontal soil layers with various heights and impedance contrasts. The height H varies from 10m to 100m and the impedance contrast β varies from 0.2 to 0.5. Concerning the site period, we notice that, when H increases for a fixed β, T s first increases. If H keeps increasing, T s goes back to lower periods (Fig. 6).

Similarly, when β decreases for a fixed H, T s increases (Fig. 7). If we suppose that the maximum amplification of the earthquake spectrum takes place for the natural period corresponding to the fundamental frequency of the soil layer, we obtain an analytical value of T s as a function of H and β. "Eq. ( 5)"

T s = 1 f 0 = 4H C s = 4ρ s ρ r C r H β = 0.0027 H β (5)
For cases for which the amplification is significant, Fig. 8 shows the evolution of T s with H/β in comparison with the theoretical model T s = 0.0027H/β

For values of H/β from 65 to 300, our set of data is very close to the theoretical model. These values of H/β correspond to soil layers with fundamental frequencies between 1.2 and 5.8 Hz, which are the frequencies included into the incident signal (Fig. 3).

For higher frequencies (lower values of H/β) which are absent from the incident signal, the amplification is very weak. It does not allow us to determine an appropriate site period T s .

Lower frequencies (higher values of H/β) correspond to the cases for which the site period T s goes back to low period, for high values of H and low values of β. These frequencies being rare in the incident signal, this phenomenon may be due to a resonance of the soil layer on a second mode. From now on we will ignore these cases.

As can be seen in Fig. 8, the correlation between our data and the theory is good enough to validate our model.

Spectral ratio

For a constant site period T s the amplitude of response spectra decreases when β increases (Fig. 9) and when H increases (since these two parameters are linked by Eq. ( 5)).

We want to represent on the same graph the spectral ratios for the different soil layers, in order to give an estimation of the spectral ratio knowing H and β.

We choose to represent the parameter H/β which controls the site period T s on the x-axis, and the parameter (SR -1)H on the y-axis (Fig. 10).

The parameter (SR -1) allows us to visualise directly whether or not there is an amplification. Indeed, when there is no amplification SR -1 = 0. This parameter is then multiplied by H to represent the fact that the spectral ratio decreases when H increases for a given period.

2D site effects

Geometrical characteristics

The studied valleys are characterised by their half-width at the surface L, their half-width at the base L 1 , their depth H and the sediment layer height H 1 . Let A be the angle between the slope and the horizontal and S 1 the surface of the section filled with sediments (Fig. 11).

The half-width L is equal to 100 m. We study triangular valleys corresponding to L 1 /L = 0 and trapezoidal valleys corresponding to L 1 /L = 0.4 . The calculations are made for different shape ratios H/L = 0.2, 0.4, 0.6, 1 and filling ratios H 1 /H = 0, 0.25, 0.5, 0.75, 1.

Spatial evolution

First, we study the evolution of the spectral ratio with the spatial location of the studied point. The spectral ratio is represented as a function of the non-dimensional offset variable x/L for different configurations.

By increasing the filling ratio, we go from empty valley behaviour with a maximum amplification at the edge of the valley to fully-filled valley behaviour with maximum amplification at the centre of the valley. For a quarter-filled valley the behaviour is the same as that of an empty valley (Fig. 12).

For intermediate filling ratios, there is a local maximum at the edge of the valley, a local maximum at the centre of the valley and a minimum at the contact point between sediments and bedrock. There is also a decrease of the spectral ratio by moving away from the valley (Fig. 13).

From now on, we will study those three specific points :

• Edge of the valley (x/L = 1)

• Contact point sediments/bedrock • Centre of the valley (x/L = 0)

Topographical effect, edge of the valley

For empty valleys, the maximum amplification takes place at the edge of the valley. It increases with the shape ratio. Trapezoidal valleys entail higher amplifications than triangular valleys (Fig. 14 and15). Gatmiri et al. (2009) 2 proposed a criteria S/A with which the spectral ratio at the edge increases (Fig. 1). Our numerical results follow this criteria (Fig. 16). We also notice that the spectral ratio at the edge is higher for empty valleys compared to partially-filled ones (Fig. 12). This is verified for every dimension and geometry of valleys and soil properties. At the contact point between sediments and bedrock of partially-filled valleys, the spectral ratio is always lower than that of the corresponding empty valley with the same dimensions. It means that from this point up to the edge of the valley, topographical effects prevail upon geological effects. So, for points which are located on the rock (i.e. outside the area filled with sediments), it is possible to consider only the topographical site effect due to an empty valley. The response spectra calculated by this method will always be higher than the ones of the real cases.

We verify that the response spectra at this specific point of deep or shallow, fully or quarterly-filled valleys are lower than the ones of equivalent empty valleys (Fig. 17 For the different configurations, all the acceleration response spectra calculated at the central point of the valleys are represented on the same graph (Fig. 21).

If there is an attenuation of some spectra for low periods due to topographical effect, it is not significant compared to the amplification of other spectra. So, we will not study the attenuation cases afterwards. Following the same method as for 1D soil layers, we want to define a criterion allowing an estimation of the site period for configurations with a significant amplification.

We propose to study the evolution of the site period with the parameter S 1 /β √ β which combines the soil properties (β) and the geometrical characteristics (S 1 ) (Fig. 22).

An example of values taken by S 1 /β √ β are given in Table 2. We study the case for which 13 000 < S 1 /β √ β < 120 000 For S 1 /β √ β < 13 000 the amplification is not significant enough to determine a representative site period T s .

For S 1 /β √ β > 120 000, we ignore the case corresponding to H/L = 1, β = 0.2, H 1 /H = 1 and β √ β = 156 000 for which the site period goes back to low periods. We observe a linear evolution of T s with S 1 /β √ β (Fig. 22). As for 1D soil layers, we observe that, for a same period, spectra have a higher amplitude for soft soils ( low β) (Fig. 23 and24). By analogy with the 1D case we propose a representation of (SR -1)S 1 as a function of the parameter S 1 /β √ β governing the site period. For triangular and trapezoidal valleys our numerical results are represented on Fig. 25. We eliminated one case (H/L = 1, β = 0.2, H 1 /H = 1, S 1 /β √ β = 156 000) for which the site period goes back to low periods.

The amplification is negligible up to a threshold corresponding to S 1 /β √ β = 13 000. Afterwards, the evolution of (SR -1) • S 1 with S 1 /β √ β shows a parabolic tendency.

From these curves, knowing the geometrical characteristics of the valley S 1 and its soil properties β,we can now calculate the coefficient S 1 /β √ β. With this coefficient we can read an estimation of the site period T s and the coefficient (SR -1)•S 1 and then deduce the value of the spectral ratio SR ( Fig. 22 and Fig. 25).

Conclusions

Although 2D complex site effects are not yet taken into consideration in building codes, many studies have pointed out the importance of their impact on seismic movement. Indeed, this site effects often leads to a higher amplification of the seismic signal compared to the 1D case. It also causes an extension of the signal length.

This study gives a method to characterise the spectral amplification of bidimensional triangular and trapezoidal valleys with various filling ratios. We considered valleys with shape ratios H/L lower than 1 and impedance contrast higher than 0.2.

We first found again known results :

• For 1D soil layers we found a good correlation with the analytical solution. Indeed, the maximum amplification of the acceleration response spectrum occurs for the natural period corresponding to the fundamental frequency of the soil layer. • For empty valleys, the maximum amplification occurs at the edge of the valley and increases with the criteria S/A. • In fully-filled valleys, the maximum amplification takes place at the centre of the valley. • When increasing the filling ratio, we go from empty valley behaviour to fully-filled valley behaviour.

Fig. 1 .

 1 Fig. 1. Definition of parameter S/A

Fig. 2 . 3 .Fig. 4 .

 234 Fig. 2. Incident Ricker signal Fig. 3. Fourier transform of the reference displacement

Fig. 5 .

 5 Fig. 5. 1D site effect of a soil layer on a bedrock

Fig. 6 .

 6 Fig. 6. Effect of height on response spectra (β = 0.2)

Fig. 7 .Fig. 8 .

 78 Fig. 7. Effect of impedance contrast on response spectra (H = 40 m)

FABCDFig. 9 .

 9 Fig. 9. Acceleration response spectra of 1D soil layers (Black: reference; blue: β = 0.2; yellow: β = 0.3; red: β = 0.4

Fig. 10 .

 10 Fig. 10. Evolution of (SR -1)H with H/β

Fig. 12 .

 12 Fig. 12. Spatial evolution of spectral ratio for various filling ratios (trapezoidal valley with H/L = 0.6 and β = 0.3)

Fig. 13 .Fig. 16 . 4 . 4 .

 131644 Fig. 13. Response spectrum evolution along a trapezoidal valley (H/L = 0.6, H 1 /H = 0.5, β = 0.3) (Reference spectrum in black)

  , 18, 19 and 20 ) 4.5. Centre of the valley 4.5.1. Envelope of response spectra

Fig. 17 .Fig. 18 .Fig. 19 .Fig. 20 . 1 c

 171819201 Fig. 17. H 1 /H = 0.25, H/L = 0.2, x/L = 0.55

FABCDFig. 21 .

 21 Fig. 21. Envelope of response spectra at the centre of the valley for various geometrical characteristics and soil properties. Black line is the reference spectrum

Fig. 22 .FABCDFig. 23 .FABCDFig. 24 .

 222324 Fig. 22. Evolution of the site period with S 1 /β √ β for trapezoidal and triangular valleys

Fig. 25 .

 25 Fig. 25. Evolution of (SR -1) • S 1 with S 1 /β √ β

Table 2 .

 2 Values taken by S 1 /β √ β for trapezoidal valleys with different shape ratios and filling ratios, β = 0.3

& 8 Avenue Blaise Pascal,

Marne la Vallée, France. solenn.le-pense@enpc.fr

B. GATMIRI, S. LE PENSE AND P. MAGHOUL

B. GATMIRI, S. LE PENSE AND P. MAGHOUL

B. GATMIRI, S. LE PENSE AND P. MAGHOUL

SEISMIC RESPONSE OF TWO-DIMENSIONAL SEDIMENTARY VALLEYS 15

The main results of this study are :

• In an alluvial valley, fully-or partially-filled with sediments, from the contact point between sediments and bedrock up to the edge and outside of valleys, topographical effects prevail upon geological effects. • At the central point of the valley, the evolution of the site period T s with the parameter S 1 /β √ β has a linear tendency. This parameter combines the soil properties and the geometrical characteristics of the valley.

• At the centre of the valley we can estimate the spectral ratio SR from the curve representing the evolution of (SR -1)S 1 as a function of S 1 /β √ β (Fig. 25).

Appendix A. Appendices


Pedro Casas
email: pcasashe@laas.fr

Johan Mazel
email: jmazel@laas.fr

Philippe Owezarski

Steps Towards Autonomous Network Security: Unsupervised Detection of Network Attacks

Keywords: Unsupervised Anomaly Detection & Characterization, Robust & Distributed Clustering, Automatic Generation of Signatures, Autonomous Security

The unsupervised detection of network attacks represents an extremely challenging goal. Current methods rely on either very specialized signatures of previously seen attacks, or on expensive and difficult to produce labeled traffic datasets for profiling and training. In this paper we present a completely unsupervised approach to detect attacks, without relying on signatures, labeled traffic, or training. The method uses robust clustering techniques to detect anomalous traffic flows, sequentially captured in a temporal sliding-window basis. The structure of the anomaly identified by the clustering algorithms is used to automatically construct specific filtering rules that characterize its nature, providing easy-to-interpret information to the network operator. In addition, these rules are combined to create an anomaly signature, which can be directly exported towards standard security devices like IDSs, IPSs, and/or Firewalls. The clustering algorithms are highly adapted for parallel computation, which permits to perform the unsupervised detection and construction of signatures in an online basis. We evaluate the performance of this new approach to discover and to build signatures for different network attacks without any previous knowledge, using real traffic traces. Results show that knowledge-independent detection and characterization of network attacks is possible, opening the door to a whole new generation of autonomous security algorithms.

I. INTRODUCTION

The detection of network attacks is a paramount task for network operators in today's Internet. Denial of Service attacks (DoS), Distributed DoS (DDoS), network/host scans, and spreading worms or viruses are examples of the different attacks that daily threaten the integrity and normal operation of the network. The principal challenge in automatically detecting and analyzing network attacks is that these are a moving and ever-growing target [START_REF] Hansman | A Taxonomy of Network and Computer Attacks[END_REF].

Two different approaches are by far dominant in the literature and commercial security devices: signature-based detection and anomaly detection. Signature-based detection systems are highly effective to detect those attacks which they are programmed to alert on. However, they cannot defend the network against unknown attacks. Even more, building new signatures is expensive and time-consuming, as it involves manual inspection by human experts. Anomaly detection uses labeled data to build normal-operation-traffic profiles, detecting anomalies as activities that deviate from this baseline. Such methods can detect new kinds of network attacks not seen before. Nevertheless, anomaly detection requires training to construct normal-operation profiles, which is time-consuming and depends on the availability of purely anomaly-free traffic data-sets. In addition, it is not easy to maintain an accurate and up-to-date normal-operation profile.

In this paper we present a completely unsupervised method to detect and characterize network attacks, without relying on signatures, training, or labeled traffic of any kind. Our approach relies on robust clustering algorithms to detect both well-known as well as completely unknown attacks, and to automatically produce easy-to-interpret signatures to characterize them, both in an on-line basis. The analysis is performed on packet-level traffic, captured in consecutive time slots of fixed length ΔT and aggregated in IP flows (standard 5-tuples). IP flows are additionally aggregated at 9 different flow levels l i . These include (from finer to coarsergrained resolution): source IPs (l 1 : IPsrc), destination IPs (l 2 : IPdst), source Network Prefixes (l 3,4,5 : IPsrc/24, /16, /8), destination Network Prefixes (l 6,7,8 : IPdst/24, /16, /8), and traffic per Time Slot (l 9 : tpTS).

The complete detection and characterization algorithm runs in three successive stages. The first step consists in detecting an anomalous time slot where an attack might be hidden. For doing so, time series Z li t are built for basic traffic metrics such as number of bytes, packets, and IP flows per time slot, using the 9 flow resolutions l 1...9 . Any generic anomalydetection algorithm F(.) based on time-series analysis [START_REF] Barford | A Signal Analysis of Network Traffic Anomalies[END_REF]- [START_REF] Cormode | What's New: Finding Significant Differences in Network Data Streams[END_REF] is then used on Z li t to identify an anomalous slot. Time slot t 0 is flagged as anomalous if F(Z li t0 ) triggers an alarm for any of the l i flow aggregation levels. Tracking anomalies at multiple aggregation levels provides additional reliability to the anomaly detector, and permits to detect both single sourcedestination and distributed attacks of very different intensities.

The unsupervised detection and characterization algorithm begins in the second stage, using as input the set of IP flows captured in the flagged time slot. The method uses robust clustering techniques based on Sub-Space Clustering (SSC) [START_REF] Parsons | Subspace Clustering for High Dimensional Data: a Review[END_REF], Density-based Clustering [START_REF] Ester | A Density-based Algorithm for Discovering Clusters in Large Spatial Databases with Noise[END_REF], and Evidence Accumulation (EA) [START_REF] Fred | Combining Multiple Clusterings Using Evidence Accumulation[END_REF] to blindly extract the suspicious flows that compose the attack. In the third stage, the evidence of traffic structure provided by the clustering algorithms is used to produce filtering rules that characterize the detected attack and simplify its analysis. The characterization of an attack can be a hard and time-consuming task, particularly when dealing with unknown attacks. Even expert operators can be quickly overwhelmed if simple and easy-to-interpret information is not provided to prioritize the time spent in the analysis. To alleviate this issue, the most relevant filtering rules are combined into a new traffic signature that characterizes the attack in simple terms. This signature can ultimately be integrated to any standard security device to detect the attack in the future, which constitutes a major step towards autonomous security: in a nutshell, our algorithm automatically produces new signatures without any previous data about traffic or knowledge about the attack.

The remainder of the paper is organized as follows. Section II presents a short state of the art in the unsupervised anomaly detection field and describes our main contributions. Section III briefly describes the unsupervised detection algorithm that we have developed. Section IV presents the automatic characterization algorithm, which builds easy-tointerpret signatures for the detected attacks. Section V presents the validation of our proposals, discovering and characterizing single source/destination and distributed network attacks in traffic traces from an operational backbone network. Section VI evaluates the computational time of the unsupervised detection approach, considering the parallelization of the clustering algorithms. Finally, section VII concludes the paper.

II. RELATED WORK & CONTRIBUTIONS

The problem of network attacks and anomaly detection has been extensively studied in the last decade. Most approaches analyze statistical variations of traffic volume-metrics (e.g., number of bytes, packets, or flows) and/or other traffic features (e.g. distribution of IP addresses and ports), using either singlelink measurements or network-wide data. A non-exhaustive list of methods includes the use of signal processing techniques (e.g., ARIMA, wavelets) on single-link traffic measurements [START_REF] Barford | A Signal Analysis of Network Traffic Anomalies[END_REF], [START_REF] Brutlag | Aberrant Behavior Detection in Time Series for Network Monitoring[END_REF], PCA [START_REF] Lakhina | Diagnosing Network-Wide Traffic Anomalies[END_REF], [START_REF] Lakhina | Mining Anomalies Using Traffic Feature Distributions[END_REF] and Kalman filters [START_REF] Soule | Combining Filtering and Statistical Methods for Anomaly Detection[END_REF] for network-wide anomaly detection, and sketches applied to IP-flows [START_REF] Krishnamurthy | Sketch-based Change Detection: Methods, Evaluation, and Applications[END_REF], [START_REF] Dewaele | Extracting Hidden Anomalies using Sketch and non Gaussian Multi-resolution Statistical Detection Procedures[END_REF].

Our approach falls within the unsupervised anomaly detection domain. Most work has been devoted to the Intrusion Detection field, targeting the well known KDD'99 data-set. The vast majority of the unsupervised detection schemes proposed in the literature are based on clustering and outliers detection, being [START_REF] Portnoy | Intrusion Detection with Unlabeled Data Using Clustering[END_REF]- [START_REF] Leung | Unsupervised Anomaly Detection in Network Intrusion Detection Using Clustering[END_REF] some relevant examples. In [START_REF] Portnoy | Intrusion Detection with Unlabeled Data Using Clustering[END_REF], authors use a single-linkage hierarchical clustering method to cluster data from the KDD'99 data-set, based on the standard Euclidean distance for inter-patterns similarity. [START_REF] Eskin | A Geometric Framework for Unsupervised Anomaly Detection: Detecting Intrusions in Unlabeled Data[END_REF] reports improved results in the same data-set, using three different clustering algorithms: Fixed-Width clustering, an optimized version of k-NN, and one class SVM. [START_REF] Leung | Unsupervised Anomaly Detection in Network Intrusion Detection Using Clustering[END_REF] presents a combined density-grid-based clustering algorithm to improve computational complexity, obtaining similar detection results.

Our unsupervised algorithm has several advantages w.r.t. the state of the art: (i) first and most important, it works in a completely unsupervised fashion, which means that it can be directly plugged-in to any monitoring system and start to work from scratch, without any kind of calibration or previous knowledge. (ii) It combines robust clustering techniques to avoid general clustering problems such as sensitivity to initialization, specification of number of clusters, or structure-masking by irrelevant features. (iii) It automatically builds compact and easy-to-interpret signatures to characterize attacks, which can be directly integrated into any traditional security device. (iv) It is designed to work on-line, using the parallel structure of the proposed clustering approach.

III. UNSUPERVISED DETECTION OF ATTACKS

The unsupervised detection stage takes as input all the IP flows in the anomalous time slot, aggregated according to one of the different aggregation levels used in the first stage. Let Y = {y 1 , .., y n } be the set of n flows in the flagged time slot. Each flow y i ∈ Y is described by a set of m traffic attributes or features on which the analysis is performed. The selection of these features is a key issue to any anomaly detection algorithm, and it becomes critical in the case of unsupervised detection, because there is no additional information to select the most relevant set. In this paper we shall limit our study to detect and characterize well-known attacks, using a set of standard traffic features widely used in the literature. However, the reader should note that the approach can be easily extended to detect other types of attacks, considering different sets of traffic features. In fact, more features can be added to any standard list to improve detection and characterization results.

The set that we shall use here includes the following m = 9 traffic features: number of source/destination IP addresses and ports, ratio of number of sources to number of destinations, packet rate, ratio of packets to number of destinations, and fraction of ICMP and SYN packets. According to previous work on signature-based anomaly characterization [START_REF] Fernandes | Automated Classification of Network Traffic Anomalies[END_REF], such simple traffic descriptors permit to describe standard network attacks such as DoS, DDoS, scans, and spreading worms/virus. Let x i = (x i (1), .., x i (m)) ∈ R m be the corresponding vector of traffic features describing flow y i , and X = {x 1 , .., x n } the complete matrix of features, refereed to as the feature space.

The algorithm is based on clustering techniques applied to X. The objective of clustering is to partition a set of unlabeled elements into homogeneous groups of similar characteristics, based on some measure of similarity. Our goal is to identify in Y the different aggregated flows that may compose the attack. For doing so, the reader should note that an attack may consist of either outliers (i.e., single isolated flows) or compact smallsize clusters, depending on the aggregation level of flows in Y. For example, a DDoS attack is represented as an outlier flow if the aggregation is done for IPdst, consisting of all the attacking IP flows sent towards the same victim. On the contrary, the attack is represented as a cluster if we use IPsrc flow-resolution. To avoid the lack of robustness of general clustering techniques, we have developed a parallel-multiclustering approach, combining the notions of Density-based Clustering [START_REF] Ester | A Density-based Algorithm for Discovering Clusters in Large Spatial Databases with Noise[END_REF], Sub-Space Clustering [START_REF] Parsons | Subspace Clustering for High Dimensional Data: a Review[END_REF], and Evidence Accumulation [START_REF] Fred | Combining Multiple Clusterings Using Evidence Accumulation[END_REF]. The particular details of the algorithm are fully documented in [START_REF] Casas | Sub-Space Clustering & Evidence Accumulation for Unsupervised Network Anomaly Detection[END_REF]; in what follows, we shall present the general idea behind the approach.

Instead of directly partitioning the complete feature space X using a traditional inter-flow similarity measure (i.e., the Euclidean distance), we do parallel clustering in N different sub-spaces X i ⊂ X of smaller dimensions, obtaining N different partitions P i of the flows in Y. Each sub-space X i is constructed using only r < m traffic features; this permits to analyze the structure of X from N (m, r) different perspectives, using a finer-grained resolution. In particular, we do clustering in very-low dimensional sub-spaces, using r = 2. To deeply explore the complete feature space, we analyze all the r-combinations-obtained-from-m sub-spaces; hence, N (m) = m(m -1)/2. The information provided by the multiple partitions P i is then combined to produce a new similarity measure between the flows in Y, which has the paramount advantage of clearly highlighting both those outliers and small-size clusters that were simultaneously identified in different sub-spaces. This new similarity measure is finally used to easily extract the anomalous flows from the rest of the traffic. Briefly speaking, if we can find a group of flows that are remarkably different from the rest of the traffic in different sub-spaces, then we have found an anomaly; if not, the flagged time-slot was just a false alarm.

IV. AUTOMATIC CHARACTERIZATION OF ATTACKS

The following task after the detection of a group of anomalous flows is to automatically produce a set of K filtering rules f k (Y), k = 1, .., K to characterize them. In the one hand, such filtering rules provide useful insights on the nature of the anomaly, easing the analysis task of the network operator. On the other hand, different rules can be combined to construct a signature of the anomaly, which can be used to easily detect its occurrence in the future. To produce filtering rules f k (Y), the algorithm selects those sub-spaces X i where the separation between the anomalous flows and the rest of the traffic is the biggest. We define two different classes of filtering rule: absolute rules f A (Y) and relative rules f R (Y). Absolute rules are only used in the characterization of small-size clusters, and correspond to the presence of dominant features in the flows of the anomalous cluster. An absolute rule for feature j has the form f A (Y) = {y i ∈ Y : x i (j) == λ}. For example, in the case of an ICMP flooding attack, the vast majority of the associated flows use only ICMP packets, hence the absolute filtering rule {nICMP/nPkts == 1} makes sense (nICMP/nPkts corresponds to the fraction of ICMP packets).

On the other hand, relative filtering rules depend on the relative separation between anomalous and normal-operation flows. Basically, if the anomalous flows are well separated from the rest of the traffic in a certain partition P i , then the features of the corresponding sub-space X i are good candidates to define a relative filtering rule. A relative rule defined for feature j has the form f R (Y) = {y i ∈ Y : x i (j) < λ or x i (j) > λ}. We shall also define a covering relation between filtering rules: we say that rule

f 1 covers rule f 2 ↔ f 2 (Y) ⊂ f 1 (Y).
If two or more rules overlap (i.e., they are associated to the same feature), the algorithm keeps the one that covers the rest.

In order to construct a compact signature of the anomaly, we have to devise a procedure to select the most discriminant filtering rules. Absolute rules are important, because they define inherent characteristics of the anomaly. Regarding relatives rules, their relevance is directly tied to the degree of separation between flows. In the case of outliers, we select the K features for which the normalized distance to the normaloperation traffic (statistically represented by the biggest cluster in each sub-space) is among the top-K biggest distances. In the case of small-size clusters, we rank the degree of separation to the rest of the clusters using the well-known Fisher Score (FS) [START_REF] Jaakkola | Exploiting Generative Models in Discriminative Classifiers[END_REF], and select the top-K ranked rules. The FS basically measures the separation between clusters, relative to the total variance within each cluster. To finally construct the signature, the absolute rules and the top-K relative rules are combined into a single inclusive predicate, using the covering relation in case of overlapping rules.

V. EXPERIMENTAL EVALUATION

We evaluate the ability of the unsupervised algorithm to detect and to automatically construct a signature for different attacks in real traffic from the WIDE project data repository [START_REF] Cho | Traffic Data Repository at the WIDE Project[END_REF]. The WIDE network provides interconnection between different research institutions in Japan, as well as connection to different commercial ISPs and universities in the U.S.. Traffic consists of 15 minutes-long raw packet traces; the traces we shall work with consist of packets captured at one of the trans-pacific links between Japan and the U.S.. Traces are not labeled, thus our analysis will be limited to show how the unsupervised approach can detect and characterize different network attacks without using signatures, labels, or learning.

We shall begin by detecting and characterizing a distributed SYN network scan directed to many victim hosts under the same /16 destination network. Packets in Y are aggregated using IPdst/24 flow resolution, thus the attack is detected as a small-size cluster. The length of each time slot is ΔT = 20 seconds. As we explained in section III, the SSC-EAbased clustering algorithm constructs a new similarity measure between flows in Y, using the multiple clustering results obtained from the different sub-spaces. Let us express this new similarity measure as a n × n matrix S, in which element S(i, j) represents the degree of similarity between flows i and j. Figure 1.(a) depicts a histogram on the distribution of inter-flows similarity, according to S. The structure of flows in Y provided by S evidences the presence of a small isolated cluster in multiple sub-spaces. Selecting this cluster results in 53 anomalous IPdst/24 flows; a further analysis of the packets in these flows reveals multiple IP flows of SYN packets with the same IPsrc address and sequential IPdst addresses, scanning primary the same TCP port. Such a behavior is characteristic of a worm in the spreading phase.

Regarding filtering rules, figures 1.(b,c) depict some of the partitions P i where both absolute and top-K relative rules were produced. These involve the number of sources and destinations, and the fraction of SYN packets. Combining them produces a signature that can be expressed as (nSrcs == 1) ∧ (nDsts > λ 1 ) ∧ (nSYN/nPkts > λ 2 ), where both λ 1 and λ 2 are obtained by separating clusters at half distance. Surprisingly enough, the extracted signature matches quite closely the standard signature used to detect such an attack in current signature-based systems [START_REF] Fernandes | Automated Classification of Network Traffic Anomalies[END_REF]. The beauty and main advantage of our unsupervised approach relies on the fact that this new signature has been produced without any previous information about the attack or baseline traffic, and now it can be directly exported towards any security device to rapidly detect the same attack in the future. Figures 1.(d,e) depict different rules obtained in the detection of a SYN DDoS attack. IP flows are now aggregated according to IPsrc resolution. The distribution analysis of inter-flows similarity w.r.t. S selects a compact cluster with the most similar flows, corresponding to the set of attacking hosts. The obtained signature can be expressed as (nDsts == 1) ∧ (nSYN/nPkts > λ 3 ) ∧ (nPkts/sec > λ 4 ), which combined with the large number of identified sources (nSrcs > λ 5 ) confirms the nature of a SYN DDoS attack. This signature is able to correctly isolate the most aggressive hosts of the DDoS attack, i.e., those with highest packet rate.

Figures 1.(f,g) depict the detection of an ICMP flooding DoS attack. Traffic is aggregated in IPdst flows, thus the attack is now detected as an outlier rather than as a smallsize cluster. Absolute rules are not applicable in the case of outliers detection. Relative rules correspond to the separation of the outlier from the biggest cluster in each sub-space, which statistically represents normal-operation traffic. Besides showing typical characteristics of this attack, such as a high packet rate of exclusively ICMP packets from the same source host, both partitions show that the detected attack does not involve the largest elephant flows in the time slot. This emphasizes the ability of the algorithm to detect attacks that are not necessarily different from normal-operation traffic in terms of volume, but that they differ in other, less evident characteristics. The obtained signature can be expressed as (nICMP/nPkts > λ 6 ) ∧ (nPkts/sec > λ 7 ).

More evaluation results can be found at [START_REF] Casas | Sub-Space Clustering & Evidence Accumulation for Unsupervised Network Anomaly Detection[END_REF], including an evaluation of true-positives/false-alarm rates, as well as a comparison against other methods for unsupervised anomaly detection. These results confirm the outperforming ability of our approach in the unsupervised anomaly detection domain.

VI. COMPUTATIONAL TIME AND PARALLELIZATION

The last issue that we analyze is the Computational Time (CT) of the algorithm. The SSC-EA-based algorithm performs multiple clusterings in N (m) low-dimensional sub-spaces X i ⊂ X. This multiple computation imposes scalability issues for on-line detection of attacks in very-high-speed networks. Two key features of the algorithm are exploited to reduce scalability problems in number of features m and the number of aggregated flows n to analyze. Firstly, clustering is performed in very-low-dimensional sub-spaces, X i ∈ R 2 , which is faster than clustering in high-dimensional spaces [START_REF] Jain | Data Clustering: 50 Years Beyond K-Means[END_REF]. Secondly, each sub-space can be clustered independently of the other sub-spaces, which is perfectly adapted for parallel computing architectures. Parallelization can be achieved in different ways: using a single multi-processor and multi-core machine, using network-processor cards and/or GPU (Graphic Processor Unit) capabilities, using a distributed group of machines, or combining these techniques. We shall use the term "slice" as a reference to a single computational entity. we separately cluster each of the N = m(m -1)/2 subspaces X i , and take the worst-case of the obtained clustering time as a representative measure of the CT in a single subspace, i.e., CT(X SSCwc ) = max i CT(X i ). Then, if N M , we have enough slices to completely parallelize the SSC algorithm, and the total CT corresponds to the worst-case, CT(X SSCwc ). On the contrary, if N > M, some slices have to cluster various sub-spaces, one after the other, and the total CT becomes (N %M + 1) times the worst-case CT(X SSCwc ), where % represents integer division. The first interesting observation from figure 2.(a) regards the increase of CT(X) when m increases, going from about 8 seconds for m = 2 to more than 200 seconds for m = 29. As we said before, clustering in low-dimensional spaces is faster, which reduces the overhead of multiple clusterings computation. The second paramount observation is about parallelization: if the algorithm is implemented in a parallel computing architecture, it can be used to analyze large volumes of traffic using many traffic descriptors in an on-line basis; for example, if we use 20 traffic features and a parallel architecture with 100 slices, we can analyze 10000 aggregated flows in less than 20 seconds. Figure 2.(b) compares CT(X) against CT(X SSCwc ) for an increasing number of flows n to analyze, using m = 20 traffic features and M = N = 190 slices (i.e., a completely parallelized implementation of the SSC-EA-based algorithm). As before, we can appreciate the difference in CT when clustering the complete feature space vs. using low-dimensional sub-spaces: the difference is more than one order of magnitude, independently of the number of flows to analyze. Regarding the volume of traffic that can be analyzed with this 100% parallel configuration, the SSC-EA-based algorithm can analyze up to 50000 flows with a reasonable CT, about 4 minutes in this experience. In the presented evaluations, the number of aggregated flows in a time slot of ΔT = 20 seconds rounds the 2500 flows, which represents a value of CT(X SSCwc ) ≈ 0.4 seconds. For the m = 9 features that we have used (N = 36), and even without doing parallelization, the total CT is N × CT(X SSCwc ) ≈ 14.4 seconds.

VII. CONCLUSIONS

The completely unsupervised algorithm for detection of network attacks that we have presented has many interesting advantages w.r.t. previous proposals. It uses exclusively unlabeled data to detect and characterize network attacks, without assuming any kind of signature, particular model, or canonical data distribution. This allows to detect new previously unseen network attacks, even without using statisticallearning. By combining the notions of Sub-Space Clustering and multiple Evidence Accumulation, the algorithm avoids the lack of robustness of general clustering approaches, improving the power of discrimination between normal-operation and anomalous traffic. We have shown how to use the algorithm to automatically construct signatures of network attacks without relying on any kind of previous information. We claim that such an approach can be used do devise autonomous network security systems, in which the SSC-EA-based algorithm runs in parallel to any standard security device, producing specific signatures to unknown anomalous events.

Finally, and contrary to previous work on clustering for detection of network attacks, we have evaluated the computational time of our algorithm. Results confirm that the use of the algorithm for on-line unsupervised detection and automatic generation of signatures is possible and easy to achieve for the volumes of traffic that we have analyzed. Even more, they show that if run in a parallel architecture, the algorithm can reasonably scale-up to run in high-speed networks, using more traffic descriptors to characterize network attacks.

Figure 1 .

 1 Figure 1. Filtering rules for characterization of attacks in WIDE.

Figure 2

 2 depicts the CT of the SSC-EA-based algorithm, both (a) as a function of the number of features m used to describe traffic flows and (b) as a function of the number of flows n to analyze.

Figure 2 .

 2 (a) compares the CT obtained when clustering the complete feature space X, referred to as CT(X), against the CT obtained with SSC, varying m from 2 to 29 features. We analyze a large number of aggregated flows, n = 10 4 , and use two different number of slices, M = 40 and M = 100. The analysis is done with traffic from the WIDE network, combining different traces to attain the desired number of flows. To estimate the CT of SSC for a given value of m and M , we proceed as follows: first, complete Feature Space Distributed Sub-Space Clustering, 190 slices (a) Time vs. nº features. (b) Time vs. nº patterns.

Figure 2 .

 2 Figure 2. Computational Time as a function of nº of features and nº of flows to analyze. The number of aggregated flows in (a) is n = 10000. The number of features and slices in (b) is m = 20 and M = 190 respectively.

ACKNOWLEDGMENTS

This work is part of the ECODE project, funded by the European commission under grant FP7-ICT-2007-2/223936.


Inken Wohlers
email: wohlers@cwi.nl

Noël Malod-Dognin

Rumen Andonov

Gunnar W Klau

CSA: Comprehensive comparison of pairwise protein structure alignments

Keywords: Protein structure alignments, alignment comparison, web server, exact algorithm, scoring function

CSA is a web server for the comprehensive comparison of pairwise protein structure alignments. Its exact alignment engine computes either optimal, top-scoring alignments or heuristic alignments with quality guarantee for the inter-residue distance based scorings of contact map overlap, PAUL, DALI and MATRAS. These and additional, uploaded alignments are compared using a number of quality measures and intuitive visualizations. CSA brings new insight into the structural relationship of the protein pairs under investigation and is a valuable tool for studying structural similarities.

CSA : Comparaison compréhensible d'alignement de paires de structures de protéines

Résumé :

CSA est un serveur web pour la comparaison compréhensible des alignements de paires de protéines. Son moteur d'alignement exact calcule, pour les scores de similarité basés sur les distances inter-résidus CMO (Contact Map Overlap maximization), PAUL, DALI et MATRAS, soit des alignements optimaux (ayant le plus haut score), soit des alignements heuristiques avec une garantie de qualité. Ces alignements, plus ceux uploadés par l'utilisateur, sont comparés en utilisant de nombreuses mesures de qualité et des visualisations intuitives. CSA apporte un nouveau regard sur les relations structurales entre paires de protéines, et constitue un outil précieux pour l'étude des similarités structurales. CSA est disponible sur http://csa. project.cwi.nl Mots-clés : Alignement de structure protéiques, comparaison d'alignements, serveur web, algorithme exact, fonctions de score

Introduction

Protein structural alignment is a key method for answering many biological questions that involve the transfer of information from well-studied proteins to less well-known proteins. Since structures are more conserved during evolution than sequences, structural alignment allows for the most precise mapping of equivalent residues. It is notably important for (i) detecting and investigating structural motifs, functional sites, and common cores and (ii) measuring similarity between proteins and bringing them in evolutionary relationship, e.g., by classification. Numerous web servers are available that offer individual methods for computing structural alignments, e.g., [START_REF] Kawabata | MATRAS: a program for protein 3D structure comparison[END_REF][START_REF] Mosca | RAPIDO: a web server for the alignment of protein structures in the presence of conformational changes[END_REF][START_REF] Margraf | The SALAMI protein structure search server[END_REF][START_REF] Holm | Dali server: conservation mapping in 3D[END_REF].

Many structure-based scoring schemes have been proposed and there is no consensus which scoring is the best [START_REF] Hasegawa | Advances and pitfalls of protein structural alignment[END_REF]. Comparative studies find that scorings have individual strengths and weaknesses and that alignments produced by different methods can differ considerably [START_REF] Mayr | Comparative analysis of protein structure alignments[END_REF]. In the context of protein classification, there are first attempts to integrate information from alignments generated by different structural alignment methods, e.g., [START_REF] Camoglu | Integrating multi-attribute similarity networks for robust representation of the protein space[END_REF][START_REF] Barthel | ProCKSI: a decision support system for protein (structure) comparison, knowledge, similarity and information[END_REF].

Here, we present CSA (Comparative Structural Alignment), the first web server for comprehensive comparison of pairwise protein structure alignments at single residue level. CSA facilitates evaluating the agreement between and advantages of alignments that maximize different established scoring schemes. It offers the computation of alignments using the scoring schemes of dali [START_REF] Holm | Protein structure comparison by alignment of distance matrices[END_REF], contact map overlap (cmo) [START_REF] Godzik | Flexible algorithm for direct multiple alignment of protein structures and sequences[END_REF], matras [START_REF] Kawabata | Protein structure comparison using the Markov transition model of evolution[END_REF], and paul [START_REF] Wohlers | Towards optimal alignment of protein structure distance matrices[END_REF]. CSA uses our own, exact algorithm [START_REF] Andonov | Maximum contact map overlap revisited[END_REF][START_REF] Wohlers | Algorithm engineering for optimal alignment of protein structure distance matrices[END_REF] that can be used with any inter-residue distance based scoring scheme. We choose cmo and paul scoring since they are tailored to the algorithm and dali and matras scoring because they are established and their programs and web servers [START_REF] Kawabata | MATRAS: a program for protein 3D structure comparison[END_REF][START_REF] Holm | Dali server: conservation mapping in 3D[END_REF] are widely used. CSA returns an optimal, i.e. top-scoring alignment, if found within the time limit, or otherwise an alignment with a quality guarantee that specifies how much improvement is at most possible. We denote this by calling our program and its alignments dalix and matrasx, in which x indicates exact.

Optimality comes at the prize of higher running time, but is especially important when comparing alignments. A top-scoring, but biologically implausible, alignment implies that the scoring scheme used is inadequate to detect the given structural relationship and a different scoring might be more advisable. In the case of pairwise structural alignment, in which primarily residue correspondences are of interest, and only secondarily the obtained similarity score, comparing alignments optimized with respect to different criteria thus brings additional insight.

In CSA, computed or uploaded alignments can be explored in terms of many inter-residue distance-, RMSD-and sequence-based scores and quality measures and with intuitive visualizations such that agreements and differences between alignments are easy to grasp. The user can thus make educated decisions about the structural similarity of two proteins and, if necessary, post-process alignments by hand. Furthermore, a comparative analysis allows to differentiate between proteins with one clear-case alignment on which various scorings agree and proteins with ambiguous alignments for which it depends on the application which alignment is preferable.

Materials and methods

Structural alignment algorithm

The exact algorithm used in CSA is based on an integer linear programming (ILP) model of the structural alignment problem as described in [START_REF] Wohlers | Algorithm engineering for optimal alignment of protein structure distance matrices[END_REF]. Solutions to the ILP are generated using the approach from [START_REF] Andonov | Maximum contact map overlap revisited[END_REF]. The algorithm combines branch-and-bound and Lagrangian relaxation, which can be seen as an iterative double dynamic programming method. The mathematical model supports a generic scoring scheme with positive and negative structural scores, sequence scores and affine gap costs. Many different scoring functions are special cases of this general scheme. Currently, CSA supports dali [START_REF] Holm | Protein structure comparison by alignment of distance matrices[END_REF], cmo [START_REF] Godzik | Flexible algorithm for direct multiple alignment of protein structures and sequences[END_REF], matras [START_REF] Kawabata | Protein structure comparison using the Markov transition model of evolution[END_REF], and paul [START_REF] Wohlers | Towards optimal alignment of protein structure distance matrices[END_REF].

Webserver implementation

The architecture of the web server is divided in a processing layer that computes (C++) and evaluates (Python) alignments and an output layer, which generates W3C-validated HTML websites, interacts with the user and displays all information (PHP and Javascript). The interface between the two layers is a MySQL database.

The alignment engine for all our four currently supported scoring schemes is identical and implemented in C++ as a stand-alone program. User-adjustable parameters are the time limit of the computation, the maximum number of branch-and-bound nodes, and the number of Lagrangian iterations in each node. Furthermore, each scoring scheme has different parameters, for example, the use of C α or C β inter-residue distances.

Computed or user-uploaded (in FASTA format) alignments are read into a Python class and subsequently written to the MySQL data base. A second Python class handles the computation of different scores. It obtains the required structural information from the PDB files with the help of the Biopython package Bio.PDB [START_REF] Hamelryck | PDB file parser and structure class implemented in Python[END_REF]. Tasks related to superpositioning are also handled by this package. Visualizations of distance and distance difference matrices are generated using the Python Imaging Library.

The website functions have been implemented in separate modules, which makes it easy to integrate additional structural alignment methods. The modularity is illustrated by the use of a tab menu. All web server functions are extensively documented, which is denoted by a question mark next to the respective section titles or table headers. Additionally, a documentation puts instructions and explanations into context. Notably, we documented all structural alignment scorings that are used within CSA and we provide the corresponding formulas and references. In the output layer, structures and their superpositions are visualized in Jmol (http://www.jmol. org) and images are generated using the PHP package pChart (http://www.pchart.net/).

Case studies

We illustrate the functionality of CSA using two case studies which are accessible from its main page via the links "Example 1" and "Example 2".

Benefits of visualization and comparison

The first case study deals with two proteins from the SISY data set [START_REF] Mayr | Comparative analysis of protein structure alignments[END_REF][START_REF] Berbalk | Accuracy analysis of multiple structure alignments[END_REF], ubiquitin-binding protein CUE2 (PDB ID 1otr, chain A, 49 residues) and the CUE domain of activating signal cointegrator 1 complex subunit 2 (PDB ID 2di0, chain A, 71 residues). The proteins belong to the SISYPHUS [START_REF] Andreeva | SISYPHUS-structural alignments for proteins with non-trivial relationships[END_REF] alignment AL00088995 of homologous proteins containing a CUE domain. The CUE domain is composed of a three helical bundle and it consists of 41 residues. It binds ubiquitin and is involved in protein degradation.

After specifying PDB IDs and chains on the main page of CSA, the user is redirected to the CSA evaluation environment. It is organized in tabs for the following tasks: overview on the protein structures, computing alignments using cmo, paul, dali or matras scoring, upload of external alignments, and the comparison of alignments.

The Structures tab lists PDB IDs, PDB file names, selected chains and their lengths and amino acid sequences. Links to the PDB [START_REF] H M Berman | The Protein Data Bank[END_REF] and to iHOP [START_REF] Hoffmann | Implementing the iHOP concept for navigation of biomedical literature[END_REF] are access points for additional Inria information concerning the proteins and their function. Protein structures are visualized in Jmol. Their C α and C β distance matrices and contact maps are visualized.

We compute cmo, paul, dalix and matrasx alignments using the default options, i.e., with a time limit of 30 CPU s. The setup of all four result pages is identical. Exemplary, we consider the CMO alignment page; parts of it are displayed in Figure 1.

Bounds on alignment score and similarity. The section Optimized score lists the resulting scores: the raw score s(A, B) of proteins A and B (here, the number of common contacts), and a similarity score that normalizes the raw score with respect to the self-similarity of the two proteins computed as 2s(A,B) s(A,A)+s(B,B) . Our exact algorithm returns lower an upper bounds (LB and UB ) on the raw and similarity scores. Based on these bounds, the relative gap in percent, 100 • UB-LB LB , quantifies by how many percent an alignment can at most be improved. Such a quality guarantee helps to quickly determine the progress of the computation as well as the similarity of the two proteins. If two proteins are dissimilar, the relative gap tends to be large, but the upper bound on the similarity score tends to be low from the beginning on. Aligning 1otrA with 2di0A w.r.t. cmo yields 125 common contacts, and the corresponding similarity score on a scale from 0 to 1 is 0.751. The relative gap is 0%, indicating that the top-scoring alignment has been found.

Structural conservation and variation. The Alignment section displays the computed structural alignment. Residues are color-coded according to either SSE (helix, sheet, coil) as assigned by DSSP [START_REF] Kabsch | Dictionary of protein secondary structure: pattern recognition of hydrogen-bonded and geometrical features[END_REF] or to residue pair score contribution. The second color-coding denotes how well single residue pairs are structurally conserved given the current alignment, cf. Fig. 1. For the two proteins containing the CUE domain, this indicates that the first identically aligned leucines are structurally conserved, and in fact this position is part of a motif for binding ubiquitin that consists of an invariant proline and two highly conserved leucines [START_REF] S C Shih | A ubiquitinbinding motif required for intramolecular monoubiquitylation, the CUE domain[END_REF]. Pairs of aligned residues with low score contribution highlight structural variations. In the cmo alignment, the N-and C-terminal regions are little structurally conserved, as well as the residues in the region of the invariant proline within the CUE domain, because the proline is located in a turn. Such a visualization of residue score contribution can hint towards a manual modification of the alignment by removing aligned residues with low score. In fact, this is what happens in the top-scoring dalix alignment of 1otrA and 2di0A, in which the four C-terminal residues with low cmo score are excluded from the alignment.

Comprehensive alignment-related data. Additional to the alignment, CSA displays the aligned segments, both using sequential and PDB residue numbering, cf. Fig. 1. Numerous raw alignment-and similarity scores are listed, for example the number of aligned residues, sequence identity and root mean square deviation (RMSD). Furthermore, some statistics concerning the alignment computation are given. These are the number of residues and inter-residue distances considered during computation. They greatly influence the memory consumption of the algorithm: the more inter-residue distances are considered, the more memory is needed and typically the larger the running time. Using default values, cmo only considers distances smaller than 7.5 Å, paul considers distances smaller than 8.5 Å (for C α distances, 9.5 Å), matras uses distances up to 50 Å and dali all distances. Because server memory is shared among users, we currently restrict computations using the dali or matras scoring to protein pairs with average length less than 150 residues. The allocation time for setting up all data structures is given, as well as the time actually spent on computing the alignment. The number of visited branch-and-bound nodes gives a good estimate on the progress of the computation. The proteins are superposed according to the alignment and visualized in Jmol. The trace of aligned residues and the distance difference matrix is plotted.

We upload an additional alignment in the tab for the first custom alignment. This alignment aligns only the 38 residues that belong to the respective CUE domain and that are structurally equivalent according to SISY. Furthermore, we upload a second custom alignment which has been generated by the dali server [START_REF] Holm | Dali server: conservation mapping in 3D[END_REF]. The dali server uses a heuristic algorithm to find a good alignment according to the dali score.

Improving, verifying optimality and assessing quality of heuristic alignments. Many different scorings and quality measures can be compared in the Comparison tab: the cmo, paul, dali, and matras raw and similarity scores, dali z-score [START_REF] Holm | Dictionary of recurrent domains in protein structures[END_REF], TM-score [START_REF] Zhang | Scoring function for automated assessment of protein structure template quality[END_REF], number and percentage of aligned residues, coordinate and distance RMSD, RMSD100 [START_REF] Carugo | A normalized root-mean-square distance for comparing protein three-dimensional structures[END_REF], and sequence identity. For 1otrA and 2di0A, all six computed and uploaded alignments differ from each other. While cmo and paul alignment were computed to optimality in less than a second, the dalix alignment has the potential to be improved by up to 12% and the matrasx alignment by up to 24%. We also observe that the alignment that was computed by the dali server and then uploaded is better with respect to dali score than the alignment computed by our exact algorithm within 30s. We thus increase the maximum running time for dalix and matrasx to 10 minutes. Now, both alignments are computed to provable optimality and our top-scoring dalix alignment slightly improves the heuristic solution returned by the dali server. dalix and matrasx alignments thus can be used to obtain quality guarantees for dali or matras alignments and in some cases also to either proof their optimality or to compute a better alignment.

Multi-criteria comparison and selecting a sound alignment. Alignment trace comparison as introduced in [START_REF] Godzik | The structural alignment between two proteins: is there a unique answer?[END_REF] gives a visual overview about agreements and differences between alignments. Here, any subset of alignments can be shown. Using this visualization, we find that all alignments (except the SISY reference, which excludes 3 residues in the center of the domain) correctly align all 41 residues of the CUE2 domain, and that they differ in aligning the neighboring N-and C-terminal residues. A radar chart compares the different scores, cf. Fig. 2. This chart helps to quantify score differences and allows to decide whether one alignment is clearly preferable, i.e., better with respect to all criteria. The chart also allows to make an intuitive decision which alignment is most appropriate in cases in which different scorings disagree as it is the case for 1otrA and 2di0A. Here, intuitively the dalix alignment is the best choice since it performs good or best according to all criteria.

Two residue pair lists show aligned residues that appear in all, resp. in the majority, of the alignments. They each constitute a consensus alignment. In the case of aligning 1otrA and 2di0A, we see that such a consensus is useful: all alignments only agree in aligning the CUE2 domain. The consensus thus highlights the structurally conserved and biologically relevant region of the alignment.

Alignment of flexible proteins

We illustrate the usefulness of comparing structural alignments in the case of protein flexibility. This is a challenge for most structural alignment methods because flexible proteins typically do not superpose well unless the flexibility is accommodated for, e.g., by explicitly introducing a hinge.

Comparing flexible and rigid scoring schemes. We align two conformations of the calmodulin protein (PDB IDs 4cln, chain A and 2bbm, chain A, with a length of 148 residues). In Inria structure 4clnA, calmodulin is bound to a ligand, in 2bbmA it is unbound. In bound conformation, a central helix is split and the components at the ends of the helix are moved towards each other. We align the two conformations using cmo and paul. We furthermore upload the alignments computed by tm-align [START_REF] Zhang | TM-align: a protein structure alignment algorithm based on the TM-score[END_REF], an algorithm maximizing the TM-score, and dast [START_REF] Malod-Dognin | Maximum cliques in protein structure comparison[END_REF], a local structural alignment method that determines the longest alignment with distance RMSD less than 4 Å. We find that both cmo and paul correctly align the two conformations over their entire length. To keep the comparison concise, we thus exclude the cmo alignment from the following visualizations. Figure 3 displays the two conformers superposed according to the TM-alignment and the alignment trace comparison. While paul aligns all residues of the two conformers correctly, tm-align aligns only the C-terminal, rigidly superposable region (except the C-terminal residue). dast also aligns the C-terminal region, but excludes and shifts further residues from the alignment. The radar chart comparing the different scores as well as the distance difference matrices displayed in Figure 3 show why: while cmo, paul, dali and matras scoring by far favor the alignment of the entire conformers, TM-score as well as RMSD100 clearly favor the TM-and dast alignment, which has a much smaller RMSD, but aligns only the C-terminal region.

Detecting flexibility and hinges. For each alignment we display the distance difference matrix. This is a symmetric square matrix with entries |d A ijd B ij | at position (i, j), where i is the i-th aligned position and j the j-th aligned position. Here, distance differences are visualized using a color gradient in which 0 Å is colored red, 2.5 Å green 5 Å blue. Regions with low inter-residue distance differences correspond to rigidly superposable fragments. For the paul alignment of 4clnA and 2bbmA, red blocks in the distance difference matrix indicate that both the N-terminal and C-terminal regions can be superpositioned very precisely. The distance differences between these two regions, however, are large, denoted by the blocks in blue color. The two regions can thus only be well superpositioned individually. A hinge is present at the residue bordering the two blocks (position 80) [START_REF] Emekli | HingeProt: automated prediction of hinges in protein structures[END_REF]. tm-align and dast align only the C-terminal region, thus avoiding any large distance differences. dast is more restrictive in excluding large distance differences, it does not align a few residues that are still aligned by the TM-alignment and which have distance differences larger than 5 Å, colored in blue.

Scores as cmo and paul, which implicitly ignore RMSD, are useful to gain information about flexible regions. While this feature is beneficial for flexible proteins it may also introduce flexibility where this is not appropriate. Protein similarities consisting in compact, well superposable fragments are therefore often better detected by maximizing scores like the TM-or the dast score.

Conclusion

Different structural alignment scoring functions have different strengths and weaknesses. Which scoring to use depends on the application and on the structural relationship of the investigated proteins. Their different focus on handling various aspects of structural similarity is one reason why there are many different structural alignment scorings and programs and no consensus which combination is best.

We therefore consider it beneficial to compute alignments using different scoring schemes and algorithms and to compare them in order to gain insight into their structural relationship. The CSA web server provides the tools for such a comparison. CSA allows to compute alignments with various scorings, returns a quality guarantee for the alignments and enables the user to additionally evaluate and compare uploaded alignments. In the most common case in which scorings and alignments disagree, it facilitates evaluating the agreement and differences between them and selecting the most suitable alignment.

Funding

This work was partly supported by the DFG [KL 1390/2-1]. The closer a point is to 1, the better the corresponding score. cmo, paul, dalix and matrasx alignments have been computed by our exact algorithm and are provably optimal concerning their respective score. The SISY reference alignment aligns 38 residues of the CUE2 domain. The dali alignment was computed by the dali server and has slightly lower dali score than the optimal dalix alignment. The reference alignment is far behind in all scores except RMSD100 and TM-score, for which it performs quite well. The matrasx alignment performs especially poor for these two measures. Intuitively, the dalix alignment is most preferable since it has optimal dali and close to optimal cmo, paul and matras scores, as well as the best TM-score and RMSD100.

Acknowledgments

Figure 1 :

 1 Figure 1: Parts of the information displayed on the website for the cmo alignment of 1otrA and 2di0A.

Figure 2 :

 2 Figure 2: A radar chart for comparison of alignment scores for six different alignments of 1otrA and 2di0A.The closer a point is to 1, the better the corresponding score. cmo, paul, dalix and matrasx alignments have been computed by our exact algorithm and are provably optimal concerning their respective score. The SISY reference alignment aligns 38 residues of the CUE2 domain. The dali alignment was computed by the dali server and has slightly lower dali score than the optimal dalix alignment. The reference alignment is far behind in all scores except RMSD100 and TM-score, for which it performs quite well. The matrasx alignment performs especially poor for these two measures. Intuitively, the dalix alignment is most preferable since it has optimal dali and close to optimal cmo, paul and matras scores, as well as the best TM-score and RMSD100.

Figure 3 :

 3 Figure 3: Top left: The two calmodulin conformers (PDB IDs 4cln and 2bbm) superpositioned according to the TM-alignment, which aligns only one of the two regions that move relative to each other. Top right: Comparison of the alignment traces. Each axis corresponds to one conformer. Black boxes denote residue pairs aligned by all three scorings, paul, tm-align and dast. Light gray denotes residue pairs aligned by only one scoring. An intermediate shade of gray denotes agreement of two scorings. paul aligns all residues of the two conformers, tm-align and dast the C-terminal region. Center: The radar chart illustrates the difference between scorings that are more in favor of a flexible alignment, i.e. cmo, paul, dali and matras, and scorings that are more in favor of a rigid superposition of low RMSD. Bottom: The distance difference matrices illustrate the difference between the flexible paul alignment, that aligns all residues in spite of large distance differences (colored blue), and the TM-and dast alignment, which exclude large distance differences, but only align the C-terminal region.

RR n°7874

We thank Maarten Dijkema for his help with setting up the server. Further, we thank Mathilde Le Boudic-Jamin for extensively testing CSA and Thomas Metzler for his idea of using a radar chart for comparison. 6.0.1 Conflict of interest statement. None declared.


Laurent Houssin
email: houssin@laas.fr

Christian Artigues
email: artigues@laas.fr

Erwan Corbel
email: erwan.corbel@thalesaleniaspace.com

Frequency allocation problem in a SDMA satellite communication system

Keywords: SDMA system, frequency allocation problem, radio resource management

Introduction

Satellite communication systems move towards greater capacity, higher flexibility (with respect to the position of the users) and better service to the end-user. SDMA (Spatial Division Multiple Access) appears to be one way to achieve these requirements at the same time [START_REF] Liberti | Smart Antennas for Wireless Communications[END_REF]. SDMA is a principle of radio resource sharing that relies on the division of the space dimension into separated communication channels. It can be used with common Frequency Division Multiple Access (FDMA), Time Division Multiple Access (TDMA) or Code Division Multiple Access (CDMA) techniques. Main future terrestrial communication standards (such that WIMAX, 3GPP, LTE) implement SDMA. SDMA basically relies on adaptive and dynamic beam-forming associated to a clever algorithm in charge of resource allocation. The satellite beam-former optimizes the antenna diagram with respect to the positions of the users in order to maximize the gain while mitigating interferences. The resource allocation algorithm carefully designs a frequency plan that • prevents or limits interferences between users, • tailors the allocated bandwidth to the user need in order to save the spectrum.

Today SDMA is currently used by IRIDIUM system in L-band, a constellation of 66 Low-Earth Orbit satellites, thanks to time beam-switching. SDMA is also foreseen as a key enabling technique to increase the capacity of future two-way satellite communications systems in lowfrequency bands (typically lower than 5-6 GHz) through the interference mitigation and high frequency reuse [START_REF] Corbel | Hybride satellite and terrestrial mobile network for 4g: Candidate architecture and space segment dimensioning[END_REF].

978-1-4244-4136-5/09/$25.00 2009 IEEE It is also expected to play an important role in future systems devoted to Public Protection and Disaster Relief (PPDR) and Global Monitoring for Environment and Security (GMES) missions which require on-demand beamforming [START_REF] Gayrard | The role of satellites in broadband wireless access for public protection, disaster relief and gmes missions[END_REF].

The satellite telecommunication system that we study in this paper aims at establishing bi-directional communications involving a gateway connected to a terrestrial network and some user terminals located in a service area. This paper studies static Frequency Allocation Problems (FAP) in this system and two scenarios are considered: one based on SDMA and the other based on "traditional" spot coverage. We propose original Integer Linear Programming (ILP) formulations and greedy allocation algorithms for these problems. The difficulty for solving the FAP is increased by considering cumulative interference constraints.

We then compare the performance of the two scenarios.

The remainder of the paper is organized as follow: Section 2 is dedicated to the telecommunication system and the description of the scenarios. In Section 3, we present a cumulative formulation of the FAP interference constraints. In Section 4, ILP formulations and greedy algorithms proposed for both scenarios are described. Section 5 presents the results obtained by the different algorithms on the scenarios. Concluding remarks are drawn in Section 6.

Telecommunication system and scenarios for frequency allocation

System description

The service area is a rectangular grid where users position is uniformly distributed and where beams are directed. The grid size is u = [-0.043980, 0.048870] and v = [-0.021152, 0.012702] (Cartesian coordinates are considered although the same study can be done with spherical coordinates) with a step of 5 × 10 -4 . The satellite orthogonal projection onto the service area defines point (0, 0) of the grid.

Beams have two particular characteristics which are the direction that we consider through the position of the beam center in the service area, and the radiation pattern. An analytic representation of the radiation pattern enables to compute the directive gain of the antennas for a considered direction. This description is such that

G Sat (u, v, u 0 , v 0 ) = G 1 ×G 2 (u, v, u 0 , v 0 )×G 3 (u, v, u 0 , v 0 ) (1)
with

G 1 = η πD λ 2 , G 2 (u, v, u 0 , v 0 ) =   2J 1 πD λ (u -u 0 ) 2 + (v -v 0 ) 2 πD λ (u -u 0 ) 2 + (v -v 0 ) 2   2 and G 3 (u, v, u 0 , v 0 ) = 2J 1 πd λ √ u 2 + v 2 πd λ √ u 2 + v 2 2 ,
where J 1 (x) is the Bessel functions of the first kind. We use the notation u, v Cartesian coordinates of the user terminal; The term G 1 corresponds to the maximum gain antenna whereas G 2 depends on the distance between the user and the beam center, i.e. (u -u 0 ) 2 + (v -v 0 ) 2 and G 3 depends on user position related to the satellite, i.e.

√ u 2 + v 2 . The product G 2 (u, v, u 0 , v 0 ) × G 3 (u, v, u 0 , v 0 )
is illustrated in fig. 1 where the left pattern is the pattern of beam centered in (0.03,-0.03) and the right one is centered in (0,0), i.e directly under the satellite. It emphasizes that N RsM odCod that depends on the modulation scheme and the code scheme (the digital communication scheme considered here is a phase key shifting). More precisely, a user will be served if

C N + I ≥ C N RsM odCod (2)
where C N +I is the user Signal-to-Interference-plus-Noise Ratio (SINR) which determines quantitatively the signal quality.

The link budget enables to compute the SINR :

C N +I -1 = C N -1 F eeder + C I -1 F eeder + C IM -1 + C N -1 U ser + C I -1 U ser , (3) in which C N U ser = K 1 × G Sat (U ser Beam → U ser) K 2 and C I U ser = G Sat (U ser Beam → U ser) j∈Interf (G Sat (U ser Beam → U ser(j))
.

Terms K 1 and K 2 involve technical parameters, such as the atmospheric loss, the antenna temperature and the Equivalent Isotropically Radiated Power (EIRP) which are considered as constant. Consequently, C N U ser and C I U ser only depends on user position and beam center position. Indeed, G Sat (U ser Beam → U ser) describes the gain for a user terminal and its beam (see equation ( 1)). The set Interf is the set of users sharing the same channel. Moreover, C N f eeder and C I f eeder which are gateway characteristics, and C IM that is a satellite characteristic are also constant. Specifications of the system are the following.

• Beams are only adaptive in direction and not in shaping (although it is technically feasible).

• For the first scenario (involving SDMA), each user has a beam directly centered on him. It leads to consider as many beams as users in the system.

• We only focus on the case of users-to-feeder (gateway) link.

Scenarios description

In the users-to-feeder link, interferences can occur in the link budget when several users share the same frequency.

We illustrate this phenomena in fig. 2 where black diamonds are users sharing the same frequency. It shows that users u 3 and u 2 are interferers for user u 1 .

Fig. 2: Users-to-feeder link.

In the first scenario, named "scenario 1", we consider that each user has a beam directly centered on him (which is possible thanks to SDMA). The number of available channels is 8. In this case,

G 2 (u, v, u 0 , v 0 ) is always equal to the maximum value (that is G 2 (u, v, u 0 , v 0 ) = 1) since (u, v) = (u 0 , v 0 ).
Contrary to scenario 1, the second scenario involves fixed beams, however the frequency assignement is variable according to the demand (the channel number can be adjusted for a beam). The service area is composed of 40 fixed beams which form a spot-based coverage. We also have 8 available channels. Consequently, (since no channel can be used more than once in a spot) we can not serve more than 8 users in a spot.

For both scenario, the FAP considered in this paper consists in finding and inteference-free frequency allocation to the users maximizing the number of served users. This problem is static: the set of users is known in advance and we do not take dynamic arrivals and departures of users into account.

A cumulative interference representation

In a FAP where channels are limited, results depend on the ability of the system to allocate the same channel to several users.

In this section, we show that it is possible to obtain a cumulative representation of interferences from equations ( 2) and (3).

Involving the link budget, constraint (2) for a user i, becomes

1 A + B i + P j∈Interf (i) G Sat (U ser Beam(i)→U ser(j)) G Sat (U ser Beam(i)→U ser(i))
≥ D

where

A = C N -1 F eeder + C I -1 F eeder + C IM -1 , B i = C N -1 U seri = K 2 K 1 × G Sat (U ser Beam(i) → U ser(i)) and D = (C/N ) RsM odCod .
Previous inequation leads to

1 ≥ AD + DB i + P j∈Interf (i) DG Sat (U ser Beam(i)→U ser(j)) G Sat (U ser Beam(i)→U ser(i))
and finally, the cumulative representation follows

j∈Interf (i) DG Sat (U ser Beam(i) → U ser(j)) ≤ G Sat (U ser Beam(i) → U ser(i)) (1 -AD -B i D).
For scenario 1, we write last constraint

j∈Interf (i) δ ij ≤ α i , (4)
where

α i = G Sat (U ser Beam(i) → U ser(i)) (1-AD -B i D) and δ ij = j∈Interf (i) DG Sat (U ser Beam(i) → U ser(j)) .
Let the system with n users, α is a n-row vector and δ is a n × n matrix.

We can deduce that α i represents the maximum level of interferences that the user i can support. In this way, δ ij describes the interference level of user j if users j and i share the same frequency. The cumulative representation is motivated by a linear representation (although constraint (2) was not linear) that enable the use of integer linear programming (see §4.2).

Concerning scenario 2, we denote the interference inequality by

j∈Interf (i) γ ij ≤ β i . ( 5
)
It is worth emphasizing that α i and β i are different since the center of the beam related to user i does not correspond to user i coordinates in scenario 2 whereas in scenario 1, it does. The same remark can be done for δ ij and γ ij .

4. Modeling and solving scenarios 1 and 2 FAP

FAP literature overview

Most approaches dealing with interference minimization FAP consider binary interference constraints, i.e. involving only two users. Because of the strong links between graph coloring and frequency allocation with binary interference constraints, most methods found in the literature are inspired by coloring algorithms. We also know unfortunately the graph coloring problems, and consequently the FAP, are NP-hard. Among the proposed methods, the constructive (greedy) algorithms are widely used since they are simple, fast and also able to solve dynamic FAP. In this category, we find the generalization of DSATUR procedure [START_REF] Brélaz | New methods to color the vertices of a graph[END_REF]. Other more sophisticated algorithms, such as local search, metaheuristics, ILP and constraint programming approaches, are frequently encountered, see [START_REF] Aardal | Models and solution techniques for the frequency assignment problem[END_REF] for a comprehensive review of the state-of-the-art methods for the FAP with binary interference constraints.

One of the difficulties appearing in the telecommunication system considered in this study (for both scenarios 1 and 2) lies in the explicit consideration of non-binary interference constraints. In terms of graph coloring, deciding whether a given coloring is feasible or not cannot be made anymore by checking pairwise user color assignments. Instead, for a given user, the cumulative interferences of the users assigned to the same color has to be computed. Then, the coloring is feasible if this cumulative interference remains under a user-dependent threshold (see Section 3). In the literature, only a few approaches take account explicitly of such interferences for frequency assignment [START_REF] Dunkin | Towards high order constraint representations for the frequency assignment problem[END_REF][START_REF] Mannino | An enumerative algorithm for the frequency assignment problem[END_REF][START_REF] Palpant | Models and methods for frequency allocation with cumulative interference constraints[END_REF]. This study is partly based on integer linear programming formulations proposed in [START_REF] Palpant | Models and methods for frequency allocation with cumulative interference constraints[END_REF].

Integer linear programming formulations (ILP)

Taking account of hypothesis and simplifications presented in Section 2, FAP corresponding to scenarios 1 and 2 can be described as coloring problems and thus formalized as the corresponding combinatorial optimization problems. Each user has to be assigned a color, representing the allocated carrier.

For the SDMA FAP, (scenario 1), the following data are considered. n denotes the number of users. U = {1, . . . , n} is the set of users. C is the number of colors (channels). α i denotes the interference threshold for user i.

δ ij is the interference component from user j on user i, if i and j are assigned the same color.

Binary decision variables x ic are defined for i ∈ {1, . . . , n} and c ∈ {1, . . . , C} with n, the number of users and C the number of available colors. x ic = 1 if color c is allocated to user i and x ic = 0 otherwise. The problem can be represented by the following ILP:

max n X i=1 C X c=1 x ic (6) C X c=1 x ic ≤ 1 i = 1, . . . , n (7) n X j=1 δ i,j x jc ≤ α i + M i (1 -x ic ) i = 1, . . . , n c = 1, . . . , C (8)
x ic ∈ {0, 1} i = 1, . . . , n c = 1, . . . , C (9)
Objective ( 6) consists in maximizing the number of accepted users. Constraints [START_REF] Liberti | Smart Antennas for Wireless Communications[END_REF] state that at most one color has to be selected for each user. Constraints ( 8) are the cumulative interference constraints. They represent, in case color c is allocated to user i, the respect of the threshold for user i taking account of users that are assigned color c, i.e. possible interferers. Constant M i has to be large enough to withdraw the constraint if i is not assigned color c (x ic = 0). More precisely, we set

M i = n j=1 δ ij -α i .
For the fixed beam FAP (scenario 2), the data are similar with additional features concerning the spots (a spot designing the area covered by a given beam): m denotes the number of spots. S = {1, . . . , m} is the set of spots. U s is set of users covered by spot s ∈ S. β i denotes the interference threshold for user i. γ ij is the interference component from user j on user i if i and j are assigned the same color.

We define the fixed beam FAP as a combinatorial optimization problem resembling the SDMA FAP preventing two users covered by the same spot from being assigned the same non-zero color. We obtain the following ILP:

max n X i=1 C X c=1 x ic (10) C X c=1 x ic ≤ 1 i = 1, . . . , n (11)
X i∈Us x ic ≤ 1 s = 1, . . . , m c = 1, . . . , C (12) n X j=1 γ i,j x jc ≤ β i + N i (1 -x ic ) i = 1, . . . , n c = 1, . . . , C (13)
x ic ∈ {0, 1} i = 1, . . . , n c = 1, . . . , C (14)
Scenario 2 ILP differs with the scenario 1 ILP via constraints (12) and (13). Constraints (12) prevent any color from being allocated more than once in a given spot. Hence a maximum number of C users may be served in the same spot. In constraints (13), values β i and γ ij differ from α i and δ ij since beams are not centered on the users (see Section 3). As for scenario 1, N i has to be large enough to ensure the constraint is verified when color c is not allocated to user i. For this purpose, we set N i = n j=1 γ ij -β i . Both above-defined ILP can be solved by an integer linear programming solver, via branch and bound.

Greedy algorithms

Solving the ILP formulations provides optimal solutions only for small problems. For large-sized problems it is necessary to use a heuristic. We propose greedy algorithms to solve scenarios 1 and 2 FAP. For both scenarios, the principle of the greedy algorithms is, first, to consider the users sequentially according to a given criterion named the user priority rule. Second, either the selected user is assigned a frequency or rejected according to a second criterion, the frequency priority rule. Let Q denote the set of users that have not been assigned a color yet. Initially, we have Q = U . At each step of the greedy algorithm, a user i is removed from Q and is either rejected or assigned a color. We describe hereafter the rules selected for scenario 1.

For both priority rules, we use the frequency margin, where the margin M (i, c) of a user i ∈ Q for a color c is given by M (i, c) = α i -j∈U \Q∪{i},c(j)=c δ ij where c(j) denote the color allocated to user j. Namely, this margin corresponds to the positive or negative slack of the cumulative interference constraint for user i terminal if it is assigned color c.

The principle of the user priority rule is to select first the most constrained users in terms of available colors, as for the well-known DSATUR algorithm for standard graph coloring problems. The proposed rule selects the user having the smallest number of available colors, a color c being available for user i ∈ Q if M (i, c) ≥ 0 and if for all users j ∈ U \ Q having already been color c, M (j, c) ≥ 0. In case of a tie, we select the user having the smallest total margin for all its available colors.

The frequency priority rule selects a frequency for the selected user with the aim to minimize the impact of the assignment on other users. To that purpose, we simulate the assignment of each available color for the selected user and we compute the margin of the users that have not been assigned a color yet. Let N (i, c) denote the number of users j ∈ Q such that M (j, c) ≥ 0 if color c is allocated to the selected user i. The proposed rule select the colors that yields the larger N (i, c). In case of a tie, we select the color c that maximizes the sum of margins M (j, c) for all users j ∈ Q.

For scenario 2, the above user and frequency priority rules can be applied with the notable difference that color availability for a user i ∈ Q in a spot s has to consider the other users in U s . As soon as a user in s has been assigned a color c, then c becomes unavailable for the other users in U s .

Computational experiments and simulations

The ILP formulations have been solved using ILOG CPLEX 11.1 [START_REF] Ilog | ILOG CPLEX C++ API 11[END_REF] and the greedy algorithms have been coded in C++. We tested the proposed algorithms with C = 8, increasing stepwise the numbers of users by 20 from 20 to 200 users (which corresponds to reuse rates from 20/8=2.5 to 200/8=25). For each number of users, a set of 100 FAP data instances was obtained by randomly generating the user positions on the service area. The results were obtained on a Intel Core 2 Duo processor with 2,33GHz. The CPU time for the ILP resolution has been limited to 60 seconds afterwhich the best obtained integer solution is returned. The CPU times for the greedy algorithms were negligible. Fig. 3 clearly shows that the best results in terms of quality of service using the proposed algorithms are obtained for scenario 1 (SDMA-based system). ILP-based algorithms obtain significantly better results than the simple greedy algorithms for both scenarios. However, the greedy algorithm for scenario 1 performs better than the ILP-based algorithm for scenario 2. This numerical example with real parameters for the simulation shows qualitatively the benefits of scenario 1. We can also note that beyond 80 users, which correspond to a reuse rate of 10, differences between algorithms increase.

Conclusion

In this paper we have developed some integer linear proramming formulations involving cumulative interferences. The cumulative approach enables to take into account the non linear characteristics of interferences. Combining the SDMA system and the cumulative approach, we proved in section 6 the efficiency of scenario 1 with the two algorithms. Even better results and system optimizations using SDMA could be obtained by allowing shifts of the beam centers around the users and adjusting the EIRP parameters, yielding as a counterpart harder FAP. These features, together with consideration of dynamic aspects, constitutes the basis for further research.

u 0 , v 0

 0 Cartesian of the beam center; η antenna gain; D antenna diameter; d primary source diameter; λ wavelength.

Fig. 1 :

 1 Fig. 1: Beam patterns

Fig. 3

 3 Fig.3 displays, for each scenario/algorithm/number of users, the average number of accepted users in the computed frequency allocation plans.

Fig. 3 :

 3 Fig. 3: Average numbers of accepted users obtained by the proposed algorithms