Yan Kong

Weiming Dong
email: wmdong@nlpr.ia.ac.cn

Xing Mei
email: xmei@nlpr.ia.ac.cn

Xiaopeng Zhang
email: xpzhang@nlpr.ia.ac.cn

Jean-Claude Paul

J.-C Paul

SimLocator: robust locator of similar objects in images

Keywords: Similar objects, Object descriptor, Object matching, Stable locations, Contour extraction

Similar objects commonly appear in natural images, and locating and cutting out these objects can be tedious when using classical interactive image segmentation methods. In this paper, we propose SimLocator, a robust method oriented to locate and cut out similar objects with minimum user interaction. After extracting an arbitrary object template from the input image, candidate locations of similar objects are roughly detected by distinguishing the shape and color features of each image. A novel optimization method is then introduced to select accurate locations from the two sets of candidates. Additionally, a mattingbased method is used to improve the results and to ensure that all similar objects are located in the image. Finally, a method based on alpha matting is utilized to extract the precise object contours. To ensure the performance of the matting operation, this work has developed a new method for foreground extraction. Experiments show that SimLocator is more robust and more convenient to use compared to other more advanced repetition detection and interactive image segmentation methods, in terms of locating similar objects in images.

Introduction

Similar image elements commonly appear in both natural and artificial scenes, which provide multiple and non-local records of the same data. These similar objects may be randomly scattered in an image, and their spatial distribution and appearance variance reinforce the image's visual effect. Moreover, using identical or similar objects that are repeated in a pattern is a standard graphic design concept [START_REF] Krages | Photography: The Art of Composition[END_REF], being a prominent compositional feature in many photographs. Therefore, locating similar objects in an image is useful in understanding the aesthetics and mechanism of a particular scene in the real world [START_REF] Thompson | On Growth and Form[END_REF]. Locating and extracting such kind of redundancy accurately and efficiently can also be exploited to consolidate image edits in a non-local fashion [START_REF] Cheng | RepFinder: finding approximately repeated scene elements for image editing[END_REF].

Previous repeated elements locating or cutout methods used solely geometric [START_REF] Cheng | RepFinder: finding approximately repeated scene elements for image editing[END_REF] or color [START_REF] Huang | RepSnapping: efficient image cutout for repeated scene elements[END_REF] similarity as criteria for matching. However, these methods fail to detect accurately similar objects with obvious variance in illumination, outer shape, inner textures, or colors (e.g., Fig. 1). Interactive image segmentation systems, such as GrabCut [START_REF] Rother | grabcut": interactive foreground extraction using iterated graph cuts[END_REF], Lazy Snapping [START_REF] Li | Lazy snapping[END_REF] and Paint Selection [START_REF] Liu | Paint selection[END_REF], are also not suitable for segmenting similar objects because excessive manual operations may be needed when the number of objects is large, or when these objects overlap each other.

In this paper, we present a new method, namely, SimLocator, to analyze similar patterns in an image. After segmenting an object from the input interactively, we incorporate the shape and color features into a template matching framework to identify possible similar object locations within the scene. We then develop a joint searching scheme to determine high confidence locations within the two candidate sets, and utilize this information to calculate the final object locations. In some cases, the result can be improved further by a matting-based method to obtain a more complete locating result. After obtaining the object locations, a new method to obtain the precise object regions is utilized to refine the contours of the objects and cut them out from the scene. Thus, our algorithm can locate similar objects that may be dramatically different in colors, textures, or even shapes, rather than the very approximately same patterns in RepFinder [START_REF] Cheng | RepFinder: finding approximately repeated scene elements for image editing[END_REF].

Related work

Various object detection methods that utilize color-, texture-, or shape-based features have been previously presented. Leung and Malik [START_REF] Leung | Detecting, localizing and grouping repeated scene elements from an image[END_REF] proposed a graph-based algorithm to extract repeated units by growing elements using a graph. Berg et al. [START_REF] Berg | Shape matching and object recognition using low distortion correspondences[END_REF] achieved recognition by means of deformable shape matching and identification of correspondences among feature points. Ahuja and Todorovic [START_REF] Ahuja | Extracting texels in 2.1d natural textures[END_REF] detected natural texels by searching within a segmentation tree. Unfortunately, this method is too slow for interactive applications, and is limited to examples imaged from a viewing direction, which is nearly along the surface normal. Usually, local feature descriptors, such as shape contexts [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF], SIFT [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] and SURF [START_REF] Bay | Speeded-up robust features (surf)[END_REF], are used for object detection, and these factors can reliably match different views of an object or a scene, but cannot capture high-level scene structure. Pauly et al. [START_REF] Pauly | Discovering structural regularity in 3d geometry[END_REF] presented a method to discover regular or repeated geometric structures in 3D shapes. Automatic colorbased affine covariant region detectors, such as MSCR [START_REF] Forssén | Maximally stable colour regions for recognition and matching[END_REF], are suitable for detecting simple objects located in a smooth background. However, stochastic object distributions, overlapping, or subtle shape and color variations render the above methods unsuitable.

Given that the same or similar objects appear in a set of images, co-segmentation methods can segment a particular object from the set of images simultaneously. Joulin et al. [START_REF] Joulin | Discriminative clustering for image co-segmentation[END_REF] combined bottom-up image segmentation tools within a discriminative clustering framework to assign all images with the labels of either foreground or background. Meng et al. [START_REF] Meng | Object co-segmentation based on shortest path algorithm and saliency model[END_REF] formulated the co-segmentation problem as the shortest path problem based on local region similarities and saliency models. Kim et al. [START_REF] Kim | A hierarchical image clustering cosegmentation framework[END_REF] solved the co-segmentation problem by using hierarchical clustering. Rather than segmenting similar objects from multiple images, the method presented in the current paper focuses on locating similar elements in a single image.

Only a few studies have focused on the detection of general repeated or similar objects. Cheng et al. [START_REF] Cheng | RepFinder: finding approximately repeated scene elements for image editing[END_REF] presented a user-assisted approach in identifying approximately repeated scene elements in an image. Using boundary band matching to locate possible elements and employing active contours to obtain object boundaries, this efficient system is convenient for reconstructing scene structures. Nevertheless, as their method relies on similarity of boundary maps, illumination and inner pattern variations limit the accuracy of their results. On the other hand, the prerequirement of foreground matting also limits the practical usage of RepFinder, as separating the foreground from the image through simple user interactions can be difficult in many cases. Schweitzer et al. [START_REF] Schweitzer | A dual-bound algorithm for very fast and exact template matching[END_REF] presented a fast template matching technique that uses both lower and upper bounds on the match measure to compute for the K best matches. However, their method, which employs the Walsh transform based matching algorithm that is very sensitive to inner texture and illumination variations, is not fit for the detection of similar objects in a natural image. Moreover, the requirement of user-indicated number of objects also limits the application of this method. Huang et al. [START_REF] Huang | RepSnapping: efficient image cutout for repeated scene elements[END_REF] presented a graphbased method to cut out repeated elements from a scene. The limitation of their method is its strong dependence on very similar colors of objects. Furthermore, their method experiences difficulties in detecting accurate location and contours of each object, especially for partially occluded objects. To the best of our knowledge, a fully automatic method that can detect similar arbitrary objects from a single natural image has not been proposed. In the object detection stage, our system allows users to select a sample of objects that have multiple instances in the input image. Those objects may be subject to deformation, overlap, illumination influence, and appearance variation. Shape and color templates of user-selected object type(s) are extracted automatically. Our system then uses template matching to obtain two sets of candidate object instances by separately using the shape and color descriptor. A joint optimization scheme is developed to decide the final locating result.

In boundary extraction stage, we use alpha matting to obtain a precise object region. Thus, in order to get a good matting result, we introduce a trimap foreground extraction algorithm to obtain better foreground.

Similar object locating

