Gildardo Lozano-Vega
email: gildardo.lozano@fh-mainz.de

Yannick Benezeth

Franck Marzani

Frank Boochs

Modular Method of Detection, Localization and Counting of Mutliple-Taxon Pollen Apertures Using Bag of Words

Keywords: pattern recognition, LBP, bag of words, apertures, palynology

Introduction

Allergic diseases due to airborne pollen affect considerably the world population. The correct estimation of the pollen concentration in the surrounding of individuals becomes important for the study of allergies. Particularly, the recognition of the taxa at genus level is relevant for the diagnosis, medication and prevention of allergic diseases. Traditional methods consist of counting manually the pollen particles under the microscope from airborne samples. These methods are time consuming and costly, limiting considerably the frequency and the extension of the region under analysis. Furthermore, this type of methods are highly susceptible to human error and inconsistency.

The recognition of pollen taxa is a complex task due to the high diversity of flowering plants. On one hand, pollen particles of the same taxonomic family usually share many similarities even when they belong to different genus branches (inter-class similarity). On the other hand, pollen are naturally formed particles whose characteristics present variations even in the same genus branch (intraclass variation).

Palynological literature describes morphological and aerobiological pollen characteristics, which are employed by specialists for the recognition of taxa [START_REF] Erdtman | An Introduction To Pollen Analysis[END_REF] [START_REF] Hesse | Pollen Terminology. An illustrated handbook[END_REF]. Average size, shape, equatorial outline, aperture type and amount, exine and intine structure, and type of ornamentation are among the main morphological characteristics, although not all can be precisely measured through common microscopes. Aerobiological characteristics refer to the period of pollination for each taxon and can be obtained by historic data from previous years. This information indicates the time slots when high levels of concentration of pollen are present in the environment.

The type and amount of apertures are particularly characteristic of the pollen taxa. Sometimes, they are the key factor to discriminate between two similar genus. Located at the pollen external wall, apertures are morphological distinctive regions which can be observed often with a common microscope due to its size. The appearance of the apertures varies according to the pollen taxon. Moreover, the orientation and focal plane of the particle with respect to the observer can also change their appearance considerably, and even make the aperture difficult to observe. Examples of allergenic pollen taxa are shown in Figure 1. In some of them, apertures are not visible at all. Common strategies for image-based classification of pollen employ global features like shape, texture and color; or local invariants [START_REF] Rodríguez-Damián | Automatic identification and classification of pollen of the urticaceae family[END_REF][4] [START_REF] Boucher | Development of a semi-automatic system for pollen recognition[END_REF][6] [START_REF] Ronneberger | 3D invariants with high robustness to local deformations for automated pollen recognition[END_REF]. Despite its importance, few works on detection of apertures exist. Boucher et al. developed specific algorithms according to the pollen type and to the point of view in order to detect the aperture, reticulum and cytoplasm of different pollen taxa [START_REF] Boucher | Development of a semi-automatic system for pollen recognition[END_REF]. Although no much detail is provided, different techniques were employed for each case such as thresholding, Laplacian of Gaussian, and region segmentation. The identified characteristic was then described by means of color, shape and size. Also oriented to the detection of aperture and colpus, Chen et al. [START_REF] Chen | Feasibility study on automated recognition of allergenic pollen: grass, birch and mugwort[END_REF] created different algorithms for Birch and Mugwort pollen. Hough transform and thresholding were used for detection of circular apertures. For apertures on the particle image border, they analyzed the intensity profile through the polar transform and applied a template matching strategy and frequency spectrum analysis to detect specific aperture types. Similar to Boucher et al., multiple algorithms were created for the diversity of pollen types, appearances and points of views. These approaches are inflexible if learning more different apertures is needed, requiring the development of new specific algorithms. In both cases, the aperture information was fed to a greater classification system for the recognition of taxa. Therefore, no individual evaluation of the aperture detection was shown, limiting the performance comparison among approaches.

We propose a different approach for the detection, localization and counting of apertures for different pollen types which overcomes the variability due to changes of appearance. The method relies on the Bag-of-Words (BOW) approach for the description of apertures in terms of primitive images. Then, multiple region images are sampled and classified to eventually create by voting, a likelihood map of the presence of apertures on the pollen image. The amount of apertures can be used eventually together with global features for the classification of pollen taxa. The proposed approach is able to learn new varieties of aperture by following a unique methodology, eliminating the need of designing new algorithms or methods. This flexibility is particularly important for growing the system with apertures with multiple appearances.

The rest of the paper is organized as follows: In section 2, the description of the apertures, detection, localization and counting methods are explained in detail. In section 3, experiments on region configuration and detection results are presented. Conclusions about the proposed approach and future work are given in section 4.

Recognition and Localization of Apertures

In this section, the methodology to represent apertures and to train the classification models is described. Additionally, the process of analysis of an unseen pollen particle for the detection of apertures is explained. First, a general overview of the method is presented in order to acquire the central idea of the proposal, and then the detailed procedure is described.

Outline of the proposed method

The proposed method starts describing the apertures by their content of primitive subimages. These primitive subimages, which form a visual codebook, are represented by a local descriptor and clustered following a similar procedure to BOW. Based on the codebook representation, a classification model learns to discriminate aperture regions form the rest of the particle. Following the same methodology, individual classifiers are created oriented to the different pollen taxa. In the evaluation stage, multiple image regions inside the particle are classified using the proposed method. Regions classified as aperture vote positively while the rest negatively. The votes are pixel-wise averaged on a confidence map which indicates the likelihood of the presence of an aperture. Finally, a local maxima algorithm allows to determine the position and amount of the apertures. The overall scheme of the method can be seen in Fig. 2. The application of individual classifiers in the same manner yields to the detection of the different apertures types, which combined results in the definitive detected aperture list.

Extraction of Apertures

Description of apertures based on Bag of Words

The foundation of the proposed strategy is the recognition of visual patterns present in the pollen apertures, that allow them to be described in terms of primitive subimages. This idea originated from the Bag-of-Words strategy (also known as Bag of Features) for text categorization, where a codebook is created from words contained in the analyzed text [START_REF] Csurka | Visual categorization with bags of keypoints[END_REF]. Then, classification is possible by observing the frequency of each word on the text. The visual counterpart and its extensions has been typically applied to object recognition [9][10]. However, to our knowledge, the proposed method is the first application of BOW to microscopic objects. First, different areas of the aperture are represented by a corresponding local descriptor. Although typical BOW indicates that such areas are determined by salient keypoints, an alternative is to compute a dense grid of fixed patches instead [START_REF] Nowak | Sampling Strategies for Bag-of-Features Image Classification[END_REF]. This technique has shown to be efficient by extracting more information from the whole analyzed object with an implied increase of computation. Because apertures regions are not expected bigger than 40 pixels, dense computation of patches is feasible. The current method employes a square grid of square patches, whose size can be determined according to the expected aperture size. Since the grid is set over the aperture region, most of the patches are informative. An example of a grid of 4x4 patches of 8 pixels by side (denoted as 4-8 region) is shown in Fig. 3. LBP as Local Descriptor. A local descriptor is needed to extract the visual pattern from each patch of the grid. Typically, descriptors are chosen based on its performance together with a keypoint detector, for example SIFT. However, since the proposed method does not need a keypoint detector, only local descriptors are concerned. A comprehensive survey of local descriptors can be found in [START_REF] Mikolajczyk | A performance evaluation of local descriptors[END_REF]. In previous experiments, local jets [START_REF] Schmid | Local grayvalue invariants for image retrieval[END_REF], a differential descriptor suitable for shapes, shown to be acceptable for the classification of aperture regions [START_REF] Lozano-Vega | Classification of Pollen Apertures Using Bag of Words[END_REF]. In the current work, Local Binary Patterns (LBP) operator is proposed and tested as local descriptor of the patches due to the gray-level and rotation invariance, in addition to the simple computation and short representation.

LBP is a simple but efficient operator that condenses the spatial pixel arrangement in a small image vicinity [START_REF] Ojala | A comparative study of texture measures with classification based on feature distributions[END_REF]. LBP has been widely employed as texture descriptor in complex tasks, for example, in face recognition [START_REF] Huang | Local Binary Patterns and Its Application to Facial Image Analysis: A Survey[END_REF]. Considering a central pixel {x,y} as the reference position, its grey level is compared one-to-one to the levels of neighbor locations. The individual comparison results are counter clockwise concatenated as bits into a binary code, which take the value 0 if the reference pixel has a greater gray level than the corresponding neighbor pixel, and 1 otherwise. This strategy makes the operator invariant to monotonic gray-level change on the image. The idea can be extended to a neighborhood of radius r and limiting the number of compared surrounding pixels (and consequently the bit length of the code) to p pixels yielding to a code denoted as LBP p,r . A graphical explanation of the obtaining of a LBP code is shown in Fig. 4. Due the bit order, similar pixel arrangements result in different codes when the area is rotated. Ojala et al. noticed this fact and reduced the possible LBP codes by grouping rotated patterns in a unique code, providing additionally rotation invariance [START_REF] Ojala | Gray scale and rotation invariant texture classification with Local Binary Patterns[END_REF]. Furthermore, they defined a set of uniform patterns (u 2 ) in which the transitions among bits (0/1 change) on the LBP code is not oftener than 2. The rest of the patterns were grouped in a single code. Ojala et al.

showed that the rotation-invariant uniform patterns, denoted as LBP riu2 p,r , are still effective since they represent on average 90% of the textural information of an image. The final amount of possible codes is reduced considerably, for example, a LBP riu2 8,r has only 10 possible binary codes. In the proposed method, the LBP descriptor of an image patch is obtained by estimating the normalized LBP code histogram of all the pixels, except those on the border according to the radius r. The LBP riu2 8,1 10-value descriptor was chosen for the proposed method due to the reduced size of the patches. Pollen images are expected to be scanned at the same scale and therefore, strong scale invariance is not required. Moreover, the change of the size of the apertures can be neglected for pollen of the same type.

Creation of the Codebook of Visual Words. Due to the natural variability and the observation point, apertures of the same type cannot show identical LBP patch descriptors. In order to classify the same visual patterns having similar LBP representation , it is necessary to create groups of similar descriptors. Each created group is the equivalent of visual words of the BOW. Assigning a word of the codebook to each LBP patch descriptor is a process that requires clustering techniques. In this work k-means clustering is employed due to its general effectivity on multiple clustering problems and to its flexibility to chose the number of clusters k, i.e., the size of the codebook. On early work, it was demonstrated that adding spatial information of the patches improves the classification performance of the regions [START_REF] Lozano-Vega | Classification of Pollen Apertures Using Bag of Words[END_REF]. For this reason the coded location of the patch was incorporated to the final patch descriptor. Taking advantage of the source grid where patches are regularly spaced, the position of the patch inside the grid is simply estimated as a number ranging from 0 to 1. Once created, the cluster model will serve to match visual words to unseen image patches. In the experimental section, variation of the size k of the codebook is evaluated with respect to the performance of classification of regions. In Fig. 5, examples of visual words are depicted with some instances of their image patches. Aperture descriptor The aperture descriptor is composed by 2 elements: the content of visual words and the spatial information relative to the particle. The first element is given by the computation of the normalized histogram of visual words of the aperture. This means that each bin of the histogram describes the frequency of each of primitive subimage inside the aperture region. Similar apertures will have a very similar histogram distribution. The codebook is directly related to the histogram: the bigger the codebook, the larger the histogram. The second element of the description consists of the relative distance of the aperture with respect to the centroid of the particle. Apertures change their appearance according to the orientation of the particle with respect to the focal plane, which also affects the position of the aperture on the image. By considering this position, information about the orientation of the apertures is also provided to the classifier. The position is computed as the relative distance from the central point of the aperture region to the centroid of the particle (RDC ), where a distance 0 corresponds to the centroid itself and a distance 1 corresponds to the outermost point of the particle. Fig. 6 shows an example of the RDC computation of some regions inside a particle. 0 .2 0 .9 0 . 5 Fig. 6. Example of some regions at different position inside the particle and its relative distance to the centroid RDC.

Individual Classifiers

An important characteristic of the proposed method is the possibility of easily add a new type of aperture or a set of aperture appearances from a single pollen type, using the same algorithm and without the need of classifying everything already learnt over again. The task is split modularly into individual classifiers for each added aperture type and they eventually will be combined into an aggregate detection. This means that when a new aperture needs to be added to the detector, only the new module has to be created and combined with the existing modules of the rest of the apertures. Since the total of detected apertures is the sum of the detections of the individual classifiers, the problem has optimal substructure property and, therefore, the optimal solution can be obtained by finding independently the optimal solution for each of the individual classifiers.

It is possible to find different appearances of the same aperture type, and therefore, a classifier able to handle such complexity is needed. In a previous work, it was shown that the Support Vector Machine (SVM) method was suitable for classification of apertures with multiple appearances [START_REF] Lozano-Vega | Classification of Pollen Apertures Using Bag of Words[END_REF]. SVM is chosen due to proven results in binary classification based on mapping data into a higher dimension to enable linear separation [START_REF] Byun | Applications of support vector machines for pattern recognition: A survey[END_REF].

Therefore, an individual classifier is trained as a binary SVM classifier, where aperture regions are regarded positive and regions not belonging to an aperture are regarded as negative. In order to guarantee that the classifier is able to discriminate both classes among all possible pollen types, samples of negative regions from all pollen types are included in the training dataset. The mapping kernel employs the radial basis function (RBF), and the parameter γ and the penalty C are tuned in a 3-fold cross-validation.

Detection of apertures from an unknown image

The SVM classifier delivers a confidence value for each analyzed region. The value can be interpreted as the level of belonging of the region to the positive class. A final binary classification is obtained by comparing the confidence value to a given decision threshold (d thr ) in order to obtain the decision for a positive (aperture) or negative class (not an aperture). Multiple regions covering an specific area can vote, according to their assigned class, to reinforce the belief of the existence or absence of an aperture in the particle. In the proposed method, several overlapping regions covering the whole particle are classified and vote to build an aperture confidence map of the particle. High confidence areas of this map represent detected apertures.

The method does not require mandatorily prior segmentation or knowledge of the particle to be evaluated. However, we have added a fast trace of the contour in order to avoid the analysis of useless regions outside the particle contour. Similarly to the thresholding of Chen et al., our proposal employs Otsu thresholding for binarizing the image which is nonparametric, unsupervised and fast. The principle of the method is to find a threshold that minimizes the intraclass gray-level variance of the background and the foreground pixels. From the resulting binary image, the contour is estimated with Suzuki's method [START_REF] Suzuki | Topological Structural Analysis of Digitized Binary Images by Border Following[END_REF] using the OpenCV implementation [20]. Excellent performance can be obtained, provided that the source image exhibits a good contrast and there is no debris stuck to the pollen.

A set of overlapping square regions of equal size are sampled from the particle's inner area defined by the contour. Then, the aperture descriptor from section 2.2 is computed on each sampled region of the particle. Classification confidence is obtained after applying the individual classifier and the decision threshold ( d thr ) is employed in order to have a binary decision.

In order to merge the aperture information of the sampled regions, a confidence map is created the same size that the particle image. Each classified region contributes with its vote pixel-wisely to build the map. Regions with positive class votes with value 1 and negative with value 0. The final value of the confidence map is the averaged contribution of all the regions votes expressed as follows in Eq. 1:

CM (x, y) = n i=1 V i n (1)
where CM (x, y) is the value of confidence map for the pixel located in (x, y), V i is the vote of the region i, and n is the total of classified regions that contain the pixel in (x, y). The greater n, the less impact of a single region on the final value, making the estimation of the map more resilient to misclassification of the regions. The resulting dense confidence map can be interpreted as the likelihood of the presence of apertures on the particle. Therefore, high values on the map are regarded as detected apertures.

The procedure to detect apertures is to find local maxima on the confidence map. Because low confidence areas can cause false detections due to low local maxima, only locations whose confidence value is above a confidence threshold (c thr ) are analyzed. Before searching local maxima, the confidence map is smoothed with a gaussian filter of one quarter the size of the regions, in order avoid noise from irregularities due to region sampling. Finally, the (x, y) positions of the local maxima are regarded as potential detection of apertures.

A final procedure is applied before considering the potential detected apertures as definite. If more than one detection are closer each other than a given detection distance (d dist), the detections are regarded of the same aperture, and therefore merged. In such case, the definitive single detected aperture is located at the middle point of the former potential apertures. This technique helps to reduce multiple detection on the same aperture. Detected apertures from individual classifiers are consolidated together to yield to the definite list of detected apertures of the system. Fig. 7 shows examples of the confidence map of different particles and the position of the detected apertures.

Evaluation of the detection

The apertures change gradually their appearance depending on the observation point, from barely visible to clear and well defined. The proposed method is expected to be trained with enough defined apertures in order to reduce the false positive rate. However, it is still possible to detect unexpected apertures that are not well defined due to their similarity to the well defined apertures.

Consequently the ground truth of apertures, which is used to evaluate the performance of the detection, is split into two sets. The ground truth set (GT) groups those enough defined apertures for which the method was trained and whose misdetections are penalized in the evaluation. The unexpected set (UE) groups those not well defined apertures which were not expected to be detected by the algorithm and whose misdetections are not penalized in the evaluation.

Examples of both sets are shown in Fig. 8 Alder Birch Mugwort The evaluation of the detection consists of quantizing the match between detected apertures and the ground truth of the location of apertures. A distance matrix is employed to compare the detections (DT) on one side, to the group defined by GT and UE on the other side. The euclidean distance is computed one-to-one between the apertures of both compared groups. An example of this matrix is shown in Fig. 9a. Each DT aperture is matched to the closest GT or UE aperture as long as the distance is below an evaluation threshold (ev dist) [Fig. 9b]. A true positive (hit) is regarded for each GT aperture that was matched to at least one detection. A false negative (miss) is regarded for each GT aperture that was not matched to any detection. An unexpected true positive (unexpected hit) is regarded for each UE aperture that was matched to at least one detection. A false positive (false alarm) is regarded for each detection that not matched to any GT or UE aperture and for each multiple detection (2 or more detections of the same aperture) [Fig. 9c]. True positive of UE apertures are not considered for the total of true positives, and their false negative are neither penalized since the method did not learn to detected them; and the statistics are computed just to gain additional insight into the performance. With these statistics; true positive, false negative and false positive; is possible to compute a curve similar to the Receiver Operating Characteristic (ROC) that compares the true positive rate (recall ) versus the false positive rate (fallout) by varying c thr. This kind of plot is useful to visualize the performance of the system under different conditions. If the specific cost function is known for the problem, a point (or segment) on the ROC curve becomes optimal [START_REF] Provost | Robust Classification for Imprecise Environments[END_REF]. Contrary to a typical ROC curve, where the threshold variation applied to a classification confidence value directly affects the trade-off between true positives and false positives, the computation of the ROC curve for of the aperture detection is result of a more complex process after the application of the c thr (local maxima, refinement and matching) and the trace of the curve could be not always smooth. Additionally, since there is no a ground truth for the negative class (actually it is represented for all the regions of the particle that are not an aperture), the total of detected apertures is considered for the computation of the false positive rate, instead of the typical total of negatives. This change makes more sensible the ROC curve to a small number of detected apertures, affecting the left-bottom side of the curve.

Experiments and Evaluation

The proposed method was evaluated with focus on the 5 most allergenic pollen types in Germany: Alder, Birch, Hazel, Mugwort and Grass. Although the 5 types possess apertures, the tested system was trained to detect apertures of three of them: Alder, Birch and Hazel. Mugwort apertures are hardly visible since they are usually hidden in a furrow. Grass apertures are small compared to the size of the pollen, fact that makes them rarely visible in most of the focal planes. For this reason, the amount of Grass apertures in our database was vey limited for training. However, the incorporation of Mugwort and Grass samples in the experiments is important in order to evaluate the robusteness to multiple pollen and to simulate similar conditions to a real-world application.

Images of pollen were scanned with a brightfield microscope Keyence BZ-9000 with magnification 40x and a resolution of 0.26 microns/pixel. Particles were stained with magenta dye. The following amount of particles for each taxa were employed: Alder 115, Birch 136, Hazel 52, Mugwort 120 and Grass 132. The total of ground truth labeled apertures is 541.

The presented methodology of aperture description, detection, localization, counting and evaluation was applied on the database and 3 individual classifier were built for Alder, Birch and Hazel. In order to test robustness to unseen data, a 4-fold cross validation was employed on the input dataset, meaning a training/testing ratio of 0.75/0.25.

Evaluation of region parameters and codebook size