We use PaintSelection [START_REF] Liu | Paint selection[END_REF] on the input image to extract an instance from the similar objects. In terms of the selection of a feasible feature to use in the template matching process, different types of features may achieve dramatically different matching results (Fig. 3). Usually, building a stable template matching system with a single feature for similar object locating is difficult because of the objects' vivid inner textures, self-occlusion, or color/lighting variations. Therefore, in order to make the object location more robust and accurate, we use two features in the template matching process and separately obtain candidate object locations for further optimization.

Object feature descriptors

For its desirable performance in a variety of tasks, its speed, robustness, adaptability to sliding window search, and popularity in the community, HOG template descriptor [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] is integrated in the proposed method as a shape descriptor to detect the possible similar object locations. For each image window, the input image is divided into small spatial regions ("cells"), and for each cell N -bins local 1-D HOG is accumulated. In our experiments, we use N H b = 9 bins and set the cell size as 8 × 8. The orientation bins are spaced over 0 • to 180 • . Each pixel in one cell is voted to its histogram based on gradient direction and magnitude. However, in contrast with the method in [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF], a normalized local histogram is not used over a larger spatial region ("block") because normalization will decrease the separating capacity of the HOG descriptor in our template matching process. Figure 4(a) shows an example of HOG features extracted from an input image. To enhance the robustness of our algorithm in complex environments, a color feature descriptor is developed to perform a separate locating process. We first segment both the template and the input image by SLIC [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]. Subsequently, we use a three-channel color histogram to represent the distribution of colors in each super-pixel. In our experiments, all colors are represented in HSV color space, and each channel is divided into N C b = 32 bins. Figure 4(b) shows an image segmentation example using SLIC method.

Object locating by matching

Similar objects are accurately located using object descriptors. For template T, let T = {p m |m = 1, . . . , N T } denote the set of feature points (HOG or color) in the template, where N T is the number of elements in T. To find similar objects efficiently, the template in the input image at each possible location is matched by scaling and rotating the image according to pre-set discrete intervals. In all our experiments, we pre-compute the object template in the discrete space of seven discrete angles {-120 • , -45 • , -10 • , 0 • , 10 • , 45 • , 120 • } and five scales {0.4, 0.8, 1.0, 1.2, 1.6}. We mainly employ histogram intersection to calculate the matching scores at each candidate location by separately using HOG and color descriptor. Two score maps are used to record the HOG and color matching scores at each location. Moreover, the feature point set (HOG or color) of the current sub-window w is denoted as W = {q n |n = 1, . . . , N w }, where N w is the number of elements in w. Thus, for a rotated and scaled template T i at a candidate location j , we formulate equations to calculate the HOG matching score S H and the color matching score S C as:

S H (T i , w j ) = 1 N H T i p m ∈T H , q n ∈W H N H b k=1 min(T H k,p m , w H k,q n ) N H b k=1 w H k,q n (1) S C (T i , w j ) = 1 N C T i p m ∈T C , q n ∈W C N C b k=1 min(T C k,p m , w C k,q n ) N C b k=1 w C k,q n , ( 2
)
where T H k,p m and S H k,q n are the corresponding HOG histograms of points p m and q n , and T C k,p m and S C k,q n are the corresponding color histograms, respectively. Both scores are normalized to [0,1].

After scanning the entire input image with scale and rotation variations, we choose the maximum score at each possible location as the final score and then record this value in the score map. Then, in both score maps, each local maximum matching score that is larger than η = 0.5 is chosen as a candidate object location. We ignore a local maximum if another local maximum with a larger matching score is seen within the distance threshold d. In our experiments, we set d to half of the bounding box circumcircle radius of the template. We take P H = {X i |i = 1, . . . , N H } to denote the set of candidate object locations acquired from the HOG score map, and P C = {Y j |j = 1, . . . , N C } to denote the set acquired from the color score map. N H and N C are the numbers of elements in P H and P C .

As shown in parts (a) and (b) of Fig. 3, the locating result obtained by using the shape and color descriptors may be evidently different due to the complex visual variations of the objects. Different features are sensitive to different variations. Moreover, wrong locations may also be included in the candidate sets because of the disturbance of the background ((a) and (b) of Fig. 3). Therefore, to obtain accurate object locations, high-confidence locations must be determined first. For the candidates in the two sets, we separately calculate the Euclidean distance of each candidate pair (X i , Y j ), and use the following metric to include the feasible pairs (stable locations) into a new set P S , which is formulated as follows:

P S = X i , Y j |X i = arg min i D X i , Y j , Y j = arg min j D X i , Y j , D X i , Y j ≤ D max (3)
where D(X i , Y j ) is the Euclidean distance between X i and Y j , X i ∈ P H , Y j ∈ P C . We set D max to half of the bounding box circumcircle radius of the template. If N p denotes the number of stable locations in P S , we then define the stable confidence values of HOG and color descriptors as c H = N p N H and c C = N p N C , respectively. Thus, the coordinates of the stable locations are calculated as follows:

Z k = β • X i + (1.0 -β) • Y j , k = 1, . . . , N p , (4)
where β = c H c H +c C . As shown in Fig. 3(c), we use the coor- dinates of the stable locations as the final location of similar objects.

Result improvement by matting

In most cases in our experiments, all the similar objects can be found accurately by using the above locating method. However, a few objects may sometimes be lost during the screening process of using the metric in Eq. ( 3). Certain objects may appear to have apparent shape or color differences compared with the template (but still belong to the same category) because of some environmental or artificial conditions. Therefore, to enhance the effectiveness of our algorithm, we combine the S H and S C as S Combine , and then use KNN matting [START_REF] Chen | KNN matting[END_REF] to complete the result set. With P S , we can use the score map, which is more similar to the template, as the main component of S Combine and define it as follows:

α H = N p N H , α C = N p N C , ( 5
)
S Combine = α H * S H + α C * S C α H + α C ,
where N p is the number of stable object locations in P. Thus, we sample randomly around the locations of P as the foreground locations, and likewise, sample randomly where scores are lower than 0.1 as background locations. Using this information, we can obtain an alpha score map S alpha by KNN matting, which describes the confidence of each pixel. Finally, for each local peak score in S Combine , we choose the locations in which the local average alpha score is larger than 0.9 as the final similar object locations. As shown in Fig. 5, all the similar objects are located after the matting process.

Contour extraction

With the information on similar object locations, we can employ alpha matting [START_REF] Levin | A closed-form solution to natural image matting[END_REF] to obtain the precise region of each location and then extract its outer contour. In using alpha matting, we first need to construct a trimap for each object to specify the foreground, background, and unknown regions. However, alpha matting may not perform as desired if the small area around the object location is used as foreground, as such information is too limited to extract a precise object region. A new method (Algorithm 1) solves this problem by extracting the foreground area for an object trimap, as illustrated in Fig. 6. With the confidence score given by S alpha (calculated from the KNN matting operation in Sect. 4.3), we can estimate a foreground area around each object location. Thus, we treated the area determined by Algorithm 1 as the foreground, the area out of the template window as the background, and the area in the middle as the unknown region. However, if the background or the inner texture is complex, alpha matting will fail to extract precise object contours. S tmpl denotes the area covered by the template contour in one object location, and A matting , the area extracted by alpha matting. If

A matting -A tmpl A tmpl > 0.3,

we use an active

Algorithm 1: Trimap Foreground Extract Input: S alpha , template area A tmpl , candidate locations set P Output: Foreground area of each candidate location 1. Convert S alpha to a binaries map B alpha using threshold 0.8, then dilate B alpha using 3 × 3 template in order to fill small holes; 2. Calculate a distance map Map dist ; for each pixel p in B alpha do Map dist [p] = the Euclidean distance between p and the nearest zero pixel; 3. Build a contour tree of Map dist . Set the lower contour as the root and the higher contour as the leaf; 4. Let L area as candidate foreground area; for each local maximum point p in Map dist do a. Find the local maximum point q which the Euclidean distance between p and q is smaller than

A tmpl , put p into set Q; b. R is a contour; if Q == ∅ then
Let R as the largest contour which area(R) < A tmpl and p in the area of R;

else Let R = largest contour; for each point q in set Q do (1) Let F as the lowest common ancestor between p and q in contour tree;

(2) Let U as the largest contour which smaller than S tmpl and contain p;

(3) Let R = U ;

Add R to L area ; 5. For the candidate locations p, the foreground area is L area ∩ TemplateWindow(p); contour-based model in [START_REF] Cheng | RepFinder: finding approximately repeated scene elements for image editing[END_REF] to refine the object contour further. In most cases, the contour extracted by alpha matting is sufficient, and is usually better than that by pure active contour method for the final object segmentation. As shown in Fig. 7, the object contours extracted by our method are more precise and smooth than the ones by only using active contour method.

Results and discussion

Our interactive system is targeted toward high level image pattern analysis. In the entire workflow, the cutout of the template is the only manual operation. We have implemented our algorithm in C++ on a computer with Intel Core i7 CPU at 3.9 GHz, 8 GB RAM, and Geforce GTX 670. Our CUDA-accelerated template matching plus the searching of stable locations typically takes 0.6 seconds to process The two leftmost occluded pandas are extracted through active contour, whereas the others are extracted by alpha matting a 900 × 900 image under seven rotations and five scalings. For our locating improvement and automatic object contour refinement through matting, the computation time depends on the image size, but usually consumes 2.0 to 4.5 seconds, thereby guaranteeing that the entire workflow of the system runs at interactive rates.

To validate the effectiveness of our SimLocator, we have tested it on a variety of input images with similar scene elements. Then, we discuss the experiments used to evaluate the performance of SimLocator qualitatively and to draw a comparison with other state-of-the-art repetition detection methods. The key ingredients of SimLocator are the similarity measurements S H (Eq. ( 1)) and S C (Eq. ( 2)), and the stable location selection metric P S (Eq. ( 3)). In the experiments conducted in this paper, we use both shape and color information. More feature descriptors, such as Gabor texture and SIFT/SURF descriptor, can also be incorporated into our method.

Comparison

We compare SimLocator with the state-of-the-art interactive repetition detection methods RepFinder [START_REF] Cheng | RepFinder: finding approximately repeated scene elements for image editing[END_REF], Dual-Bound [START_REF] Schweitzer | A dual-bound algorithm for very fast and exact template matching[END_REF] and RepSnapping [START_REF] Huang | RepSnapping: efficient image cutout for repeated scene elements[END_REF]. To attain a fair comparison, we also run RepFinder on the original image without segmenting the entire foreground first, as in their paper. In fact, for many examples, segmenting the foreground by using a few strokes is difficult, and thus we treat this pre-processing requirement as a strong limitation for a similar objects detection system. By contrast, our algorithm can accurately and efficiently locate similar objects directly on the original image, which allows for more practical uses than RepFinder. Figure 8 shows the results of the comparison. Our method is more robust than the other template- RepFinder is sensitive to the defocused effect (dandelions), moderate shape changes (chocolate cakes and strawberries), and complex textural backgrounds (green lanterns and purple flowers). In our algorithm, however, the integration of the color feature and stable location searching scheme effectively solve those problems. On the other hand, Dual-Bound is somewhat too accurate for a matching method. In most cases in our experiments, Dual-Bound can only find the template itself, which is not suitable for locating similar objects in natural images. Finally, RepSnapping can simultaneously cut out almost all the repeated elements with limited user interactions. However, the method is very sensitive to color variations (green lanterns), and is also unable to acquire information on individual object locations. Figure 9 illustrates more comparison results. SimLocator is more robust than RepFinder when the image has a perspective effect (meatballs and strawberries), 3D transformation (flowers), or severe shape deformation (fried balls and strawberries).

Template selection

In our experiments, the selection of template does not largely affect the setting of parameters and the accuracy of the detection result, if the user does not choose an object that is too different from the other objects. As shown in Fig. 10, our algorithm is robust when using variant templates. In this respect, SimLocator outperforms RepFinder [START_REF] Cheng | RepFinder: finding approximately repeated scene elements for image editing[END_REF] for both locating operations.

Limitations

SimLocator provides an efficient computational framework to locate similar elements in an image. Our algorithm is robust to intra-class appearance variations, illumination changes, and complex background, as well as certain levels of deformations and occlusions. However, like all other locating approaches, our algorithm fails in some cases, especially in images with 3D transformations, severe occlusion, or dense clusters, which may lead to some errors, as shown in the two examples of failed applications in Fig. 11. Moreover, the object contours may not be precise because of the influence of background gradients. In these cases, we can integrate more feature descriptors in template matching or specify more strokes to classify the objects from the background manually.

Conclusion and future work

We have developed a novel framework called SimLocator for locating and extracting similar elements in an image scene with very little user interaction. Unlike previous methods that only considered single features in matching, we incorporate shape and color information in the computation of our dense correspondences across objects, thus making our method more robust than other state-of-the-art repeat/similar pattern analysis systems. In the future, we plan to test more feature descriptors in the matching process, for example using Gabor texture descriptors to help to locate objects with similar inner textures but different color patterns. Developing new applications that can benefit from the information of images with similar content is another direction. Moreover, automatic similar object detection and extraction remains a very challenging research direction that merit further investigation.