For the first evaluation step, the correct configuration for patch size, grid size and number of visual words of the codebook needs to be selected based on the performance of the SVM classification of the regions.

Since the variation of the decision threshold d thr is evaluated together with the confidence threhshold c thr in the matching phase, the selected combination has to evidence the best performance at different levels of d thr instead of just a fixed value. The Area Under the Curve (AUC) of the ROC curve of the region classification was employed as comparison metric in order to chose the best combination, robust to different levels of the d thr. For this metric, a perfect performance for all d thr obtains an AUC = 1, and an average random perfomance obtains an AUC = 0.5.

Different region sizes, ranging from 24 to 38 pixels, along with their combination of grid size-patch size and number of visual words were tested. The evaluation was repeated for each individual classifier. Fig. 10(a-c) shows the behavior for each classifier of AUC for different combinations. For the 3 individual classifiers, 25 visual words is generally the best and most consistent selection of codebook size. Although more visual words may improve the performance, it also increases the computation time, especially for a great amount of regions. Fig. 10d compares AUC for the 3 classifiers using 25 visual words. Best results seem related to bigger region sizes. Furthermore, for greater grid sizes (5 to 6), performance seems more stable. Particularly, the 6-6 region (36 pixels) shows better performance and this was the configuration used for the aperture detector. Similar combinations, like the 5-7 region could have been good election too.

Evaluation of classification performance

Using a 6-6 region configuration, the pollen dataset was sampled, ranging from 30 to 60 overlapping regions per particle, depending on the size. A SVM model was created from the training set for each fold of the cross validation scheme. The Radial Basis Function (RBF) was employed as kernel, the SVM parameter γ and penalty C were optimized using a grid approach with AUC of the ROC as performance metric. A classification confidence value for each sample was computed after applying the SVM model to the testing set. The learning process was reapeated for each individual classifer.

In order to select the values of d thr and c thr, each classifier was evaluated on the detection of apertures and a family of ROC curves was plotted. Fig. 11 shows example of the curves. Distances d dist and ev dist were set experimentally to 27 and 35 respectively. Assuming a balanced cost function, where the cost of missing a detection is equal to the cost of a false positive, the optimal performance is found in the closest point to the upper left corner of the ROC, where the optimal d thr and c thr can be determined. Due to the imbalanced datasets, the d thr is offset towards zero, being the optimal value relatively low.

As a final step, the output of the optimized individual classifiers were combined to create the final list of detected apertures from the 5-taxon dataset. The final list was evaluated against the ground truth of all the pollen types and the following performance statistics were computed as average of the 4-fold evaluation with confidence level of 0.95. For apertures of the 3 learnt pollen types, joint performance showed a true positive rate of 0.78 ± 0.03 , a false positive rate of 0.23 ± 0.05, and a precision of 0.76 ± 0.05. Considering all the 5 pollen types, the true positive rate remains 0.78 ± 0.03, but the false positive rate moved to 0.29 ± 0.06, and the precision to 0.70 ± 0.06. The difference is caused to false positive detections on Mugwort and Grass. There is no evidence that other detectors of aperture have been evaluated, and therefore there is no baseline for comparison. However, the resuts are promising considering the high difficulty of learning 3 different aperture types with multiple appearances and adding 2 pollen types which apertures are not easily visible.

Conclusions and future work

We have proposed a flexible image-based method to detect, localize and count apertures of pollen. The method is capable of learning independently apertures from different pollen types employing individual classifiers and with the same algorithm. Furthermore, the BOW approach makes the method robust to changes in appearance of the same aperture type. Additionally, the proposed region sampling strategy allows to adjust the configuration for each aperture type.

Experiments confirmed suitability of the method in similar conditions to a real-world application by adding extra pollen taxa with no visible apertures. The cross-validation evaluation showed that the method is robust to unseen particles. The method succeeded under total lack of previous knowledge about location and amount of the apertures and under challenging conditions such as intra-class variation, inter-class similarity, multiple input taxa and analysis on a single focal plane.

Because there are not performance results of the existing methods, there is no baseline for comparison among methods. However, aperture counts, computed using the proposed method, have been recently employed in a pollen classification system together with shape, texture and size features [START_REF] Lozano-Vega | Analysis of Relevant Features for Pollen Classification[END_REF]. The performance was comparable to that of Chen et al. [START_REF] Chen | Feasibility study on automated recognition of allergenic pollen: grass, birch and mugwort[END_REF] and superior to Boucher et al. [START_REF] Boucher | Development of a semi-automatic system for pollen recognition[END_REF], who have also employed aperture detection on their pollen classification experiments.

Since apertures change their appearance when observed at different focal planes, next step is the extension of the method for multiple images gathered at these planes. One path is to apply the method to multiple focal planes independently. The advantage would be the confirmation of apertures already detected and the detection of missed apertures. A second path is the application of the lo- cal descriptor involving multiple focal planes at the same time, maybe employing a volume LBP operator instead, which describes the volumetric aperture.

Fig. 1 .

 1 Fig. 1. Examples of most important allergic pollen in Germany. Apertures of Alder, Birch and Hazel pollen are highlighted by a green circle. Note the different appearances for the same pollen type. Although Mugwort and Grass also own apertures, they are rarely visible on the image due their size or position on the particle.

Fig. 2 .

 2 Fig. 2. Method of representation of apertures and detection on unknown particles. a) Aperture descriptor: LBP codes computed on a dense grid of aperture regions are clustered to create the codebook of visual words. Then, apertures are described by the histogram of visual words combined with the relative distance to the centroid (RDC ) as spatial information. b) Detection scheme of an aperture: Multiple overlapping regions are evaluated and classified to create the confidence map from which the aperture position is estimated.

Fig. 3 .

 3 Fig. 3. Example of a 4x4 grid employed for the computation of local descriptors. The patch size is a square of 8 pixels.

Fig. 4 .

 4 Fig. 4. Example of computation of the LBP code from a single pixel. a) Values of the gray levels of the 8 neighbors are shown overlapped to each pixel location. The value of the reference pixel is highlighted in red at the center. b) Individual binary results of the comparison to the reference value. c) Final concatenated LBP code.

Fig. 5 .

 5 Fig. 5. Example of 6-pixel patches grouped according to the assigned visual word. Note the similarity of the patches of the same visual words. Patches are amplified for visualization purposes.

Fig. 7 .

 7 Fig. 7. Example of confidence maps overlapped over the original pollen image Alder (left), Birch (center) and Hazel (right) pollen. Warm colors (red-orange) indicates high values of the map, regarding a high probability of the presence of an aperture. On the contrary, cold colors (blue-green) indicate low values of the map. The position of the detected aperture after finding local maxima is shown by a blue mark.

Fig. 8 .

 8 Fig. 8. Examples of the aperture sets of Alder, Birch and Mugwort used for evaluation. Well defined apertures marked by a red arrow pointing left belong to the ground truth set (GT). Not well defined apertures marked by a blue arrow pointing right belong to the unexpected set (UE).

Fig. 9 .

 9 Fig.9. Example of evaluation of 4 detected apertures (DT) against 3 ground truth (GT) and 1 unexpected (UE) apertures. The distance matrix show the euclidean distance between the groups (a). Minimum distances are highlighted in bold. If the distance is below an ev dist = 27 pixels, the DT aperture is matched as 1 in the matrix of matches (b). Finally, evaluation statistics are extracted form the matches (c).

Fig. 10 .

 10 Fig. 10. AUC performance of region SVM classification for indivudual classifiers: a) Alder), b)Birch and c)Hazel. d) A comparative of the same classifiers just for 25 visual words.

Fig. 11 .

 11 Fig. 11. ROC curves of evaluation on aperture detection for indivudual classifiers: a) Alder) b)Birch c)Hazel. For clarity, only curves for fold 4 are shown as an example. Different points of the curve are created by variating c thr. Variation of d thr is represented by different curves.

Acknowledgements. The authors are grateful to Celeste Chudyk, Yann Ryann and Morad Larhriq for scanning and preparing the pollen datasets and to the Max Plank Institute for Chemistry for permitting the use of their facilities.


Ignacio Salas
email: ignacio.salas@mines-nantes.com

Gilles Chabert
email: gilles.chabert@mines-nantes.com

Alexandre Goldsztejn
email: alexandre.goldsztejn@univ-nantes.fr

The non-overlapping constraint between objects described by non-linear inequalities

Packing 2D objects in a limited space is an ubiquitous problem with many academic and industrial variants. In any case, solving this problem requires the ability to determine where a first object can be placed so that it does not intersect a second, previously placed, object. This subproblem is called the non-overlapping constraint. The complexity of this non-overlapping constraint depends on the type of objects considered. It is simple in the case of rectangles. It has also been studied in the case of polygons. This paper proposes a numerical approach for the wide class of objects described by non-linear inequalities. Our goal here is to calculate the non-overlapping constraint, that is, to describe the set of all positions and orientations that can be assigned to the first object so that intersection with the second one is empty. This is done using a dedicated branch & prune approach. We first show that the nonoverlapping constraint can be cast into a Minkowski sum, even if we take into account orientation. We derive from this an inner contractor, that is, an operator that removes from the current domain a subset of positions and orientations that necessarily violate the non-overlapping constraint. This inner contractor is then embedded in a sweeping loop, a pruning technique that was only used with discrete domains so far. We finally come up with a branch & prune algorithm that outperforms Rsolver, a generic state-of-the-art solver for continuous quantified constraints.

Introduction

The goal of this article is to calculate the set of all positions and orientations that can be given to an object so that it does not overlap a second one (see the left graphic in Figure 1, which shows the simpler case where no orientation is considered). We address the general case of objects described by nonlinear inequalities. Calculating this set is a key task for solving packing problems, that consists in placing a set of objects in a bounded space so that they do not overlap pairwise.

In this introduction, we first define objects in our context and give a precise statement of the non-overlapping constraint. Then, we explicit the type of objects we consider and explain what we mean by calculating a set. We finally mention related works. In Section 2, we show that our problem can be cast into the calculation of a Minkowski sum. Based on this observation, we propose a branch & bound algorithm in Section 3. Experimental results are shown in Section 4 and a conclusion follows.

Object definition

For clarity, let us first assume that the orientation is fixed, that is, objects can only be translated.