Fig. 1

 1 Fig. 1 Comparison of different methods. Results show that the Sim-Locator method outperforms RepFinder and Dual-Bound in both accuracy (detecting the correct object locations) and completeness (finding

Fig. 2 Fig. 3

 23 Fig. 2 Algorithm workflow. User scribbles indicate an object template. We use Histogram of Oriented Gradients (HOG) as the shape feature descriptor and SLIC (Simple Linear Iterative Clustering) to generate the color feature descriptor (histogram-based). Similar object

Fig. 4

 4 Fig. 4 Shape (HOG) and color (super-pixel histogram) features used in our template-matching algorithm

Fig. 5 Fig. 6

 56 Fig. 5 Optimizing the locating result by stable locations and KNN matting. Wrong locations can be removed in this process

Fig. 7

 7 Fig. 7 Object contour extraction. (a) Object contour is extracted only using active contour. (b) Object contour is extracted using our method. The two leftmost occluded pandas are extracted through active contour, whereas the others are extracted by alpha matting

Fig. 8

 8 Fig. 8 Comparison of different methods locating similar objects. Templates are shown in the leftmost column

Fig. 9

 9 Fig. 9 Comparison of our method with RepFinder. Templates are placed beside the image samples

Fig. 11

 11 Fig. 11 Examples of failed runs. The school of fish in the left image has a dense cluster and severe occlusion. One cup in the right image suffers from 3D transformation and has a semantically different appearance from the template and the other objects

Acknowledgements

We thank anonymous reviewers for their valuable input. We thank Fuzhang Wu for making some Dual-Bound results. We thank the flickr members who kindly share their images under Creative Commons License: Bestfriend (tea cups), RenateEurope (purple flowers), acaffery (balloons), Luko Gecko (fish school), Bold Huang (petunia), Swami Stream (chocolate cake). We thank the users of 500px.com who have shared their images through public domain: Ricky Marek (pomegranates), Nitin Prabhudesai (dandelions), Joram Huyben (green lanterns). We also thank the users of pinterest.com who have put the media in their spaces: Gus's Mom (panda cakes), Rhonda E. Peterson (meatballs), Cindy Loo (fried balls). The two fruit cake images are borrowed from [10]. This work is supported by National Natural Science Foundation of China under project Nos. 61172104, 61271430, 61201402 and 61202324, by Beijing Natural Science Foundation (Content Aware Image Synthesis and its Applications, No. 4112061), and by SRF for ROCS, SEM.


Marc Bonnet

Aliénor Burel
email: alienor.burel@inria.fr

Patrick Joly

Effective Transmission Conditions for Thin-Layer Transmission Problems in Elastodynamics

Introduction

This research is motivated by the numerical modeling of ultrasonic non-destructive testing experiments. Some tested media feature thin layers (made e.g. of resin), which are difficult to handle in numerical computations due to the very small element size required for meshing them. Similar issues arise in the treatment of thin coatings, see e.g. [START_REF] Haddar | Stability of thin layer approximation of electromagnetic waves scattering by linear and non linear coatings[END_REF], [START_REF] Bendali | Lemrabet The effect of a thin coating on the scattering of a time-harmonic wave for the Helmholtz equation[END_REF]. To overcome these difficulties, one idea consists in using effective transmission conditions (ETCs) across the two interfaces bounding the layer. This work aims at establishing such ETCs by means of a formal asymptotic analysis with respect to the (small) layer thickness, in the spirit of [START_REF] Chun | Highorder accurate thin layer approximations for time-domain electromagnetics, Part II: Transmission layers[END_REF] for Maxwell's equations.

Problem Setting

We consider the case of a thin layer of an isotropic elastic material occupying the strip

Ω i η = R×[-η 2 , η 2
]. In addition, let Ω + η (resp. Ω - η ) denote the remaining portions of the propagation domain situated above (resp. below) the thin layer. The upper and lower boundaries of the layer are denoted ∂Ω + η and ∂Ω - η ; they are connected to Ω + η and Ω - η , respectively. We assume that the layer material is isotropic and homogeneous, with mass density ρ i and Lamé coefficients λ i , µ i . We denote by ρ, λ, and µ the (possibly heterogeneous) coefficients of the material in Ω ± η .

-x

Ω + η Ω i η Ω - η ∂Ω + η ∂Ω - η 6 y 0 - η 2 + η 2 6 6 n n
The displacement field u i η in Ω i η , as well as the displacement fields u ± η inside Ω ± η satisfy the elastodynamics equations:

ρ i ∂ 2 t u i η -div σ i (u i η ) = 0, in Ω i η , (1)
ρ ∂ 2 t u ± η -div σ(u ± η ) = 0, in Ω ± η . ( 2
)
where σ(u) and σ i (u) respectively denote the stress tensor in the surrounding and layer material, respectively, given by Hooke's law applied to a given displacement u. Equations ( 1) and ( 2) are coupled with transmission conditions on the interfaces ∂Ω ± η (we omit the time variable t for simplicity):

     u ± η x, ± η 2 = u i η x, ± η 2 t(u ± η ) x, ± η 2 = t i (u i η ) x, ± η 2 , (3)
where t(u) := σ(u)n and t i (u) := σ i (u)n are the traction vectors relative to Ω ± η and Ω η and n is the normal vector to ∂Ω ± η (see the figure above).

Eliminating formally u i η , we can write a transmission problem for u η := (u + η , u - η ). For any function f : R 2 → R d , using the notation

{f } η = [f ] η , f η : R → R 2d , [f ] η (x) := f (x, η/2) -f (x, -η/2), f η (x) := f (x, η/2) + f (x, -η/2) /2,
this transmission condition can be written in the form

{t(u η )} η + T η {u η } η = 0
where T η is a (nonlocal) DtN transmission operator that can easily be defined implicitly from the solution of the interior Dirichlet problem in the strip Ω i η . The next idea is that, when η tends to 0, T η becomes local and that one can get explicit analytical approximations of it.

Principle of construction of ETCs

This construction is based on an ansatz for the interior solution u i η of the form

u i η (x, y) = U 0 x, y η + η U 1 x, y η + η 2 U 2 x, y η + ... (4) where U k : Ω i 1 → R 2 .
This implies in particular a similar expansion for the traces

u i η (x, ± η 2 ) = u 0 ± (x) + η u 1 ± (x) + η 2 u 2 ± (x) + ... (5)
Substituting (4) into (1) allows us to compute explicitly the U k from the u k ± by induction on k: these are polynomial functions in y/η. These expressions lead us to introduce a family of differential operators of order , A (∂ x , ∂ t ), ≥ 0, such that

A 0 (∂ x , ∂ t ){u 0 } = 0 (6)
and

t i (u i η ) x, ± η 2 = t 0 ± (x)+η t 1 ± (x)+η 2 t 2 ± (x)+... ( 7
)
where

{t k } = k+1 j=0 A k+1-j (∂ x , ∂ t ){u j } (8)
and where we have defined

{t k } = [t k ], t k , [t k ] = t k + -t k -, t k = t k + + t k - 2
and the same for {u k }. Note that from (5) and ( 7)

   {u i η } η = {u 0 } + η {u 1 } + η 2 {u 2 } + ... {t i (u i η )} η = {t 0 } + η {t 1 } + η 2 {t 2 } + ... (9)
We rewrite the transmission conditions (3) as

{t(u η )} η = {t i (u i η )} η , {u η } η = {u i η } η
so that, using (9) and (8), we get

η {t(u η )} η = k≥0 η k+1 k+1 j=0 A k+1-j (∂ x , ∂ t ){u j }
which, thanks to (6), can be rearranged as

η {t(u η )} η = η A (∂ x , ∂ t ) j η j {u j } = k =0 η A (∂ x , ∂ t ) {u η } η + O(η k+1 )
The transmission condition of order k + 1 is then obtained formally by dropping the O(η k+1 ) term.

A stable transmission condition

Applying the above method with k = 2 leads to the following transmission conditions

A[u η ] η = η t(u η ) η -η BJ ∂ x u η η , [t(u η )] η = η (ρ ∂ 2 t -JAJ ∂ 2 x ) u η η -JB ∂ x [u η ] η , (10)
where A, B and J are the following 2 × 2 matrices:

A = µ i 0 0 λ i +2µ i , B = µ i 0 0 λ i , J = 0 1 1 0 . (11)
A fundamental point is the well-posedness and uniform stability in η of the transmission problem (2), (10). This is in fact a consequence of an energy conservation result: any smooth enough solution of ( 2), (10) satisfies:

d dt E η + E i η = 0 (12)
where

E η = ρ 2 Ωη |∂ t u η | 2 dx + 1 2 Ωη σ(u η ) : ε(u η ) dx E i η = η ρ 2 R | ∂ t u η | 2 dγ + η 2 R Q ∂ x u η , [u η ] η dγ
and Q(x, y) : R 2 × R 2 → R is the symmetric quadratic form (J T = J and (BJ) T = JB) Q(x, y) := x T JAJ x + 2 x T JB y + y T A y (13)

The well-posedness and stability result is then a consequence of the Theorem. The quadratic form Q(x, y) is positive.

At the conference, we shall present various numerical simulations to illustrate the accuracy and the efficiency of our approximate model. Moreover, some insights about the error analysis will be given.


Nicolas Barnier
email: barnier@recherche.enac.fr

Pascal Brisset
email: brisset@recherche.enac.fr

Thomas Rivière
email: riviere@recherche.enac.fr

Slot Allocation with Constraint Programming: Models and Results

Keywords: Slot Allocation, Constraint Programming, Combinatorial Optimization, Modelling

Current European Air Traffic Control system is far exceeded by the demand and the resulting delays are a financial and psychological burden for airlines and passengers. The Central Flow Management Unit, in charge of regulating the flights to respect en-route capacity constraints of Air Traffic Control Centres, uses a greedy algorithm to allocate departure slots which features several drawbacks concerning soundness, interpretation of the constraints and optimization. A new model taking advantage of Constraint Programming (CP) has been proposed within the SHAMAN system (from CENA) but still suffers from unevenly distributed workload and capacity violation. This paper presents two alternative models, also implemented using CP technology, of the slot allocation problem focused on the controllers workload: an extension of the SHAMAN model with a standard formulation, and a novel approach involving the sort constraint. Both of them can maintain workload constantly below a given capacity and the latter also provides efficient failure proof on over-constrained instances. The behaviours of the different models are discussed and results are presented with partial and full instances from real French air traffic data set. We eventually describe the potential operational improvement supplied by these continuous models.

Introduction

Airspace congestion is today the most critical issue European Air Traffic Management (ATM) has to face. French Air Traffic Control Centres (ATCC) capacities are far exceeded by a constant growth in air traffic demand, resulting in ever increasing flight delays. These time and management costs are such a nuisance for all airlines and passengers that the European Commission has released a special statement [IP/99/924, 1999] acknowledging that current ATFM systems are unable to support high traffic loads and unscalable for the predicted growth.

The Central Flow Management Unit (CFMU) in Brussels is in charge of reducing these congestion costs by, among several other strategic or tactical measures, delaying departure slots for the flights involved in overloaded sectors. The purpose of delaying is to respect the en-route capacity constraints provided by each ATCC as a number of planes per hour according to their daily schedule.

The CFMU currently solves this problem in two stages [START_REF] Cfmu | Basic CFMU Handbook -General & CFMU Systems[END_REF]:

A "pre-tactical" tool (PREDICT) identifies overloaded sectors and helps (an experienced human operator) to place regulations by simulating their impact.

Then a "tactical" system (TACT) including the CASA (Computer Assisted Slot Allocation) tool uses a greedy algorithm to dynamically allocate regulated departure slots on a "first planned first served" basis while flight plans data are received.

Basically, CASA computes a list of slots by dividing the time length of the regulated period by its capacity and tries to fill each slot according to the estimated arrival time of the flights in the concerned sector. Several revision processes occur when new flight plans are received that may shift all the flights already scheduled during the same period to respect the "first planned first served" principle.

The solution obtained obviously depends very much on the resolution order of the constraints and no optimization is possible 1 . Moreover, the regulation for a given flight is computed with respect to the sector inducing the greatest delay (in case of multiple regulations), which may lead to capacity constraints violation for the other sectors crossed by the same flight. Eventually, a former study [START_REF] Gotteland | Building and evaluating a minimal regulation scheme[END_REF] estimates that CFMU local regulations have little effects on ATCCs overall overload. However, CASA is an operational tool which must comply with other additional constraints (not taken into account in the models presented in this paper) and well integrated in a complex ATM process.

To overcome the drawbacks and lack of efficiency of CASA, the CENA (French Air Navigation Study Centre) has developped a slot allocation module based on Constraint Programming within the SHAMAN 2 platform [START_REF] Plusquellec | Description du module d'allocation de crneaux utilisant la programmation par contraintes implant dans shaman[END_REF]]. This tool states capacity constraints on all of the open sectors by dividing their overall lengths in contiguous 30 min long slices and constraining the number of flights arriving in the sectors over each slice (e.g. suppose sector AIX is open from 12h00 until 13h30 with capacity 32 flights per hour, three constraints ensuring that no more than 16 flights will enter the sector between respectively 12h00 and 12h30, 12h30 and 13h00, 13h00 and 13h30, are posted). In this approach, the use of Constraint technology provides a global view of all the capacity constraints such that solutions (if found) are consistent (all the constraints are satisfied).

However, the SHAMAN model does not prevent traffic peaks exceeding the sectors capacities during periods which do not begin at multiples of 30 min from the start time of their constraints, such that controllers have to face too high workloads. Furthermore, even if SHAMAN does not perform any optimization on its solution 3 , the search strategy is guided by the minimization of the sum of the delays. Hence delayed flights have a chronic tendancy to gather at the beginning of the 30 min periods, yielding traffic peaks above the requested capacities.

As the statements of capacity constraints are obviously not very satisfactorily interpreted by current ATFM systems from an operational point of view, a better model of these requirements would be to maintain controllers workload below the specified capacities continuously. We provide for this aim several new models endowed with a more realistic, thus harder, satisfaction of capacity constraints:

Sliding discrete windows which allow to smooth the profile of controllers workload from the weakest constrained model (the SHAMAN model) to the hardest (the satisfaction of capacity constraints over any time period of a given length);

The "sorting" model which states continuously the constraints over the ranks of the entry time of flights within a given sector.

Constraint Programming allows to express a very straightforward formulation of these two latter models, as well as for a "twin" of the SHAMAN model presented for comparison purpose. The versatility of CP technology eases fast implementation and alternatives testing. An all-purpose constraint library, FaCiLe 4 [START_REF] Barnier | FaCiLe : A Functional Constraint Library[END_REF] written in Objective Caml [Leroy, 2000] has been used to experiment with the various models, leading to drastically different allocation schemes and abilities to obtain proof of optimality or unfeasibility.

The paper is divided as follows: we first give a precise description of the problem and some hints about the size and complexity of the instances we are interested in. The next section introduces briefly Constraint Programming technology, thereafter we present several models for the capacity constraints, starting with standard ones and refining them with a more continuous formulation. Eventually, section 5.2 is devoted to the results for simplified and real instances where behaviour of the different models are compared. We conclude with an overview and some indications about possible future work.

Slot Allocation Problem Description

Air traffic flow management is a daily pre-tactical filter intended to regulate scheduled flights across controlled airspace. Its main aim is to limit the number of planes in a given space over a given period. This planning will be precisely scheduled by the airtraffic controllers in real-time. The ATFM problem is described in terms of Flights: A flight starts from one airport at a specified time, follows a predefined route at a fixed speed and arrives at another airport.

1 Except at the very final stage when the slots allocated to lately cancelled flight may be redistributed. The partition of the airspace, i.e. the number and the shape of the sectors, varies during the day. The capacity of the same sector may also changes at given times. We call sector-period a sector for particular period of time and capacity.

The constraints of the problem are the capacities of the sectors. There are several degrees of freedom to satisfy these constraints: choosing different routes for flights, delaying departures, changing speed of the aircrafts during flight, asking aircrafts to hold their position... In this paper, we focus only on the ground-holding policy: each flight may be delayed at departure. ATFM is in charge of solving the problem on a day-by-day basis.

The difficulty of the problem does not come from the complexity of its constraints but from its size. We have used real data archived by the French civil aviation tool COURAGE to solve the ATFM problem. E.g. for May 20th, 1999:

The French airspace is concerned by 7375 flights entering between 0h00 and 23h59. 140 sectors are active during the day, some of them with various capacities (up to 6 different).

More than 700 flights may enter a single sector-period.

Capacities vary from 19 to 52 flights per hour.

The objective of the slot allocation problem is to reduce the delays. There are several ways to achieve this purpose: reducing the total sum of the delays (utility), the maximum delay (equity), the average delay, etc. [Maugis, 1996] gives an extensive description of what could and should be a cost function and finally concludes with results for the simplest one, the total sum of delays. It is also the choice of [START_REF] Bertsimas | The air traffic flow management problem with enroute capacities[END_REF] with a slightly different model. In this study, we choose to minimize the maximum of all delays, but we are more interested in the qualitative properties of the solution than in its numerical cost.

Constraint Programming

Constraint Programming (CP) is an emergent modelling technology for declarative description and effective solving of large combinatorial optimization problems. The idea is to solve problems by stating constraints (conditions, properties on the decision variables) which must be satisfied by all the solutions while avoiding the tedious and error-prone task of explicitely maintaining the consistency of these constraints as search goes on.

Constraint solvers efficient enough to compare well with Operational Research technology emerged in the early 90s, based on the CLP (Constraint Logic Programming) paradigm and integrated in Prolog systems [Van Hentenryck, 1989]. Nowadays, CP is a challenging technology addressing a wide range of combinatorial academic and industrial hard problems in many fields (planning, scheduling, crew rostering, configuration, circuit design, DNA sequencing etc). State-of-the-art constraints solvers like Ilog Solver [START_REF] Solver | ILOG Solver 4.4 user's manual[END_REF], Eclipse [START_REF] Ecl | ECL ✔ PS⑨ user manual[END_REF], CHIP [START_REF] Dincbas | The constraint logic programming language CHIP[END_REF], Prolog IV, SICStus Prolog or Mozart are used by more and more important companies (British Airways, British Telecom, SNCF, Hong Kong International Terminals, Michelin etc [Barták, 1998]).

CP is a very active research field and solvers technology evolves very quickly. One of the most promissing areas is the integration in constraints solvers of other paradigms of combinatorial optimization technologies, e.g. local search, stochastic search, linear programming. A current working version of a constraint program solving the slot allocation problem written with FaCiLe uses a generic invariants library [START_REF] Michel | Localizer: A modeling language for local search[END_REF] typically designed towards the implementation of local search algorithms [Rivière, 2001] 5 . The next three sections attempt to describe the main underlying ideas of CP.

Constraints

A constraint is a logical relation among several variables (or unknowns), each taking a value in a given domain (e.g. an integer interval). A constraint thus restricts the possible values that variables can take by maintaining its consistency. Basically, constraints remove from domains values which cannot appear in solutions: for instance an inequality constraint ✁ ✂ ✄ between two variables ✁ and ✄ with respective domains ☎ ✆ ✝✝✞✟✠ and ☎ ✟ ✝✝✞✡✠ will reduce the variables domains to respectively ☎ ✆ ✝✝✞☛✠ and ☎ ☞ ✝✝✞✡✠ .

Each time a constraint is posted to the constraint store of a solver, the involved variables domains are refined and trigger the awakening of other constraints involved in one of the formerly reduced variables. This process is repeated until a fix point is reached where no more reduction can be deduced. Constraints are declarative statements, close to natural description of combinatorial problems, which describe what relationship must hold without specifying a computational procedure to enforce it. This feature confers to CP solvers the status of modelling tools. They allow to express straightforwardly practical problem thanks to a rich language of constraints from linear or non-linear arithmetic to global constraints on sets of variables (difference, cardinality, sequence etc) and higher-order constraints (or logical constraints on constraints6 ). Modern solvers integrates a wide variety of builtin constraints whose consistency schemes are borrowed from state-of-the-art OR techniques, thus offering the versatility of CP modelling features while preserving the efficiency of specialized and wellproven algorithms.

Search

Posting constraints alone does not reduce domains enough to find a solution (at least in non-trivial problems or when the problem has several solutions) and search has to be performed. Basically, a variable is chosen as well as an instantiation value in its current domain, then related constraints are waken and the propagation process starts, reducing the domains to prune the search tree as much as possible. Thereafter another variable is chosen and instantiated until all variables are bound to a value and a solution is found or a failure occurs. In the latter case, the search algorithm backtracks to reconsider the last choice point. An optimization process can be elegantly integrated to this branching mechanism by dynamically adding a constraint that enforces the cost of the problem to be strictly less than the one of the last solution found.

CP languages are very well suited to express readily complex search strategies based on dynamic ordering (i.e. which depends on the state of the system) of variables and instantiation values. Such nondeterministic search procedures are implemented with Prolog-like goals. E.g. one of the most popular strategy, called the "first fail principle"7 , is to choose the uninstantiated variable with the smallest domain and/or involved in the greatest number of constraints. Most constraint solvers provide built-in goals to implement such heuristics.

The orthogonal strengths of a rich declarative constraint language and ability to specify search strategies with a high level of abstraction help to separate the problem modelling from the search for solutions. Hence they improve the correction of programs and allow fast prototyping, developpement and alternatives testing.

Formal Description

A Constraint Satisfaction Problem (CSP) consists of:

A set of variables ✌ ✍ ✎✁ ✏ ✑ ✝✝✝ ✑ ✁ ✒ ✓.
For each variable ✁ ✔, a finite set ✕ ✔8 of possible values, called its domain; we note

✕ ✍ ✕ ✏ ✖ ✝✝✝ ✖ ✕ ✒ . An assignement ✗ ✍ ✘✙ ✏ ✑ ✝✝✝ ✑ ✙ ✒ ✚ of ✌ is a valuation of each vari- able of ✌ from a value of its domain: ✗ ✛ ✕ . A set of constraints ✜ ✍ ✎✢ ✏ ✑ ✝✝✝ ✑ ✢✣ ✓.
Each constraint ✢✔ restricts the values that a subset of variables ✎✁ ✔✤ ✑ ✝✝ ✑ ✁ ✔✥ ✓ ✦ ✌ can simultaneously take. Therefore, constraint ✢✔ can be represented as a subset ✕ ✧★ of the cartesian product ✕ ✔✤ ✖ ✝✝✝ ✖ ✕ ✔✥ and a function:

✢✔ ✩ ✕ ✪ ✫ ✎true✑ false✓ such that ✬ ✗ ✍ ✘✙ ✏ ✑ ✝✝✝ ✑ ✙ ✒ ✚ ✛ ✕ ✑ ✘✙ ✔✤ ✑ ✝✝ ✑ ✙ ✔✥ ✚ ✛ ✕ ✧★ ✭✮ ✢✔ ✘ ✗ ✚ A solution to a CSP is an assignment ✗ of ✌ such that every constraint is satisfied, i.e. ✯ ✣ ✔✰ ✏ ✢✔ ✘ ✗ ✚.
Depending on the problem, to solve a CSP means to find: a proof of failure, i.e. there is no solution;

just one solution, with no preference as to which one; all solutions; an optimal, or at least a good solution, given some objective function defined in terms of some or all of the variables.

Actually, real-life problems (huge size and/or high complexity) sometimes cannot yet be solved optimaly. In such cases, "good solutions" obtained with a limited amount of optimization (e.g. with a computation time bound) or satisfying a subset of the constraints are searched for.

Three Models

In this section, we describe three models for the slot allocation problem. These models are not equivalent and differ according to the interpretation of the sector load constraints.

The models are described with the following data:

✾ : the set of sector-periods, each with a ✿❀❁❂ ❀ and an ❃❄ ❅; ❆ : the set of flights; We first present the Basic model implemented within SHAMAN in section 4.1 [Maugis, 1996, Plusquellec and[START_REF] Plusquellec | [END_REF]]. Then two novel models are introduced in section 4.2.

Non Overlapping Windows

The SHAMAN model (refered afterward as the Basic model) discretizes the sectors with ❊ ✍ •❍ min periods to enforce the capacity constraints and attempts to reduce the sum of all delays. Several unsuitable side-effects arise from this model, leading to irregular workload profiles.

The Basic model

We consider the load constraint in successive contiguous periods, i.e. non overlapping time windows:

■ ❇ ✍ ✎❉ ❇ ❏ ✑ ❉ ❇ ✏ ✑ ✝✝✝ ✓:
successive periods of length ❊ , each with a ✿❀❁❂ ❀ and a ❃❄ ❅. The first period starts with the sector-period:

✿❀❁❂ ❀ ✘ ❉ ❇ ❏ ✚ ✍ ✿❀❁❂ ❀ ✘ ✿✚.
In an other way:

■ ❇ ✍ ✎ ☎✿❀❁❂ ❀ ✘ ✿✚ ✑ ✿❀❁❂ ❀ ✘ ✿✚ ❑ ❊ ☎✑ ☎✿❀❁❂ ❀ ✘ ✿✚ ❑ ❊ ✑ ✿❀❁❂ ❀ ✘ ✿✚ ❑ ▲ ❊ ☎✑ ✝✝✝ ✓
The Basic model is written with auxiliary boolean variables:

▼ ✔ ◆ ❖ ◗ : flight ❈ enters the sector ✿ during the period ❉ ❇ | .
A first constraint relates the delay variables with the boolean variables. A second one sets the sector load capacity.

✬ ✿ ✛ ✾ ✬ ❉ ❇ | ✛ ■ ❇: ▼ ✔ ◆ ❖ ◗ iff ✿❀❁❂ ❀ ✘ ❉ ❇ | ✚ ❙ ❀❇ ✔ ❑ ❋ ✔ ✂ ✿❀❁❂ ❀ ✘ ❉ ❇ | ✚ ❑ ❊ ❚ ✔❯❱ ▼ ✔ ◆ ❖ ◗ ❙ ✢❁❉ ❁ ❇
This model needs a huge amount of boolean variables:

❲ ✎ ▼ ✔ ◆ ❖ ◗ ✓ ❲ ✍ ❲ ❆ ❲ ❳ ❇ ❯( ❃❄ ❅ ✘ ✿✚ ✪ ✿❀❁❂ ❀ ✘ ✿✚ ❊

Drawbacks

The Basic model suffers from strong discontinuity. It only states that the sector load limit must be respected at the beginning of each period, usually 48 times a day (for 30 min windows). Such that workload may climb much higher for time windows that do not start at multiples of ❊ min from the beginning of the constraint, as illustrated in figure 1 (with ❊ ✍ • ❍ min): here, for an overall amount of flights that does not exceed the capacity either in the first period or in the second one, an in-between period of the same length may have almost twice as much flights. Furthermore, the constraints of this model does not prevent traffic peaks for time windows shorter than the ❊ -length periods, as shown in figure 2: the two workload profiles correspond (roughly) to the same amount of flights (supposedly below the capacity for the duration of the period), though the first one (allowed by the Basic model) is obviously not well suited to ease the task of controllers. Because the main objective of the cost function is to reduce the delays, the expected side-effect of this model is a concentration of flights entering at the beginning of the periods: postponed flights are scheduled as early as possible in the next time window with the Basic model. Figure 3 illustrates this behaviour: the dashed curve is the workload pro-file before any regulation, then the Basic model scheme tends to shift the traffic peak of the first time window at the beginning of the second one. This chronic side-effect is confirmed by the experiments (see section 5.2.1).

Continuous Models

We propose two new formulations which attempt to smooth the workload with a more constrained modelling of capacity constraints. The first model is an extension of the Basic one. The second one relies on a sorting global constraint.

Sliding Periods

Keeping the same idea of the Basic model, it is possible to get a more continuous model with overlapping periods, illustrated in figure 4. Here we use another parameter ( ) ) which is the time step between periods:

■ ❇ ✍ ✎ ☎✿ ❀ ❁ ❂ ❀ ✘ ✿ ✚ ✑ ✿ ❀ ❁ ❂ ❀ ✘ ✿ ✚ ❑ ❊ ☎✑ ☎✿ ❀ ❁ ❂ ❀ ✘ ✿ ✚ ❑ ) ✑ ✿ ❀ ❁ ❂ ❀ ✘ ✿ ✚ ❑ ) ❑ ❊ ☎✑ ✝ ✝✝ ✓
Of course this model supersedes the Basic one and is equivalent for ) ✍ ❊ . It can be implemented with the same kind of boolean variables.

The workload profile is expected to become smoother as ) decreases, such that more ❊ -length time periods are taken into account.

The Sorting Model

The Sorting Constraint Let ❋ be a totally ordered set. The sort constraint is the relation associated with the standard sort function. A sort function takes a sequence of length ❄ containing elements in ❋ and returns another sequence containing the same elements ordered. The sort constraint "relates" two sequences containing finite domain variables taking values in ❋ .

For example, let

❬ ✍ ✎ ☎❍ ✪ ✞ • ✠ ❭ ☎ ✟ ✪ ✞ ❍ ✠ ❭ ☎ ✞ ❍ ✪ ✞✞✠ ❭ ☎ ☛ ✪ ✞✟✠ ❭ ☎ ☛ ✪ ✟✠ ✓
and

( ✍ ✎ ☎ ✞ ✪ • ✠ ❭ ☎ ✡ ✪ ✞ ❍ ✠ ❭ ☎ ✟ ✪ ☞✠ ❭ ☎ ✞✞ ✪ ✞)✠ ❭ ☎ ✞ ❍ ✪ ✞✡✠ ✓
be two sequences of interval variables. Posting the constraint

✿❴❂ ❀ ✘❬ ✑ ( ✚
must lead to the following refinements: [START_REF] Guernalec | [END_REF] proposed an efficient filtering algorithm for this constraint. Surprisingly enough, the complexity of this complete (i.e. the propagation refines the domains as much as possible) narrowing algorithm has an optimal complexity of ❵ ✘ ❄ ❛❜ ❝ ❄ ✚, the same than classic sort algorithms. However, the narrowing is done only on the bounds of the variables such that no propagation is performed if "holes" appear in the domains.

❬ ✍ ✎ ☎ ✞ ✪ • ✠ ❭ ☎ ✟ ✪ ☞✠ ❭ ✞ ✞❭ ☎ ✞✞ ✪ ✞✡✠ ❭ ☎ ✡ ✪ ✟✠ ✓ ( ✍ ✎ ☎ ✞ ✪ • ✠ ❭ ☎ ✡ ✪ ✟✠ ❭ ☎ ✟ ✪ ☞✠ ❭ ✞✞ ❭ ☎ ✞✞ ✪ ✞✡✠ ✓ [Guernalec

A Dual Model

For our problem, we use one sorting constraint for each sector-period. The constraint is set on auxiliary variables corresponding to the entry times and the sorted entry times of the flights: The constraints relate the auxiliary variables to the main delay variables and set the sector load limit with a dual model of the standard one. We no longer count the number of flights within a given period, but we rather consider the rank of the entry times and specify that two flights must be separated by a given duration whenever their ranks are too close.

Two flights with entry times ❡ ✔ and ❡ | in the ordered sequence such that ❈ ❑ ✢❁❉ ❁ ❙ ❢ must be separated by at least ❊ minutes from each other (c.f. figure 5): The previous constraint for the load limit is not totally correct because it is set even if one of the two concerned flights enters the sector before or after the corresponding period. A solution is to relax the constraint using auxiliary boolean variables:

✬ ✿ ✬ ❈ ✛ ❆ ❞ ❇ ✔ ✍ ❀❇ ✔ ❑ ❋ ✔ ✬ ✿ ✿❴❂ ❀ ✘ ❞ ❇ ✑ ❡ ❇ ✚ ✬ ✿ ✬ ❢ ✛ ✎ ✞ ✑ ✝✝ ✑ ❲ ❆ ❲ ✪ ✢❁❉ ❁ ❇ ✓ ❡ ❇ | ❑ ❊ ❙ ❡ ❇ | ❣ ✧❤ ❖ ❤ 30th 31st 10th 11th ✐ •❦ • ♠♥ ♦ ✐ ♣❦ ✐ •• • ♠♥ ♦ ✐ ♣ •
✬ ✿ ✛ ✜ ✬ ❢ ▼ ❇ | iff ✿❀❁❂ ❀ ✘ ✿✚ ❙ ❡ ❇ | ✂ ❃❄ ❅ ✘ ✿✚ ✬ ✿ ✛ ✜ ✬ ❢ ✛ ✎ ✞ ✑ ✝✝✝ ✑ ❲ ❆ ❲ ✪ ✢❁❉ ❁ ❇ ✓ ▼ ❇ | q ▼ ❇ | ❣ ✧❤ ❖ ❤ ✍✮ ❡ ❇ | ❑ ❊ ❙ ❡ ❇ | ❣ ✧❤ ❖ ❤
The boolean variable ▼ ❇ | states that the sector-period ✿ is concerned by the flight ❢ . If one of the two flights is out of the period, the constraint is not posted.

The Sorting model states the capacity constraints continuously over a given sector and thus is expected to maintain the workload constantly below the fixed limit. The Sliding model should be equivalent with ) reduced down to the time unit, i.e. all possible ❊ -length time windows are constrained.

Experiments

Implementation

The three models have been implemented with FaCiLe, a finite domain library we have written using the Objective Caml system [Leroy, 2000]. This strongly typed functional language provides well documented libraries for various data-structures. A fast compiler, ported on various OS and processors, produces efficient native code.

Our constraint library includes standard integer finite domain variables, linear and non-linear arithmetic constraints, reification, higher-oder constraints and some global constraints: difference, cardinality, sorting, etc. The search is controlled with goals in a Prolog way.

Results

We give in this section some results obtained for the different models applied on real data, for May 20th, 1999.

All the experiments are done with a precision r (time unit) of 5 min and a maximum delay of 60 min, except where mentioned otherwise.

A Single Sector

In order to analyse the solutions of the different models, we first looked at a very simplified problem with a single sector and no variation of capacity. We chose the sector (of May 20th, 1999) concerned by the maximum number of flights (644). The hourly capacity is ☛ ❍. The flights are expected to enter in the sector between 0h52 and 23h39.

It was first possible to check the equivalence and order between models (M1 ❙ M2 means that M2 is more constrained than M1, i.e. that a solution for M2 is a solution for M1):

Basic ❙ Sorting Basic s Sliding with ) ✍ ❊ Sliding ❙ Sorting Sliding with ) ✍ r s Sorting Table 1 reveals the corresponding numerical results: a more constrained model gets a higher cost. For various models and parameters settings are presented respectively in the last three rows the sum of all delays, the number of on-time flights (no delay) and the number of flights with delays less than 15 min. We see in this table that most of the flights are not delayed flights or get a small delay (compared to 644 concerned flights). Note that the average delay for all the solutions is smaller than the time precision (5 min). Numerical cost is a poor indicator about a solution. The expected differences between the three models are qualitative. In figures 6 and 7 is plotted the instantaneous load of the sector, i.e. the number of flights which will arrive during the next ❊ period. Dots correspond to the unregulated flow (the initial data) while solid curves show the solutions.