Translating an object means fixing the position of a particular point that we call the origin. This origin point can be arbitrarily chosen. For instance, in rectangle packing, the origin of a rectangle can be a vertex or its center point. Once this convention for the origin is made, the shape of the object is just a regular constraint. To illustrate this, let us consider again rectangle packing. If the origin is the lower-left corner, then a rectangle of dimensions l 1 and l 2 is the set of all p ∈ R 2 satisfying c(p) ⇐⇒ 0 ≤ p 1 ≤ l 1 ∧ 0 ≤ p 2 ≤ l 2 .

(1)

Alternatively,

c(p) ⇐⇒ - l 1 2 ≤ p 1 ≤ l 1 2 ∧ - l 2 2 ≤ p 2 ≤ l 2 2 . ( 2
)
defines a rectangle with the center point as origin, which, of course, is just a shift of the previous constraint. So, the shape of an object can be expressed as a constraint, the latter containing an implicit convention for the origin. As a last example, a circle of radius r is the set of all p ∈ R 2 such that

c(p) ⇐⇒ p ≤ r (3)
the origin being, in this case, the center of the circle.

This definition of c(p) corresponds to an object with no translation nor rotation. The general constraint corresponding to an object translated by x and rotated by α is easily obtained as follows. In the case of translation only, the part of the plane covered by an object placed at some position x is the set of all p such that c(px) is satisfied. Let us introduce orientation. By a classical geometric argument, an object placed at x and turned by some angle α is the constraint

c R -α (p -x) (4)
where R α is the rotation matrix with angle α:

R α = cos(α) -sin(α) sin(α) cos(α) . (5)
Equivalently, we can say that c(p) ≡ c(R 0 (p-0)) represents the object placed at 0 and rotated by the angle 0. To conclude:

an object is a constraint on the plane (i.e., with two variables), -placing an object means fixing the coordinate of its implicit origin, -orienting an object means fixing the angle of the rotation around its implicit origin.

In this paper we consider objects described by nonlinear inequalities c(p) ⇐⇒ f (p) ≤ 0. E.g., a circle of radius 1 which origin is the center point is the set of all p ∈ R 2 such that f (p) ≤ 0 with f : p → p -1. For the clarity of presentation, we will assume each object to be described by a single inequality, but our results can be easily extended to the more general case of objects described by disjunctions of conjunctions of inequalities. This is possible by introducing min and max operators. For instance, (f 1 (x) ≤ 0 ∨ f 2 (x) ≤ 0) is equivalent to max{f 1 (x), f 2 (x)} ≤ 0. Differentiability is not required by our algorithm. In fact, there is also no assumption on the functions involved, except that they are defined by mathematical expressions based on usual operators (+, ×, √ . , exp, etc.). In particular, there is no convexity assumption on the input objects.

The non-overlapping constraint

We can now focus on the non-overlapping constraint. The non-overlapping constraint involves two objects, one being fixed, the other representing the unknowns of the problem. For this reason, we will call "reference object" the first one and denote by c R its describing constraint. The second one will be called "moving object" and its constraint denoted by c M .

We will show at the end of this section that the general case, where both objects are translated and/or rotated, can actually be obtained from the simpler case where the reference object has no transformation. Intuitively, the frame where the overlapping constraint is stated can be centered on the reference object and aligned with its orientation, albeit the exact formula is not so trivial.

So, we shall only introduce in our definition below the position x M and the rotation angle α M of the moving object. The non-overlapping constraint is the negation of the overlapping constraint that can be stated as follows.

Definition 1 (Overlapping constraint). Given two constraints c R and c M , a vector x R ∈ R 2 and α R ∈ [0, 2π]:

overlap (c R ,c M ) (x M , α M ) ⇐⇒ ∃p ∈ R 2 , c R (p) ∧ c M R -α M (p -x M ) . (6)
In the case of translation only, this simplifies to

overlap (c R ,c M ) (x M ) ⇐⇒ ∃p ∈ R 2 c R (p) ∧ c M (p -x M ). ( 7
)
Our goal is to calculate the overlapping constraint. By calculating, we mean here that an explicit (numerical but verified) representation of the solution set

S ′ := {(x M , α M ), overlap (c R ,c M ) (x M , α M )} (or S := {x M , overlap (c R ,c M ) (x M )}
in the case of translation only) has to be returned by our algorithm.

We show now that the overlapping constraint in the case where the reference object is given a position x R and an orientation α R can be tested using S ′ (and hence using its explicit representation). More precisely: Proposition 1. The moving object with position x M and orientation α M overlaps the reference object with position x R and orientation α R iff

(R -α R (x M -x R ), α M -α R ) ∈ S ′ (8)
Proof. By definition, (x M , α M ) satisfies the overlapping constraint with the reference object at position x R and orientation

α R iff exists ∃p ∈ R 2 such that c R (R -α R (p -x R )) and c M (R -α M (p -x M )). This is equivalent to ∃q ∈ R 2 , namely q = R -α R (p -x R ) ⇐⇒ p = R α R q + x R , (9)
such that c R (q) and

c M (R -α M (R α R q + x R -x M )) ⇐⇒ c M (R -α M +α R (q + R -α R (x R -x M ))). ( 10
)
This is equivalent to (8).

Remark 1. The non-overlapping constraint is usually considered for packing applications. Its description is easily obtained from the one of the non-overlapping constraint. The former is preferred here because it simplifies the constraints expressions and its link with the Minkowski sum presented in Section 2.

Intervals, Boxes and Paving

The representation we use is called a paving. This representation is a natural choice in the context of constraint programming where the large majority of algorithms dedicated to continuous variables, if not to say all, assume domains of variables to be intervals (see, e.g., [START_REF] Benhamou | Continuous and interval constraints[END_REF][START_REF] Jaulin | Applied Interval Analysis[END_REF]). In the following definition, we call box a Cartesian product of d intervals where d is either 2 in the case of S or 3 in the case of S ′ .

Definition 2 (Paving). A paving of a set S ⊂ R d is a triplet (I, B, O) where I (for "inside"), O (for "outside") and B (for "boundary") are three sets of boxes verifying

∪I ⊂ S, (∪O) ∩ S = ∅ and ∪ (B ∪ I ∪ O) = R d . ( 11
)
An example of paving is shown in Figure 2.

Contribution and Related Works

This paper proposes an algorithm that computes a paving of the overlapping constraint. One contribution is on the modeling aspect of this problem: we show that the overlapping constraint can be expressed as a Minkowski sum. On the one hand, this generalizes the approach and simplifies the description of the algorithms. On the other hand, it allows handling the simple translation case and the more involved translation and rotation case homogeneously by casting the rotation to a translation into an augmented space (see Proposition 2). The other contribution is algorithmic. We propose an original inner contractor for this problem, that is, an operator that identifies parts of the solution set. This operator implements a sweeping loop and exploits the properties of the Minkowski sum. The second operator is an outer rejection test based on classical constraint propagation. Both are interleaved in a branch and prune algorithm that computes the desired paving.

From another point of view, definition 1 means that our problem falls into the category of existentially-quantified constraints. A state-of-the-art algorithm for calculating a paving with existentially-quantified inequalities is given in [START_REF] Ratschan | Efficient Solving of Quantified Inequality Constraints over the Real Numbers[END_REF] and is implemented in the Rsolver tool [Rat] (Rsolver implements a general algorithm, somehow a numerical version of CAD, for solving quantified constraints).

General techniques [START_REF] Goldsztejn | Inner and Outer Approximations of Existentially Quantified Equality Constraints[END_REF][START_REF] Ishii | Interval-Based Projection Method for Under-Constrained Numerical Systems[END_REF] for quantified equality constraints could also be used: either by adapting [START_REF] Goldsztejn | Inner and Outer Approximations of Existentially Quantified Equality Constraints[END_REF] to compute an over approximation of the boundary of the overlapping constraint, or using a necessary condition for the boundary of the overlapping constraint expressed as an under-constrained system of equations and using [START_REF] Ishii | Interval-Based Projection Method for Under-Constrained Numerical Systems[END_REF]. However, the overlapping constraint naturally involves inequality constraints and using costly techniques dedicated to equality constraints turns out to be counterproductive.

Finally, we shall mention that an exact formula for the overlapping set S has been given in [START_REF] Beldiceanu | Non-Overlapping Constraints between Convex Polytopes[END_REF] in the case where objects are polytopes. The set, in this case, is the convex hull of the points obtained by summing one vertex of the first polytope to one vertex of the second one.

Overlapping as a Minkowski Sum

In this section, we show that the overlapping constraint can be reformulated as a Minkowski sum. This relation underlies our branch & bound solver, that will be presented further. Let us first recall the definition of the Minkowski sum of two sets:

Definition 3 (Minkowski sum). Given two equi-dimensional sets S 1 , S 2 ⊆ R d , the Minkowski sum is S 1 + S 2 = {x 1 + x 2 ∈ R d : x 1 ∈ S 1 , x 2 ∈ S 2 }. ( 12
)
The Minkowski difference is defined accordingly by:

S 1 -S 2 = {x 1 -x 2 ∈ R d : x 1 ∈ S 1 , x 2 ∈ S 2 }. ( 13
)
The right graphic of Figure 1 (page 2) shows an example of two sets and their Minkowski sum.

Considering S 1 as a constraint c 1 (i.e., x ∈ S 1 ⇐⇒ c 1 (x)) and, similarly, S 2 as c 2 , we have, equivalently:

S 1 + S 2 = {x ∈ R d : ∃p ∈ R d , c 1 (p) ∧ c 2 (x -p)}. ( 14
)
where d is the number of variables in the constraints. Comparing ( 14) and (7), we immediately see that S = S R -S M , i.e. the overlapping constraint can be represented as a Minkowski sum in the case of translation only.

We show now that S ′ is also a Minkowski sum, a less trivial result. To this end, we embed the moving object S M ⊆ R 2 into R 3 encoding its rotation within the additional dimension:

S ′ M := {(v, β) : c M (R β v)} = {(v, β) : R β v ∈ S M }. (15)
Now, the following proposition states that the overlapping constraint with rotation S ′ can be written as a Minkowski difference of two such "augmented" sets (see Figure 3).

Proposition 2.