Model

❊ ) t ✉ ✈ ✇①✉ ✈ ② ③ ④ ✇ ✇①✉ ✈ ⑤ ⑥⑦ ④ ✇
At first glance the Basic model does not seem to regulate much. Actually, it only insures that the curve goes under the capacity bound every ❊ minutes (every hour starting from 0). The figure 7 corresponding to the Sorting model shows more interesting and expected results. The solution gives an almost constant load from 5h00 to 20h00.

Figure 8 shows the influence of the ) parameter on the Sliding model. With ) ✍ ❊ ✍ ✟ ❍, we get the solution of the Basic model. With a smaller ) (15), we observe that the load goes under the capacity every ) minutes but has still chances to go up 10% over the capacity in between.

Figure 9 shows the side-effects of the Basic model. In this experiment, we reduced the capacity (10%) and augmented the max delay (120 min) to force more flights to be delayed. In figure 9 is plotted the instant number of flights inside the sector (one can notice that a plane does not stay for a long time in a sector). Peaks occur at the beginning of this periods: a delayed flight has to be scheduled during the next non-full period; then a minimum delay corresponds to the beginning of the period. The Sorting model does not suffer from this side-effect and ensures an even load. 2 gives some indication of the ability of the models to prove that a problem has no solution. With a reduced capacity of 36, the failure could not be proved in a "finite time" (⑧ in the table) except using the Sorting model. Note the problem seems to have a phase transition since a solution is easily found with a slightly greater capacity (37).

Model

Full Problem

We were able to prove that there is no solution for the Sorting model with the given capacities. But solutions can be obtained with allowing an overload of the capacity of each sector. We found a solution with a capacity overload of 25%.

The regulation has no spectacular qualitative effect like in the simplified case; the whole problem is too much complex to have simple local properties. However, in figure 10 are compared the effects of the Basic and Sorting model on sector J around noon. The solid curves correspond to the flights entering the sector during the next 60 min and the dotted ones to the instantaneous number of flights present in the sector. This sector is very crowded and the regulation produced by the Basic model has little effect on the overload, whereas the Sorting model yields a load that levels out along the capacity line. This confirms the expected regular behaviour of the Sorting model.

Conclusion

Time slot allocation in ATFM is a hard combinatorial problem yet hardly solved by current CFMU systems. Capacity constraints are subject to interpretations and Constraint Programming provides an efficient and straightforward way to implement various models addressing control workload with different operational points of view. Non overlapping windows hardly cope with a realistic semantics of capacity constraints and would probably lead Control Centres to underestimate their capacities because of periodic traffic peaks whereas r-sliding windows and sort models ensure an even regulation inducing a much more constrained CSP with possibly higher overall delay costs. The sort model however needs fewer control parameters than its time windows counterpart. It is also far more efficient on fail proofs because of very powerful constraint propagation while offering similar computation time performances.