S ′ = S R × {0} -S ′ M . ( 16
) Proof. By definition, (x M , α M ) ∈ S ′ holds if and only if ∃p ∈ R 2 such that c R (p) and c M (R -α M (p -x M )). Equivalently, there exists u R ∈ S R and u M ∈ S M such that u R = p and u M = R -α M (u R -x M ) ⇐⇒ x M = u R -R α M u M .
Finally, the vector (x M , α M ) is proved to be the sum of (u ′ by an inner contractor C in , that is, an operator that guarantees:

R , 0) ∈ S R × {0} and (R α M u M , α M ) ∈ S ′ M .
[x] ′ ⊆ [x] ∧ [x]\[x] ′ ⊆ I. (17) If [x] ′ = ∅ then [x] ′ is contracted to a box [x]
′′ by an outer contractor C out that guarantees:

[x] ′′ ⊆ [x] ′ ∧ [x] ′ \[x] ′′ ⊆ O. (18) Finally, if [x] ′′ = ∅ then [x]
′′ is bisected in two new boxes that are pushed in the list of boundary boxes B. The recursion stops when the total surface of boxes in B is less than ε% of the initial box surface (s 0 ), ε being a user-defined parameter.

The originality of our approach lies in the inner contractor that we describe now. It is dedicated to the handled problem as it makes use of the specific structure of the quantified constraint (14).

Inner contractor

Our inner contractor is based on the inner arithmetic and a sweep loop: The former builds small inner boxes, and the latter makes the union of these boxes in order to remove slices of the initial box, leading to a so called inner contraction. The inner arithmetic is a variant of the classical interval arithmetic that allows to build a sub-box of a box [x] that is inside a given set S described by inequalities. This technique was first introduced in ➜3 of [START_REF] Chabert | Sweeping with Continous Domains[END_REF] and used in [START_REF] Araya | Upper Bounding in Inner Regions for Global Optimization under Inequality Constraints[END_REF] in the context of global optimization. This inner arithmetic can also be used with a initial point (or initial box) that is inflated, that is to say, given a box [x] and x ∈ [x], it produces a box [x] such that

Algorithm 1: (I, B, O) = pave([x], ε) s0 ← surface([x]); // initial
[x] ′ ← Cin([x]); I ← I ∪ ([x]\[x] ′ ); // inner contraction 8 [x] ′′ ← Cout([x] ′ ); O ← O ∪ ([x] ′ \[x] ′′ ); // outer contraction 9 if ([x] ′′ = ∅) then 10 ([x]1, [x]2) ← bisect([x
x ∈ [x] ⊆ [x] ∧ [x] ⊆ S, (19)
or an empty box if x ∈ S. This arithmetic has similar properties to its classical counterpart: the time complexity is in the length of the constraint expression and gives an optimal box [x] (i.e., of maximal size in every dimension) if no variable occurs twice in the expression. Before describing how to contract a box [x] with this new arithmetic, let us first address a simpler question: how to find a subbox of [x] that is inside S ?

A possible answer is to look for two boxes [x] 1 and [x] 2 such that

[x] 1 ⊆ S 1 , [x] 2 ⊆ S 2 and ([x] 1 + [x] 2 ) ∩ [x] = ∅ (20) because, in this case, ([x] 1 + [x] 2 ) ⊆ [x] ∩ S.
To find such boxes, one can calculate in parallel two pavings, one of S 1 and one of S 2 , and stop the process as soon as two boxes satisfying (20) are found. By combining boxes of the first paving with boxes of the second one, this approach amounts to run a branch & bound in a (2 × d)-dimensional space. Note that this branch & bound is a sub-solver embedded in the main one, the one for the x variable. Our goal is to reduce the sub-solver to d dimensions only, which is, by the way, the incompressible price to pay for handling d existentially-quantified parameters.

To this end, we use the same idea as above, but this time based on Relation (14) (see Figure 4). To build an inner box in [x], let us first assume that we have picked some point x inside [x] (this point is, in fact, automatically yield by the sweep loop, as this will be explained in Figure 5). Then we look for another point p such that

c 1 (p) ∧ c 2 (x -p). ( 21
)
Finding this point is the task of the subsolver.3 Once p is found, it is "inflated" to a subbox [x] 1 of (p + [x]x) that is inside S 1 , which is possible with the inner arithmetic. The point (xp) is also inflated to a sub-box [x] 2 of ([x]-p) that is inside S 2 . However, if x turns out to be outside S, the last inflation cannot succeed and the process is interrupted in this case. Otherwise, the resulting box satisfies (

[x] 1 + [x] 2 ) ⊆ S and ([x] 1 + [x] 2 ) ∩ [x] = ∅. Note that ([x] 1 + [x] 2 ) ∩ [x] = ∅ is just a consequence of x ∈ [x]
. So the initial boxes used for both inflations are somehow arbitrary, but fixing them as we did is an heuristic that tends to maximize the surface of the final inner box

([x] 1 + [x] 2 ) ∩ [x].
Now that we have a technique to build an inner box inside [x] that contains a specific point x, we can use this service inside a sweep loop. The sweep loop can be simply viewed as a way to contract a box by "piling up" boxes until some face is entirely covered. This is quickly depicted in Figure 5. The interested reader may refer to [START_REF] Chabert | Sweeping with Continous Domains[END_REF] for further details.

Outer contractor

The outer contractor is less sophisticated than the inner one and acts as a simple rejection test: the box is either entirely discarded or kept intact.

Rejecting a box [x] means proving

[x] ⊆ S 1 + S 2 , that is ∀x ∈ [x] ∀p ∈ R d , ¬(c 1 (p) ∧ c 2 (x -p)). ( 22
)
This assertion can be checked by running the same subsolver we used for the inner contractor, except that the point x is replaced by the current box [x]. If the subsolver finds no solution, the previous assertion is proven. Note that only the coordinates of p are bisected, so the subsolver actually proves a stronger assertion if the contraction with respect to c 2 is not optimal (the actual assertion depends on the consistency level enforced by the contraction with c 2 ). Note also that the precision used in the subsolver is dynamically set to the width of [x] in order to have a somewhat uniform time spent by the subsolver throughout the global search (on x). This dynamic precision also ensures that the outer contractor is convergent, that is, it rejects any small enough boxes outside S.

One may be surprised by the simplicity of this rejection test and expect a more elaborated contractor for the outer region, inspired by what we did for the inner region. But the situation could be interpreted in the other way around. Since the overlapping constraint is in two dimensions only, an inner satisfiability test would probably fits, as long as it is fast and convergent. However, such a test amounts to prove for [x] ⊆ S 1 + S 2 the following assertion

∀x ∈ [x] ∃p ∈ R d , (c 1 (p) ∧ c 2 (x -p)).
and, contrary to (22), the ∀ and ∃ quantifiers are involved, which means that the problem is much harder. So, the inner contractor can be seen here as a way to make up for the lack of inner test.

Our previous argument is only based on running time. It is clear that an outer contractor could also lead to a more compact paving, but the size of the paving is anyway conditioned by the representation of the boundary so that a drastic gain on this aspect is not really expectable.

Experimental results

Experimental setup

The algorithm proposed in this article calculates a paving (I, B, O) of the overlapping constraint. The difficulty of this task mainly depends on three criteria:

variable occurrences: the number of times each variable appears in the expressions of the inequalities directly affects the efficiency of the contractors; the more occurrences a variable, the less efficient contractors. This is a wellknown drawback of the classical interval arithmetic that carries over the inner arithmetic (used by the inner contractor). -convexity: if objects are non-convex, the boundary of the non-overlapping constraint will be less smooth. So the paving will be more complicated (hence more time consuming), especially near the boundary. -degrees of freedom: that is, whether we take into account rotation or not.

Allowing rotation gives a problem of much higher difficulty for multiple reasons. First, the size of the paving is exponential in the number of degrees of freedom so we cannot expect to get a 3D paving within the time scale of a 2D paving. Second, the angle in the inequalities creates a lot of multi-occurrences (see Equations ( 4) and ( 5)) and considerably increase non-convexity by the introduction of trigonometric functions.

Our benchmark is based on these criteria. We have made two types of experiments. The first one is with translation only. We have considered three objects of increasing difficulty. Object ➏1 is a simple ellipsis:

Object ➏1 : (p 1 /2) 2 + p 2 2 ≤ 1. ( 23
)
Object ➏2 is an ellipsis rotated by some fixed angle. Objects ➏1 and ➏2

are obviously of equal complexity if rotation is a degree of freedom, but not if Fig. 6. Objects of increasing complexity. From left to right: objects ➏1, 2 and 3.

we limit ourselves to translation. This is because the rotated object introduces multi-occurrences for p 1 and p 2 :

Object ➏2 : 1.5 × p 2 1 + 1.5 × p 2 2 -p 1 × p 2 -0.2 ≤ 0. ( 24
)
Finally, the third object has a "peanut" shape. It cumulates multiple occurrences and non-convexity, as depicted in Figure 6:

Object ➏3 : (p 2 1 + p 2 2 ) 2 -2 × (p 1 × p 2 ) -0.02 ≤ 0. (25)
The non-overlapping constraint involves two objects: the "reference" one (which coordinates are fixed at the origin of the frame) and the "moving" one that represents the unknowns. We have considered all possible combinations with the three types of objects above, that is, the 6 first cases in Table 1. In the second set of experiments, we have introduced rotation and tested with two ellipsis and two "peanuts" (cases 7 and 8).

In each experiment, the paving process is interrupted when the total surface of the boundary B is less than ε% the surface of the initial box (initial domain for the variables), where ε is a user-defined precision parameter. We have applied the same policy with Rsolver, the tool we are comparing to.

Since the precision is in proportion of the initial domain surface, it should be noted that the quality of the paving depends also on the width of the initial enclosure. The larger the initial domain, the less precise the resulting paving. For this reason, to give ε a meaningful value, we have set in the experiments the initial box to a fairly accurate enclosure of the overlapping set S or S ′ (as it can be seen in Figure 7 and subsequent).

Results (without rotation)

We first compare in Table 2 the running times obtained by Rsolver and our algorithm for the 6 first cases, with a precision ε set each time to 3.25%. This choice for ε corresponds to the minimal value that gives no timeout with Rsolver. The pavings obtained are depicted in Figure 7.

We then provide in 8 a more detailed analysis for the two extreme cases (case 1 and 6), where ε varies from 10% downto 1%. They show some significant absolute performance gain, as well as some better asymptotic behavior with respect to RSolver. The results first confirm the presumed complexity levels of the different cases, since the "harder" instances indeed require more time to be solved. They also show that our algorithm is more competitive than Rsolver. But, of course, Rsolver is a generic solver that does not take advantage of the specific structure of the handled problem. It should also be noted from the graphics that the gap between our approach and RSolver, in a given case, increases as we use smaller values for the precision.

Results (with rotation)

We only present here preliminary results with rotation.

Figures 9 and 10 show 2D sections of the 3D pavings we have obtained in cases 7 and 8 with a precision set to 3.25%. A 2D section is obtained by fixing the angle to some value and selecting the boxes in the 3D pavings for which the dimension of the angle contains this value. The two other dimensions are plotted.

We have set the domain of the angle to the interval [0, 0.7] for the case 7, and to the interval [0, 0.3] for the case 8. The only paving that was possible to obtain within the time limit was with our algorithm and for the case 7. This paving has been calculated in 9 minutes whereas RSolver does not return after 80 minutes. Both algorithm do not terminate after 100 minutes in the case 8.

When a program timeouts, only a partial result is displayed. A partial result means that the surface of B exceeds ε% the initial width.

The main purpose of this experiment is to show that our approach is still valid in the case where rotations are taken into account. Rotations only transform the expressions that describes the objects, without requiring any change in the algorithm itself. It is clear however that calculating a full 3D paving is a heavy task, whatever the algorithm is.

Conclusion

In the case of objects defined by non-linear inequalities, the non-overlapping constraint can only be handled numerically. In this paper, we have given an efficient way to generate verified pavings approximations for this constraint. These pavings represent explicitly the constraint, that is, the set of all acceptable positions and orientations of an object with respect to another. Our preliminary experiments have shown strong efficiency gains with respect to the state of the art solver for quantified numerical constraints RSolver, in particular in the case where the orientation of the objects is taken into account. In our future work, these pre-computed pavings will be incorporated within a packing algorithm, as explicit descriptions of the overlapping constraints. However, with rotation, we have seen that full 3D pavings are clearly too big so that a more adaptative approach will probably have to be considered as well. The idea would be to calculate on-the-fly sub-sets of the overlapping constraint, depending on the actual domains of the objects positions assigned by the packing solver.

Fig. 1 .

 1 Fig. 1. Non-overlapping constraint. Left: two objects SR and SM , the region in red represents the set of all positions xM for SM that violate the non-overlapping constraint. Right: two objects SR and -SM and their Minkowski sum, which coincides with the overlapping constraint (the negation of the non-overlapping constraint).

Fig. 2 .

 2 Fig. 2. Paving of an ellipsis.The interior red boxes belongs to I, the unknown blue boxes in the boundary belongs to B and the green outer boxes belongs to O.

Fig. 3 .

 3 Fig. 3. Sets whose Minkowski difference gives the overlapping constraint of two ellipsis, when rotation is taken into account. Left: the augmented reference ellipsis, with rotation coordinate set to 0. Right: the augmented moving ellipsis S ′ M .

Fig. 4 .

 4 Fig. 4. Steps of the inner Contractor.

Fig. 5 .

 5 Fig. 5. The sweep loop. The sequence of pictures illustrates a contraction for the lower bound of x1. At each step, the point x to be inflated is the lower-left corner of the gray box. The inner box [x] is painted in white. The lower bound of x1 can be reduced as soon as the projection of the white boxes on x2 spans the face [x2], which is the case at step e).

Fig. 7 .

 7 Fig. 7. Pavings obtained for cases 3-6. The precision is set to 3.25%. Left: with RSolver. Right: with our algorithm.

Fig. 8 .

 8 Fig.8. Time vs Precision (left: case 1; right: case 6). Each curve represents the paving time (vertical axis) with respect to the precision ε (horizontal axis). Both axis are in logarithmic scale. The blue curve is RSolver and the red curve is our method.

Fig. 9 .

 9 Fig. 9. 2D section of the 3D paving obtained for the case 7. From left to right: the rotation angle are 0, 0.4 and 0.7. Top: with Rsolver. Bottom: with our algorithm.

Fig. 10 .

 10 Fig. 10. 2D section of the 3D paving obtained for the case 8. From left to right, the rotation angle are 0, 0.1 and 0.3. This paving is calculated with our algorithm.

  ] ′′ );

	11
	12	push [x]1 and [x]2 in B;
	13	s ← s+surface([x]1)+surface([x]2);	// update the surface
	14	end
	15	end
	16	return (I, B, O);

Table 1 .

 1 Cases of study.

	Case Reference object Moving object Rotation
	1	1	1	no
	2	1	2	no
	3	1	3	no
	4	2	2	no
	5	2	3	no
	6	3	3	no
	7	1	1	yes
	8	3	3	yes

Table 2 .

 2 Running time (in s) for the 6 first cases (the precision is set to 3.25%).

	Case Rsolver Our algorithm
	1	4,07	0,37
	2	51,55	2,67
	3	611,85	7,90
	4	132,46	5,58
	5	656,11	12,00
	6	771,00	26,82

This subsolver is implemented with a standard branch & prune based on Hc4[START_REF] Benhamou | Continuous and interval constraints[END_REF]. Since only one solution is sought, at each node of the search, we check inequalities with a point p picked randomly in the current domain. If both are satisfied, the search is interrupted and p is returned. The depth of the search is also controlled by a precision on the domain width, which ensures that the subsolver terminates in bounded time. In case of normal termination, no p hence no inner box has been found.


Dante Kalise

Axel Kröner

Reduced-order minimum time control of advection-reaction-diffusion systems via dynamic programming

We present a numerical approach for the timeoptimal feedback control of an advection-reaction-diffusion model. Our approach is composed by three main building blocks: approximation of the abstract system dynamics, feedback computation based on dynamic programming and state observation. For the approximation of the abstract dynamics, we consider a finite element semi-discretization in space, leading to a large-scale dynamical system, whose dimension is reduced by means of a Balanced Truncation algorithm. Next, we apply the dynamic programming principle over the reduced model, and characterize the value function of the optimal control as a viscosity solution of a Hamilton-Jacobi-Bellman equation, which is numerically approximated with a semi-Lagrangian scheme. Finally, the computation of the corresponding feedback controls and its insertion into the control loop is performed by implementing a Luenberger observer.

I. INTRODUCTION