Large size and high dimensionality of air traffic input data make the slot allocation problem hard to optimize with respect to the cumulated delays criterion and future work should address this issue more efficiently. But the integration in the objective function of other factors (which may also be stated at the constraint modelling stage) like the regularity of workload distribution seems to be an important operational requirement and would penalize the basic models.

  ❀❇ ✔ : the time the flight ❈ enters the sector ✿ according to the original timetable; ✢❁❉ ❁ ❇: the capacity of the sector-period ✿ (flights/hour); ❊ : time base for the capacity constraint (minutes).All the models use the main decision variables ❋ ✔: departure delay for flight ❈.

Figure 1 :

 1 Figure 1: Overloaded in-between period with the Basic model.

Figure 2 :

 2 Figure 2: The Basic model does not prevent the first workload profile: peaks may arise for short time periods.

Figure 3 :

 3 Figure 3: Traffic peaks at the beginning of the period with the Basic model.

  Figure 4: Sliding windows

  time of the ❢ th flight entering the sector-period ✿.

Figure 5 :

 5 Figure 5: Capacity constraint applied to sorted flights ( ❊ ✍ • ❍, ✢❁❉ ❁ ✍ ▲❍).

Figure 6 :

 6 Figure 6: Regulation with the Basic model (flights entering during the next ❊ min). After regulation (solid curve), the load falls below the capacity every ❊ ✘ ✍ ✟ ❍ ✚ min.

Figure 7 :Figure 8 :

 78 Figure 7: Regulation with the Sorting model (flights entering during the next ❊ min). The Sorting model ensures that the load remains constantly below the capacity.

Figure 9 :Figure 10 :

 910 Figure 9: Instantaneous load (number of flights in the sector) for the Basic and the Sorting models. The Basic model tends to yield traffic peaks at the beginning of each period, whereas the Sorting model keeps a low average load.

  2 System to Help Analysis and Monitoring of Acc resources and Air route Network 3 Usually done in Constraint programs with a standard Branch & Bound algorithm.

	Sectors: The airspace is divided into sectors
	crossed by the routes of the flights. A sec-
	tor is a 3D polyhedra, usually a vertical cylin-
	der endowed with a capacity expressed as the
	maximum number of flights entering the sec-
	tor during a time period (usually one hour).
	4 Functional Constraint Library
	2

Table 1 :

 1 Delays for different models and parameters

	Basic	60	690	602	624
	Basic	30	1960	532	595
	Sliding 60 30 1760	549	597
	Sliding 60 15 2480	504	565
	Sorting, Sliding	60 r 3660	467	553

Table 2 :

 2 Proof of infeasibilityTable

		Capacity	Result
	Basic	36	⑧
	Sliding r	36	⑧
	Sorting	36	Failure proved
	Sorting	37	19775

In our application, it helps to maintain efficiently complex data structures involved in a dynamic search heuristic.

Like implication e.g. ✱ ✲ ✳ ✴ ✳ ✵ ✶ ✷ or disjunction e.g. ✱ ✸ ✹✺ ✻ ✳ ✼ ✳ ✸ ✹✽ ✻ ✱.

"To succeed, try first where you are most likely to fail."[Barták, 1998]

As formerly said, CP main contribution field is combinatorial problems, most of the time expressed with variables ranging over integer domains. However, some commercial and research solvers integrate efficient constraints over real, rational, finite set, boolean and even tree variables. This paper and our constraints librairy FaCiLe only consider integer domains.

Authors

Nicolas Barnier