In this paper, we consider a minimum time problem with dynamics governed by a one-dimensional advection-reactiondiffusion model; such a problem arises, for instance, in minimum time stabilization of fluid flow passing through slender structures, or in the control of chemically reactive processes. The abstract evolutionary system is given by      ∂ t y = Ay + Bu, y(0) = x, z = Cy , (I. [START_REF]Morlab package for matlab[END_REF] where y = y(ξ, t) represents the state of the system (here ξ and t stands for the spatial coordinate and time, respectively), u = u(t) is a time-dependent scalar control signal, and z = z(ξ, t) represents an observation. The system is spatially confined to ξ ∈ Ω = (0, L) ⊂ R, L > 0, and time t ∈ [0, T ]. Using the usual notations for Lebesgue and Sobolev spaces, the free-dynamics operator A : H 1 0 (Ω) → H -1 (Ω) is defined as

Ay = a 1 ∂ 2 ξ y + a 2 ∂ ξ y + a 3 y,
with a 1 ∈ R + , a 2 , a 3 ∈ R, whereas the control operator B and the observation operator C are both indicator functions over open subsets Ω c ⊂ Ω and Ω o ⊂ Ω, respectively. For an initial state x ∈ H 1 0 (Ω) and u ∈ L ∞ (0, T ) the 1 D. Kalise is with the Johann Radon Institute for Computational and Applied Mathematics (RICAM), Austrian Academy of Sciences, Altenbergerstraße 69, A-4040 Linz, Austria. He acknowledges support from the FWF-START project: Sparse Approximation and Optimization in High Dimensions and the ITN -Marie Curie Grant n. 264735-SADCO. E-mail: dante.kalise@oeaw.ac.at 2 A. Kröner is with the Johann Radon Institute for Computational and Applied Mathematics (RICAM), Austrian Academy of Sciences, Altenbergerstraße 69, A-4040 Linz, Austria. E-mail: axel.kroener@oeaw.ac.at system has a unique solution in C((0, T ), C( Ω)), see [9, p. 360] and [21, p.19]. Furthermore, the uncontrolled system is exponentially stable for a 3 < π 2 a 1 . For a target T = { z ∈ C( Ω) | z C( Ω) ≤ ε }, ε > 0, we address the following minimum time problem:

           min u∈U T, s.t. ∂ t y = Ay + Bu , y(0) = x, y(T ) ∈ T (I.2) for U = L ∞ ((0, T ), U ) with U = { u ∈ R | u a ≤ u(s) ≤ u b a.e. } ,
and bounds u a , u b ∈ R. For further comments on the continuous problem for the case a 2 = 0, we refer to [START_REF] Kunisch | Time optimal control of the heat equation with pointwise control constraints[END_REF].

Our aim is to derive approximate controllers in feedback form for the aforementioned problem. A first step towards this goal is to perform a finite element semi-discretization in space of the advection-reaction-diffusion equation. The resulting control problem replaces the abstract dynamics (I.1) by a large-scale, finite dimensional dynamical system. At this point, it would be theoretically feasible to apply the Dynamic Programming Principle (DPP) in order to characterize the value function of the optimal control problem as the viscosity solution of a Hamilton-Jacobi-Bellman (HJB) equation, as presented in [START_REF] Ferretti | Internal approximation schemes for optimal control problems in Hilbert spaces[END_REF]. However, the fact that such an equation is defined over a state space of the same dimension as the semi-discretized dynamics, renders it computationally intractable. Therefore, it is necessary to include an additional approximation step in order to reduce the dimensionality of the semi-discrete system. A standard tool for this purpose is the application of a Balanced Truncation algorithm, which synthesizes the relevant information in terms of controllability and observability of the system, into a lower dimensional model of prescribed accuracy (with respect to a measure of the input-output mapping). Once the reduced dynamics have been obtained, the application of the DPP over the resulting system leads to an approximate, low-order feedback controller. However, the control law assumes a full knowledge of the reduced state of the system which, according to (I.1), is not available; instead, we receive a measured output z. Therefore, the final step in the control loop consists in the implementation of a reduced order Luenberger observer in order to generate a reliable estimate of the internal dynamics of the system, for an accurate computation of the feedback mapping.

To set this paper in perspective, we relate our approach to the previous works [START_REF] Alla | An adaptive POD approximation method for the control of advection-diffusion equations[END_REF], [START_REF] Kunisch | HJB-POD based feedback design for the optimal control of evolution problems[END_REF], [START_REF] Kunisch | POD-based feedback control of the Burgers equation by solving the evolutionary HJB equation[END_REF], [START_REF] Kröner | Optimal feedback control for the undamped wave equation by solving a HJB equation[END_REF], which also consider optimal feedback control problems for infinite-dimensional dynamics, using either proper orthogonal decomposition or spectral elements to obtain a low dimensional semidiscrete system. A numerical solution for feedback control problems for nonlinear parabolic equations is considered in [START_REF] Benner | Numerical solution of optimal control problems for parabolic systems[END_REF]. Numerical implementations and approximation results for feedback problems of (second-order) hyperbolic equations using Riccati equations can be found in [START_REF] Hernández | Numerical approximation of the LQR problem in a strongly damped wave equation[END_REF], [START_REF] Hernández | A locking-free scheme for the LQR control of a Timoshenko beam[END_REF], and references therein. Similar approaches, based on model predictive control, which can be interpreted as a relaxed version of dynamic programming, have been presented in [START_REF] Dubljevic | Predictive control of parabolic pdes with state and control constraints[END_REF], [START_REF] Altmüller | Distributed and boundary model predictive control for the heat equation[END_REF]. Minimum time problems were first considered in [START_REF] Bellman | On the bang-bang control problem[END_REF]. For time optimal control problems of parabolic equations see [START_REF] Fattorini | Time and norm optimal controls: A survey of recent results and open problems[END_REF], [START_REF] Kunisch | Time optimal control of the heat equation with pointwise control constraints[END_REF], [START_REF] Lasiecka | Finite element approximation of time optimal control problems for parabolic equations with Dirichelt boundary conditions, System Modeling and Optimization[END_REF], [START_REF] Raymond | Pontryagin's principle for time-optimal problems[END_REF], [START_REF] Micu | Time optimal boundary controls for the heat equation[END_REF], [START_REF] Zuazua | Approximate controllability for semilinear heat equations with globally Lipschitz nonlinearities[END_REF]. The novelty of this article resides in the combination of HJB techniques for the computation of time-optimal feedback controllers with model reduction and state observation algorithms providing thus a consistent approach for closed-loop control of infinitedimensional systems. The paper is organized as follows. In Section II , the semidiscretization of the dynamics and its reduction via Balanced Truncation is presented. In Section III we introduce the dynamic programming framework and a semi-Lagrangian scheme for the solution of the resulting HJB equation. In Section IV, we propose a Luenberger observer for state estimation in the closed-loop implementation. Section V reports numerical examples assessing the performance and capabilities of the proposed approach.

II. DISCRETIZATION AND MODEL REDUCTION

In this section, we present a semi-discretization of the reaction-diffusion-advection equation (I.1) resulting in a finite-dimensional dynamical system. At this level, the application of a HJB-based approach is still unfeasible from a computational point of view, since the PDE characterizing the value function of the problem is posed over a state space of the same dimension as the semi-discretization. To reduce the dimensionality of the problem, we apply a Balanced Truncation algorithm to the semi-discrete dynamical system. In this way, we obtain a low-dimensional model for which the corresponding HJB equation is numerically solvable. The building blocks of this procedure are illustrated in Figure 1. Balanced Truncation for optimal control of evolution equations has been considered in [START_REF] Reyes | A balanced truncation-based strategy for optimal control of evolution problems[END_REF]. Classical Balanced Truncation is limited to linear problems, but it has desirable properties, such as stability preservation. To discretize (I.1) in space we use continuous, piecewise linear finite elements. For a given

N ∈ N, let G = { j∆x } N +1 j=0 , ∆x = L N +1 , be a regular subdivision of the in- terval [0, L]. For the corresponding standard nodal basis func- tions { ϕ i } N i=1 , we introduce the mass matrix M ∈ R N ×N with M i,j = (ϕ i , ϕ j ), the stiffness matrix K ∈ R N ×N with K i,j = (∂ ξ ϕ i , ∂ ξ ϕ j )
, and for the discretization of the advection term D ∈ R N ×N with D i,j = (∂ ξ ϕ i , ϕ j ), where (•, •) denotes the usual inner product in L 2 (Ω). Further, we define the observation matrix

C N = M -1 No χ No,N M ∈ R No×N , N o ≤ N , for matrix χ No,N ∈ R No×N of rank N o with entries χ No,N i,j
∈ { 0, 1 }, according to the observation set Ω o and the mass matrix M No ∈ R No×No which we choose according to the observed components. For the ansatz y ∆x = N i=1 y i ϕ i , where y i are time-dependent nodal values of the finite element functions y ∆x , we introduce y N = (y 1 , . . . , y N ) T which satisfies the semi-discrete system given by

     ∂ t y N = A N y N + B N u, y N (0) = x N , z = C N y N (II.1) for x N ∈ R N with x N i = (x, ϕ i ), and
A N = M -1 (a 1 K + a 2 D + a 3 M ) , B N = M -1 χ N c , [χ N c ] i = (χ c , ϕ i ) , with χ c the characteristic function of Ω c .
Error estimates in L 2 -norm for this type of semi-discrete approximation are presented in [START_REF] Wheeler | A priori L 2 error estimates for Galerkin approximations to parabolic partial differential equations[END_REF].

For the convenience of the reader, we shortly recall the basic idea of Balanced Truncation, see [START_REF] Volkwein | Proper Orthogonal Decomposition: Applications in Optimization and Control[END_REF], [START_REF] Reyes | A balanced truncation-based strategy for optimal control of evolution problems[END_REF]. In the following we assume that A N is stable, (A N , B N ) is controllable, and (A N , C N ) is observable. To simplify the notation we drop the index N . First, note that under the aforementioned assumptions, there exists a controllability Gramian P and an observability Gramian Q solving the Lyapunov equations

PA + A T P = -C T C, QA T + AQ = -BB T .
Since the Gramians are positive-definite, We can compute their Cholesky factors R and L, i.e.

P = RR T , Q = LL T ,
and solve the singular value decomposition of L T R. We fix a reduced-order model dimension r ∈ N, and split

L T R = (U 1 U 2 ) Σ 1 0 0 Σ 2 V 1 V 2 , with orthonormal matrices U = (U 1 U 2 ) and V = (V 1 V 2 ) T ,
and diagonal matrices

Σ 1 = diag(σ 1 , . . . , σ r ), Σ 2 = diag(σ r+1 , . . . , σ l ) , with σ 1 ≥ • • • ≥ σ r ≫ σ r+1 ≥ • • • ≥ σ l > 0 ,
and l = rank(L T R). The singular values of L T R, known in this context as Hankel singular values, provide a simultaneous measure of both controllability and observability energies of the corresponding balanced state. By setting

W = LU 1 Σ -1 2 1 , S = RV 1 Σ -1 2 1
we can compute the following transformation

A r = W T AS, B r = S T B, C r = CW ,
and obtain the reduced problem

     ∂ t y r = A r y r + B r u, y r (0) = W T x, z r = C r y r (II.2) with y r ∈ R r .
For the error between the solution of the reduced model (II.2) in comparison to the solution of the full model (II.1) on (0, ∞) there holds, see, e.g., [START_REF] Reyes | A balanced truncation-based strategy for optimal control of evolution problems[END_REF],

z r -z L 2 (0,∞) ≤ 2(σ r+1 + • • • + σ l ) u L 2 (0,∞) .
This estimate gives an upper bound for the error in the observation depending on the dimension of the reduced model.

III. THE DYNAMIC PROGRAMMING APPROACH

By applying a Balanced Truncation algorithm over the semi-discrete model, we obtain a reduced dynamical system. This leads to the following approximated minimum time problem

           min u∈U T, s.t. ∂ t y r = A r y + B r u , y r (0) = W T x, y r (T ) ∈ T , (III.1)
This latter problem is solved by applying the dynamic programming principle, which in this case is given by

T (x) = inf u∈U (τ + T (y x (τ, u))) (III.2) for all x ∈ R = { x ∈ R l | T (x) < ∞ }, τ ∈ [0, T (x)) with x ∈ T .
Here y x (τ, u) denotes the state at time τ for initial state x at time zero and control u. The DPP leads to the following HJB equation characterizing the minimum time

     H(x, DT ) = 0 in R \ T , T = 0 on ∂T , T (x) → +∞ as x → x 0 ∈ ∂R (III.3) with Hamiltonian H(x, p) = sup u∈U (-f (x, u) T p) -1 , and f (x, u) = A r x + B r u .
Applying the Kruzkov transform

v(x) = 1 -e -T (x) x ∈ R 1 x ∈ R we further obtain v + H(x, Dv) = 0 in R n \ T , v = 0 on ∂T . (III.4)
In order to solve the HJB equation numerically, we apply a semi-Lagrangian scheme following [START_REF] Alla | An efficient policy iteration algorithm for dynamic programming equations[END_REF]; for a general introduction to semi-Lagrangian schemes we refer to [START_REF] Bardi | Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations[END_REF], [START_REF] Falcone | Semi-Lagrangian Approximation Schemes for Linear and Hamilton-Jacobi Equations[END_REF]. Be begin by discretizing the dynamical system in time with stepsize h, and apply the DPP for the discrete-time dynamical system. Then, for spatial mesh parameter k ∈ R + , we introduce a regular mesh G = { x I | I ∈ Z r , Ik ∈ S } and denote the set of all multi-indices I with x I = Ik ∈ G by J. We restrict the computations to a domain S ⊂ R r , and impose an artificial Dirichlet boundary condition on ∂S which we set to the value 1. The fully discretized HJB equation then reads

                     v h,k (x I ) = min u∈U (βI[v h,k ](x I + hf (x I , u))) + 1 -β, for I ∈ J in , v h,k (x I ) = 0 for I ∈ J T , v h,k (x I ) = 1
for I ∈ J out (III.5) with β = e -h and (III.7)

J out = { I | x I + hf (x I , u) ∈ S for any u ∈ U } , J T = { I | x I ∈ T ∩ S } , J in = I ∈ J \ (J out ∪ J T ).
This latter expression requires the knowledge of the whole state, whereas our system considers an observation equation.

In the next section, we introduce an observer implementation for an accurate trajectory computation.

IV. LUENBERGER OBSERVER In order to generate a coherent link between the observation and the feedback law, it is necessary to implement a state observer. More specifically, the approximation of system (I.1) generates a dynamical output z, whereas the feedback mapping (III.7) assumes that a full knowledge of the internal state y is available. By implementing a Luenberger observer, we are able to circumvent this difficulty and to close the control loop. The interplay of the observer and control blocks is illustrated in Figure 2. The observer is built upon the reduced order dynamics, corresponding to the "true" state from which the feedback mapping can be computed. The estimated state ŷ is governed by

∂ t ŷ = A r ŷ + B r u + L(C r ŷ -z) , ŷ(0) = ŷ0 ,
where ŷ0 ∈ R r and L ∈ R r×No needs to be computed such that A r -LC r is asymptotically stable and with a decay rate faster than the free system dynamics A r . In our particular setting, we recall that A r is asymptotically stable, and therefore our problem reduces to find a suitable L which accelerates the decay of the estimation error. The observer gives us now a corresponding state estimate ŷ which is inserted in the expression (III.7), yielding a feedback u = u(ŷ) to be connected to the full-order system. Plant (real system) Initial state HJB equation Observer We illustrate the application of the proposed approach and the interplay between the different building blocks. First, we describe specific details concerning common settings for the numerical implementation. System configuration and parameters. In both tests, the physical domain corresponds to Ω = [0, 1], the operator B is an indicator function over the subset Ω c = (0.4, 0.6), and the observation operator C is an indicator function over Ω o = (0.3, 0.7) for Test 1 and Ω o = Ω for Test 2. In general, the main issue when prescribing a setting relates to abstract controllability and observability assumptions, which play a role in both the reduction and in the HJB step; the aforementioned configuration has been selected in order to fulfill these requirements. Nevertheless, we report that the presented approach can be implemented also for point control/observation configurations. The control set corresponds to U = [-1, 1]. Simulations are shown with the initial condition x(ξ) = 0.2 sin(πξ), and for a time frame of 1[s]. Semi-discretization and model reduction. Semidiscretization in space is performed with piecewise linear finite elements. In both tests, the results were obtained with a discretization of 200 elements in space. The resulting large-scale system is reduced by means of a Balanced Truncation algorithm; we use the numerical implementation provided by the MORLAB package [START_REF]Morlab package for matlab[END_REF]. Hamilton-Jacobi-Bellman equation. The numerical approximation of the resulting HJB equation by means of a semi-Lagrangian schemes follows the general guidelines presented in the literature as in [START_REF] Bardi | Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations[END_REF], [START_REF] Falcone | Semi-Lagrangian Approximation Schemes for Linear and Hamilton-Jacobi Equations[END_REF]. The equation is solved in a subdomain of the reduced state-space, in our case the computational domain is given by Q = [-1, 1] 3 considering only three reduced states in both tests. The statespace discretization is given by the space parameter k = 0.025 and the pseudo-time parameter h = k/ A r , which gives a good balance between accuracy and computation time (as arrival points are computed in close neighborhood of their departure). The semi-Lagrangian discretization leads to the system (III.5) which is solved by means of a value iteration algorithm defined upon the grid point values V = {v h,k (x I )},

V n+1 = S(V n ) [S(V )] I =                      min u∈U (βI[V ](x I + hf (x I , u))) + 1 -β, for I ∈ J in , 0 for I ∈ J T , 1 for I ∈ J out which in our case is stopped when two consecutive iterations hold V n -V n+1 ≤ k 2 .
The control set U is discretized into 10 equidistant points. The target set T is a ball around the origin of radius ǫ = h. State observer and trajectory computation. The implementation of the observer, in this particular case, reduces to determining L ∈ R r×r such that the eigenvalues of A r -LC r decay sufficiently fast compared to the system dynamics. For this problem, we prescribe eigenvalues 10 times faster than the slowest eigenvalues of A r , and the initial observation estimate is set ŷ(0) = 0. For trajectory computation, at every time step, once the current state of the system has been estimated, the feedback rule is computed by means of (III.7), and the system advances in time with a second-order implicit integrator. Let us stress that, because of the computational effort performed in the HJB step, where the value function is computed for the whole state space of interest, the cpu time related to trajectory and feedback computation is negligible, as it mainly reduces to a single-node evaluation of arrival points.

We present two numerical tests, first a purely reactiondiffusion model, and in a further step we also include advective effects.

A. Test 1: reaction-diffusion

In this first test we set for the free-dynamics operator a 1 = 0.05, a 2 = 0, and a 3 = 0.1, i.e., advection is not considered. The resulting reaction-diffusion system is exponentially stable, and thus our control goal reduces to achieve faster regulation of the initial state to the origin. After discretizing the system in space by finite elements we balance the resulting large-scale system. Figure 3 shows the evolution of the Hankel singular values in balanced coordinates; the error of the truncation step is governed by the largest Hankel singular value which has been neglected, and in our case, to consider 3 reduced states leads to an acceptable error. Note that this figure also assess our choice of Balanced Truncation for the reduction procedure, as the fast decay of the singular values is fundamental if a good approximation with a few number of states is required. For a system with 3 reduced states, computation time of the HJB step is of approximately 50 seconds. CPU time increases dramatically for high-dimensional problems, however, by using accelerated iterative algorithms as in [START_REF] Alla | An efficient policy iteration algorithm for dynamic programming equations[END_REF], it is possible to solve similar problems in four dimensions in the same amount of time, and to keep a cpu time for five dimensional problems on the order of 10 minutes (for coarse meshes). The next step is illustrated in Figure 4, where the value function obtained via the DPP approach is shown. Such value function shape often appears in minimum time control of stable linear systems, and is related to switching curves splitting the associated control space into a discrete number of sets, i.e., a bang-bang type of control. Dynamical response of the system is presented in Figures 5 and6, where the performance of the minimum time controller can be observed compared to the free dynamics, the fast decay of the estimation error, and the respective control signal, which emulates the expected bang-bang behavior. Differences with an exact bang-bang controller, namely the chattering after the first switching, are observed due to the approximate character of the discrete HJB equation and the computation of trajectories over a discrete state-space grid (where switching curves are also approximated), which is a natural limitation of discrete dynamic programming-based approaches. Aiming at a more realistic online implementation, in [START_REF] Falcone | A semi-lagrangian scheme for Lp-penalized minimum time problems[END_REF], the authors have considered an ℓ-1 penalized minimum time problems where the bang-bang structure is replaced by a bang-zero-bang behavior which reduces chat- tering in a considerable way. The controller exhibits a good behavior with respect to output noise (we consider a white noise amplified up to a 30% of the maximum output value), despite the assignment of fast poles for the observer. This is due to the bang-bang structure of the system, where only strong enough output perturbations will cause an estimate to be located on the other side of a given switching surface.

B. Test 2: advection-diffusion-reaction

We perform a second test including advective effects; we set a 1 = 0.05, a 2 = 1.0, and a 3 = 0.1. Note that despite the inclusion of an advective term, such a setting does not deteriorate the stability properties of the original continuous system. A first difference with respect to the previous test can be observed in Figure 7, where it can be seen that the scaling of the decay of the Hankel singular values change, thus suggesting that a larger number of reduced states should be considered for a similar level of accuracy as in the previous example. However, in this test, we will consider three reduced states, as the error associated to the fourth singular value is still within an acceptable range.

The approximation of the associated HJB equation is illustrated in Figure 8, where a cut of the value function along x 3 = 0 is shown. The dynamic response of the controlled system is presented in Figures 9 and10. In Figure 9, it can be seen that although the presence of advection translates into a regulation of the state in finite time, there is a considerable improvement in the performance of the system when a minimum time controller is considered, as both advection and regulation via the control action are combined. In this test case we assume full observation of the state, i.e. our estimation error is uniquely related to the guess of the initial condition ŷ(0). Note that, although there is a decay in the estimation error, the effect of advection can be clearly observed in Figure 10. Finally, in a similar way as in the Top: Free-system dynamics (position v/s time). Bottom: Full-order controlled dynamics with reduced-order, minimum time controller.

previous example, Figure 10 shows an approximate bangbang control signal, with inaccuracies due to the discrete resolution of the HJB equation.

VI. CONCLUDING REMARKS

We have presented a computationally feasible framework for the implementation of approximate feedback controls in infinite-dimensional systems. Our proposed approach consists in three main steps: semi-discretization and model order reduction of the abstract system dynamics, the application of DPP techniques in order to achieve a characterization of the value function for the lower dimensional optimal control problem, and the implementation of a state observer in order to close the control loop. Every block has a solid computational framework where different numerical methods are available. Our results suggest that the approach is able to yield reasonably robust controllers (in the sense that they can perform well when connected to the large-scale, semidiscrete dynamics), already with a reduced number of states. The proposed approach is flexible as it supports variations in the application of reduction techniques, the optimal feedback to be sought, the type of observer considered, and therefore it can be used for a wide class of dynamics and control design problems. Future extensions of this work shall address the formulation of a precise error quantification connecting the different approximation errors introduced in every step, as also a precise characterization of the robustness properties observed in the numerical essays.

Fig. 1 :

 1 Fig.1: Building blocks for the low-order optimal control problem.

(III. 6 )

 6 Here I[•] denotes a linear interpolation operator. Having approximated the value function by means of the aforementioned scheme, we can recover the feedback controller for a given state x by u(x) = arg min u∈U βI[v h,k ](x + hf (x, u)) .

Fig. 2 :

 2 Fig. 2: Overview over the observer mechanism

Fig. 3 :

 3 Fig. 3: Test 1 (without advection). Decay of the Hankel singular values in balanced coordinates.

Fig. 4 :

 4 Fig. 4: Test 1 (without advection). Value function: cut at x3 = 0.

Fig. 6 :

 6 Fig.6: Test 1 (without advection). Dynamical response of the system. Dynamic discrepancy between internal state of the system and its estimate, and control signal. Tests without noise and with white noise in the output of the plant (top and bottom, respectively).

Fig. 7 :

 7 Fig. 7: Test 2 (with advection). Decay of the Hankel singular values in balanced coordinates.

Fig. 8 :

 8 Fig. 8: Test 2 (with advection). Value function: cut at x3 = 0.

Fig. 9 :

 9 Fig.9: Test 2 (with advection). Dynamical response of the system. Top: Free-system dynamics (position v/s time). Bottom: Full-order controlled dynamics with reduced-order, minimum time controller.

Fig. 10 :

 10 Fig.10: Test 2 (with advection). Dynamical response of the system. Top: Dynamic discrepancy between internal state of the system and its estimate. Bottom: Control signal.