Antoine Laborde Agroparistech

Benoît Jaillais

Jean-Michel Roger

Maxime Metz

Delphine Jouan

Rimbaud Bouveresse

Luc Eveleigh

Christophe Cordella

Subpixel detection of peanut in wheat flour using a matched subspace detector algorithm and near-infrared hyperspectral imaging

Keywords: Hyperspectral imaging, near-infrared spectroscopy, detection algorithm, spectral simulation, matched subspace detector Hyperspectral Imaging NIR: Near-Infrared LMM: Linear Mixing Model MSD: Matched Subspace Detector PCA: Principal Component Analysis

The detection of adulterations in food powder product represents a high interest especially when it concerns the health of the consumers. The food industry is concerned by peanut adulteration since it is a major food allergen and it is often used in transformed food products. Near-infrared hyperspectral imaging is an emerging technology for food inspection. It was used in this work to detect peanut flour adulteration in wheat flour. The detection of peanut particles was challenging for two main reasons: the particle size is smaller than the pixel size leading to impure spectral profiles; peanut and wheat flour exhibit similar spectral signature and variability. A Matched Subspace Detector (MSD) was designed to take into account those difficulties and detect peanut adulteration at the pixel scale using the associated spectrum. Defatted peanut flour and wheat flour were mixed in eight different proportions (from 0.02% to 20%) to test the performances of the hyperspectral measurement technique with the detection algorithm. Additionally, a set of simulated data was generated to overcome the lack of reference value at the pixel scale. The application of MSD on simulated data showed to be relevant with experimental observations and helped to design the detector. The most performant designs for the detector were compared using the simulation by estimating the sensitivity as well as the experimental data by comparing the detection maps. Finally, the number and the position of detections on experimental data were investigated to show the relevancy of the results. This proved that the use of hyperspectral imaging and a fine-tuned MSD enables to detect a global adulteration of 0.2% of peanut in wheat flour using a single hyperspectral image acquisition.

Introduction

Hyperspectral Imaging (HSI) is a technique combining spectroscopy with spatial imaging to obtain both spatial and spectral information from a sample. When associated with the nearinfrared (NIR) spectral range, HSI is a powerful technique to provide a fast, non-destructive and cost-effective control method. Since many samples may be chemically heterogeneous, the spatial information provides a high added value for many application fields. Indeed, HSI has been investigated for several decades in remote sensing for earth observation [START_REF] Xie | Remote sensing imagery in vegetation mapping : a review[END_REF], food industry [START_REF] Elmasry | Principles and Applications of Hyperspectral Imaging in Quality Evaluation of Agro-Food Products : A Review[END_REF], agriculture [START_REF] Dale | Hyperspectral imaging applications in agriculture and agro-food product quality and safety control: A review[END_REF] or medical uses [START_REF] Lu | Medical hyperspectral imaging : a review[END_REF]. Particularly, HSI is an emerging technology for food inspection since it provides non-destructive analysis of heterogeneous samples [START_REF] Feng | Application of Hyperspectral Imaging in Food Safety Inspection and Control : A Review[END_REF] and hyperspectral images allow the visualization of chemical maps of the samples. As an example, Elmasry et al. use this technology to provide water, fat and protein distributions on beef samples [START_REF] Elmasry | Chemical-free assessment and mapping of major constituents in beef using hyperspectral imaging[END_REF]. The study of the estimation of food nutrients in samples lead to the assessment of the quality of such products [START_REF] Gowen | Hyperspectral imagingan emerging process analytical tool for food quality and safety control[END_REF].

Hyperspectral imaging technique enables to obtain a spectral measurement for each individual detector of the sensor. As a consequence, each spectrum is representative of a small surface of the sample. Spatial imaging makes the technology more sensitive to minor components since they have more influence on the field of view of one pixel than on the entire sample. This offers the opportunity to detect adulterants in food by characterizing each pixel of the image [START_REF] Kamruzzaman | Assessment of Visible Near-Infrared Hyperspectral Imaging as a Tool for Detection of Horsemeat Adulteration in Minced Beef[END_REF], [START_REF] Kamruzzaman | Rapid and non-destructive detection of chicken adulteration in minced beef using visible near-infrared hyperspectral imaging and machine learning[END_REF], [START_REF] Su | Fourier Transform Infrared and Raman and Hyperspectral Imaging Techniques for Quality Determinations of Powdery Foods: A Review[END_REF]. In practice, it consists of classifying each pixel spectrum of the hyperspectral image as a target or a background spectral signature. As NIR spectroscopy has shown to be a powerful technique for characterizing organic matter [START_REF] Ellis | Fingerprinting food: current technologies for the detection of food adulteration and contamination[END_REF], [START_REF] Manning | Developing systems to control food adulteration[END_REF], HSI appears to be a promising tool for adulteration detection in the food industry. The literature shows plenty of such applications: Vermeulen et al. studied the detection of ergot bodies in cereal flour [START_REF] Vermeulen | Online detection and quantification of particles of ergot bodies in cereal flour using nearinfrared hyperspectral imaging[END_REF], Fernández Pierna et al. investigated the detection of melamine in milk powder [START_REF] Fernández Pierna | Use of a multivariate moving window PCA for the untargeted detection of contaminants in agro-food products, as exemplified by the detection of melamine levels in milk using vibrational spectroscopy[END_REF] [15], Verdú et al. proposed to study the adulteration of wheat products [START_REF] Verdú | Detection of adulterations with different grains in wheat products based on the hyperspectral image technique: The specific cases of flour and bread[END_REF]. Some adulteration cases are linked to serious public health problems like allergic reactions.

The case of peanut allergy is interesting as a typical example. Whereas this allergy is a major issue worldwide [START_REF] Burks | Peanut allergy[END_REF] [START_REF] Hefle | Food allergy and the food industry[END_REF], peanuts are still widely used in the food industry [START_REF] Zhao | Potential use of peanut by-products in food processing: A review[END_REF] and maybe also linked to adulteration cases. HSI seems to be a relevant technology to tackle the problem of adulteration in the food industry as a detection tool. Such detection has already been made using crushed peanut in wheat flour [START_REF] Mishra | Detection and quantification of peanut traces in wheat flour by near infrared hyperspectral imaging spectroscopy using principal-component analysis[END_REF] [START_REF] Mishra | Application of independent components analysis with the JADE algorithm and NIR hyperspectral imaging for revealing food adulteration[END_REF]. Defatted peanut flour contains smaller peanut particles without fatty acid components. Detecting such a product in wheat flour may be of high interest since powder samples are more frequently used in the food industry.

However, the use of HSI for detection may be a challenge for several reasons. Indeed, the purpose is to identify different materials based on their spectral signature to label each pixel as a target or background pixel. However, there is not a unique spectral signature for each material for two main reasons [START_REF] Manolakis | Detection algorithms for hyperspectral imaging applications[END_REF]: the reflectance value at each wavelength of a given material is not deterministic but is a random variable. Its variability is linked to the lighting conditions, the material surface, the sample heterogeneity, and many other factors.

Moreover, two different materials may have very similar spectra. Especially for food industry products, the shape of NIR spectra is often similar since they are the result of the mix of the main nutrients. This is particularly the case for wheat and peanut as they are similar products once transformed into flour. More specifically, as peanuts are defatted, the fatty acid NIR fingerprint does not appear in the spectrum. Thus, the ambiguity of spectral information and the spectral variability issue are two main challenges for detection purposes in the food industry.

Another difficulty arises when dealing with powder samples. Indeed, the particle size may be smaller than the pixel size. For a hyperspectral camera, a pixel integrates the radiance signal from all materials in its field of view. Thus, if this field of view contains particles with different spectral signatures, the pixel is considered as mixed and the resulting spectrum does not correspond to any pure chemical defined as target or background. This problem is known as the subpixel detection [START_REF] Manolakis | Detection algorithms for hyperspectral imaging applications[END_REF].

The detection problem can be seen as a binary hypothesis test for which each pixel is assigned to the target or the background class. In this condition, no spatial pattern is taken into account and each pixel is considered as an independent situation. For the null hypothesis, the pixel is assumed to contain only background particles whereas, for the alternative hypothesis, the pixel contains some target particles. Detection algorithms have been designed to take into account both variability and subpixel issues using spectral modeling [START_REF] Manolakis | Detection algorithms for hyperspectral imaging applications[END_REF].

The mixed pixel issue is tackled using the Linear Mixing Model (LMM) [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regressionbased approaches[END_REF] which assumes the radiance measured by a pixel is the sum of the chemical radiances weighted by their surface contribution in the pixel field of view. The spectral variability issue can be addressed using subspace modeling that is integrated into the LMM formulation. Finally, a detector can be designed using the likelihood ratio and by comparing its value to a threshold. Particularly, the Matched Subspace Detector (MSD) [START_REF] Manolakis | Detection algorithms for hyperspectral imaging applications[END_REF] [24] is derived using the assumption that both target and background variabilities are modeled. Such a detector requires to be fine-designed to give high performances. In the case of MSD, the way the subspaces are designed as well as their dimensions highly influence its performances. Such detectors are often evaluated on reference hyperspectral images that are already labeled at the pixel scale (example of Hydice data set for remote sensing applications [START_REF] Manolakis | Detection algorithms for hyperspectral imaging applications[END_REF]). Thus, detector performances can be calculated on reference images so that the design can be optimized. For the inspection of food samples, the field of view of a pixel is much smaller (1 mm × 1 mm) [START_REF] Geladi | Hyperspectral imaging : calibration problems and solutions[END_REF] than for remote sensing applications. In this situation, having chemical reference information for each pixel would require performing reference analysis for each pixel field of view surface. This process is not realistic when dealing with food samples. As a consequence, no reference data can be used at the pixel scale for designing the detector.

Spectral data simulation consists of generating new spectral data from a known statistical distribution. This technique can be used to compare different modeling techniques [START_REF] Jörgensen | A comparison of methods for analysing regression models with both spectral and designed variables[END_REF]. Such a simulation can take the spectral variability of a pure sample into account by using Principal Component Analysis (PCA) [START_REF] Mishra | Detection and quantification of peanut traces in wheat flour by near infrared hyperspectral imaging spectroscopy using principal-component analysis[END_REF]. Additionally, the LMM may be used to perform mixed spectra resulting from an interaction between several materials. The use of a near-infrared hyperspectral imaging system for the detection of adulterant in food has been studied in plenty of applications. However, the issues regarding the variability of the samples, the spectral ambiguity between species, the mixed pixels and the lack of reference data for the detector design make the detection difficult. To our knowledge, no study has been proposed to tackle detection for such samples with particle size smaller than the pixel. The purpose of this study is to evaluate how the MSD approach using the LMM and the modeling of spectral variability can provide performant detection for such a detection problem. As no reference values for the detector design are available, a spectral simulation method is proposed. After studying the performances of the detector on simulated data, we propose to study the detection using real hyperspectral measurements of flour mix at graduated concentrations of peanut.

Materials and methods

Samples

White wheat flour (Grands Moulins de Paris, Francine batch number ER510 -FTU104, France) was used as the background sample. Samples were taken from two different packs.

Defatted peanut flour (KoRo Handels GmbH, batch number C170151, Germany) was used for the target sample. Flour samples were mixed together to obtain 8 different mass concentrations of peanut flour: 20 %, 10 %, 5 %, 2 %, 1 %, 0.5 %, 0.2 % and 0.02 % for a total mass of 13.75 g. Mass measurements were performed using a precision balance (Sartorius Entris, 0.01 g precision). Additionally, pure peanut flour and pure wheat flour were prepared. Each sample was prepared in triplicate. Mixed samples were put in a container to be shaken and mixed with a spatula. Samples were put in a rectangular sample holder (30 mm width x 70 mm length) made of a 7 mm depth cavity. The top of the sample holder was skimmed to remove excess powder without affecting the packing density.

Hyperspectral imaging system

A line-scan pushbroom Specim SWIR camera (SPECIM, Oulu, Finland) was used for the image acquisition. The hyperspectral camera acquired 288 spectral bands from 900 nm to 2500 nm with a 5.6 nm step. The camera acquired 392 pixels per line and the pixel size was 250 µm × 250 µm. Six halogen lamps were used for the measurement and heated up for 30 minutes before the acquisition. A white reference measurement was performed before each acquisition using a white diffuse reflectance standard (Spectralon®, SRS-99-010, Labsphere).

Additionally, the dark reference image was acquired after closing the shutter of the camera.

Each sample was measured independently leading to 30 data cubes.

Data processing

Each image was cropped to focus on the sample in the central cavity of the sample holder leading to data cubes of size 200 × 320 × 188. The white reference image was averaged along the perpendicular direction of the sensor array to obtain one spectrum for every pixel of the sensor line ( ).The reflectance image is calculated using (Eq. 1):

(Eq. 1)
where corresponds to the black measurement for each pixel and is the raw intensity measurement of the sample. First (under 1200 nm) and last (over 2200 nm) wavelengths were removed as spectra were too noisy. Spectra were processed using a Savitsky-Golay filter to reduce the noise for the remaining wavelengths (2 nd order polynomial, 7-points window, and no derivative). A Standard Normal Variate (SNV) transformation was applied to compensate for scattering effects.

Spectral simulation using Principal Component Analysis

A spectrum can be considered as a vector lying in a -dimensional space. Each dimension is defined by one variable, namely a wavelength of the spectral range. The variability of spectral data is the variance of the reflectance for each wavelength in this -dimensional space. NIR spectral data exhibits a high correlation between the variables. As a consequence, defining the spectral variability independently for each variable is not efficient.

PCA is a method for dimensionality reduction that decomposes the data matrix according to orthogonal sources of the highest possible variance. The data matrix can be decomposed as follow:

(Eq. 2)
where is the score matrix, is the loading matrix, is the matrix of residuals, and the upper script symbol T refers to the transposed matrix.

Under this representation, the distribution of scores for each component can be considered as Gaussian with mean -which is, in practice, null -and variance . Each component of the matrix provides a unit vector spanning the subspaces of the sample variability. As a consequence, a new spectrum can be simulated by generating its coordinates in the principal components space. In this representation, the coordinates are generated using Gaussian distributions.

The Linear Mixing Model (LMM) describes the linear interaction between spectral signatures of pure materials in a mixture context [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regressionbased approaches[END_REF]. In this model, a spectrum is described by (Eq.

3).

(Eq. 3)

where is the pure spectrum of the i-th material, is the associated concentration, is the residual vector and is the number of species in the model.

The PCA and LMM are combined considering that may represent a source of variability as well as a pure spectral signature. The process of data simulation used in this study is described by the following procedure (Figure 1):

1) PCA is performed on the data matrix of pure wheat data and pure peanut data distinctly. The average spectra of both matrices and are calculated and considered as the pure spectral signatures of the materials.

2) For every dimension , the distribution of scores is assumed to be Gaussian of mean (which is equal to 0) and standard deviation . These parameters are estimated on the pure sample image for wheat and peanut distinctly.

3) For a given peanut percentage , the average spectrum is simulated using the LMM:

(Eq. 4)

being a scalar varying between 0 and 1.

4) The simulation is completed by adding the random variability to the average simulated spectrum . The variability is obtained by multiplying the simulated scores with the principal component loadings. The simulated scores are generated from Gaussian distributions with the parameters calculated in step 2) for peanut and wheat separately. The total variability attributed to is a balance between peanut and wheat controlled by the proportion .

(Eq. 5)

where the tilde symbol designates the simulated matrices.

[Insert Figure 1] The simulation procedure is applied to obtain 100 spectra for each peanut concentration 5%, 10%, 15% and 20%.

Matched subspace detector design

According to the LMM, a given pixel has a spectrum described by: (Eq. 6)

where the vectors define the variability space of . These vectors can be obtained from statistical techniques like PCA or Non-Negative Matrix Factorization as shown in the section 2.5. The coefficients are the abundance coefficients associated to the . is the noise vector coming from the sensor or the measurement conditions. Despite the fact Eq. 6 is mathematically equivalent to Eq. 3, it does not hold the same physical sense in the context of this work. Eq. 3 is the model used for representing the interaction between peanut and wheat flour in a simulation context. Eq. 6 is the model used to represent the mixture context of a pixel for a detection purpose.

In the context of a subpixel detection problem, two competing hypotheses are tested. For the null hypothesis associated with the background class, the pixel is assumed to contain only the background sample. Consequently, the LMM decomposes the pixel spectrum according to the sources of variability associated to the background sample . For the alternative hypothesis associated with the target class, the pixel contains background particles as well as target particles. Thus, is modeled using the LMM with the sources of variability associated with the background and the target . The MSD is thus based on the following statistical test:

:

:
where and define the dimensionality associated to the variability for the background and the target respectively. Two matrices are defined corresponding to these hypotheses: contains the vectors in columns, and S contains the vectors and in columns as shown below.

and

The generalized likelihood ratio approach gives the detection statistic for the MSD [START_REF] Manolakis | Detection algorithms for hyperspectral imaging applications[END_REF]:

(Eq. 7)

where and are the projection matrices on the orthogonal subspace of and respectively. These projectors are obtained from the and matrices by the following formula:

(Eq. 8)

After calculating the detection statistic for each pixel spectrum of a sample using (Eq. 7), a threshold must be applied to classify between both classes: target (the pixel contains target and background) and background. This threshold is chosen using the Neyman-Pearson approach that consists of maximizing the detection rate by keeping the false alarm rate under a given limit. For this purpose, the threshold is fixed as the maximum value of the detection statistic obtained on one replicate of the pure wheat sample:

(Eq. 9)

with .

The two other pure wheat replicates were used in order to assess the robustness of the thresholding method.

For both data matrix and , a PCA on non-centered data was performed and the first loadings were extracted to obtained vectors on wheat flour and on peanut flour.

and correspond to the number of extracted loadings for wheat and peanut flours respectively.

As Eq. 7 and 8 show, the design of the MSD highly depends on the design of and . The way of constructing these matrices depends on the choice of the spectral profiles and .

In the method presented in this work, those spectral profiles are obtained by selecting consecutive components from PCA performed on the pure samples. However, the number of components to choose to obtain and : and are tunable parameters that change the calculation of . The design of the MSD consists of finding optimal values for and . In the following, the performance of the detector is qualified using its sensitivity. In this paper, the sensitivity refers to the minimum local concentration (at the pixel scale) for which the detection rate is over 99%.

Software

The data processing is performed using Python 3.7. For data simulation, the PCA is performed using the Scikit-Learn 0.18.1 implementation consisting of a Singular Value Decomposition (SVD). For the design of MSD, the PCA on non-centered data is performed using the eigenvalue decomposition of performed with Numpy 1.16.4. Figure 2 shows that the simulated data exhibit less variability on PC 2 compared to pure peanut and, even less so, to pure wheat. As Eq. 5 shows, the amount of variability added to the model depends on the parameter which represents the simulated target concentration. Consequently, when the simulated concentration is low, the amount of variability is more similar to the one of wheat flour. Moreover, the difference observed on the variability level of pure wheat flour and simulated data is explained by the Gaussian simplification of the real distribution of the data. Indeed, the real distribution of data points in the PCA may not fit a Gaussian perfectly. More particularly, extreme points on PC 2 seem to be more present on real data. These results show that simulated data can be relevantly used for the estimation of the sensibility of detectors.

Results and discussions

Data simulation for the detector design

[Insert Figure 2] Table 1: Design parameters of three detectors of interest. shows the details of the design (the values for and ) for three detectors. These detector designs are selected because they show the best and most interesting results among all those which have been tested. The next section shows the results for other design and focuses on the choice of the parameters.

Table 1: Design parameters of three detectors of interest.

Name of the detector

Detector 1 Detector 2 Detector 3

1 2 2 1 1 2
The design of the detector requires to evaluate its sensitivity to optimize the choice of the parameter values ( and ). Figure 3 shows the detection rate of the three different detector designs described in Table 1: Design parameters of three detectors of interest.. The detection rate indicates the fraction of detected targets for a given peanut concentration of the simulated data. For zero peanut concentration, spectra from real wheat flour images were used. The graph shows that all detectors do not have any false alarm on real wheat measurements. This means that the thresholding method is robust for all three detector designs. Detectors 2 and 3 reach a detection rate of 100% for a simulated peanut concentration of 20% and they both have similar detection rates for smaller peanut concentrations. Detector 1 exhibits a smaller detection rate for every concentration, and it does not reach a 100% detection rate for 20% of peanut adulteration. This shows that, according to the simulated data, the detectors 2 and 3 have a similar sensitivity which is higher than the one of the detector 1.

[Insert Figure 3] reference value at the pixel scale, the real position of the targets is unknown. However, the fact that neighbor pixels are simultaneously detected strengthens the probability that there is effectively peanut in these pixels. In other words, the detection of a neighborhood of pixels is more credible than the detection of an isolated pixel. This argument is developed in the section Detection position. On the other hand, detector 1 only detects one pixel exclusively. The comparison of the detection maps tends to show that detector 1 is less sensitive than detector 2.

The maps in Figure 4 show only a few colored pixels which means that detectors have approximately the same performances. The detector 3 detects 13 more pixels than detector 2 and only one (in the middle of the map) can be considered suspicious since it is isolated.

On the other hand, the detector 2 shows only one exclusive detection (in the middle of the map) which is isolated.

These observations made on real data are relevant to the results obtained with the simulation method (Figure 3): the detectors 2 and 3 have similar performance and are more sensitive than detector 1. This shows that simulated data can be relevantly used for designing and analyzing matched subspace detectors. Additionally, detection map comparison shows that detector 3 seems to be more sensitive than detector 2 despite the fact that they exhibit almost equal performances on the simulated data.

[Insert Figure 4]

Figure 4: Detection map comparison on real data (focus on a sample with 2% of peanutreplicate A). The map above shows the comparison between the detectors 1 and 2, below is the comparison between the detectors 2 and 3. Each pixel is colored according to the output of both detectors (see the legend).

Matched subspace detector

Design of the detector

Figure 5 shows the sensitivity of several detectors calculated on the simulated data.

The first graph on the left shows the effect of varying with . This means the background dimensionality is varying whereas the dimensionality of the target model is fixed.

gives poor performances since no spectra are detected even for a peanut concentration of 20%. The best performances are obtained for . Then, increasing leads to lower detection rates. The graph on the right shows the evolution of the sensitivity when fixing . In this condition, the performances of the detectors are identical for and . Then, choosing a higher value for decreases the detection rate. We also show the design with and (see Figure 3) provides a lower sensitivity than the detector 3. Consequently, the results show that there is an optimal design , for the MSD regarding the performances on the simulated data.

[Insert Figure 5]

Figure 5 : Detector sensitivity for varying and . On the left, is fixed and goes from 1 to 4. On the right, goes from 1 to 4 and is fixed. Detection rate are calculated on simulated data for concentration from 5% to 20% and on real wheat measurement data for 0%.

Figure 6 depicts the geometrical interpretation of the matched subspace detector.

For simplicity, we assume the spectrum vector belongs to a 3-dimensional subspace defined by three virtual wavelength bands , and . The target and the background subspaces are assumed to be 1-dimensional and are represented by vectors and . The common subspace is a 2-dimensional subspace represented by the plane . The MSD compares the residuals of the decomposition of under hypothesis and .

Geometrically, the squared residuals under corresponds to which is the norm of the projection of onto the orthogonal subspace of . Similarly, under the alternative hypothesis , the residuals are calculated as the norm of the projection of onto the orthogonal subspace of . As a consequence, the quantities involved in the matched subspace detector formulation can be translated in terms of vector norms on the graph of the Figure 6: and . Then, the matched subspace detector metric can be rewritten as follow:

(Eq. 10)

The last equality of (Eq. 10) can be deduced from the fact that is the orthogonal projection of on the subspace leading to the relationship: .

[Insert Figure 6] ). The target and background subspaces are represented as 1-dimensional and the common subspace is the red plane defined by the union. For more clarity, the triangle OWA is rectangle in W and the triangle OAP is rectangle in P.

The geometrical interpretation is useful to visualize that the matched subspace detector compares two models with two different dimensionalities. If the pixel spectrum only contains wheat, the modeling of on is sufficient. Consequently, is small because adding vectors from the target subspace should not significantly improve the model fitting. However, if the pixel spectrum contains some peanut, adding vectors from the target variability to the model should significantly improve it. The role of the detector threshold is to define the limit for which the distance between both models is large enough to consider that a peanut particle is present or not.

From this interpretation, the results from different detector designs can be explained. When , the models for and become highly unbalanced. For instance, for the design where , and , is modeled using a 2-dimensional subspace under compared to a 5-dimensional subspace under . This unbalanced situation always holds whatever the design of the MSD. For this reason, the matched subspace detector metric does not compare directly the residuals under both hypotheses but a ratio. The role of the thresholding phase is to take into account this unbalanced problem. One method consists of tuning the threshold on the detector statistic obtained on a pure background sample. This method works when the histograms of the detector statistic between the pure wheat flour and the adulterated pixels can be separated. However, when the design is too highly unbalanced, the model fitting for and becomes very competitive because additional vectors from the peanut subspace can be used to improve the fitting on even if it only contains wheat. This is a consequence of the fact the spectral signatures of peanut and wheat flours are similar. In these conditions, it becomes difficult to find a threshold value that meets the needs for a high detection rate.

When and are too high, each model takes into account a large variability. However, this is not a good strategy since peanut and wheat spectral signatures are similar and their variabilities are high. Thus, increasing the parameters and leads to more ambiguity for the detector and a high-dimensional model will be preferred. Finally, when and are too small ( and ), the model cannot consider the background variability which is detrimental for this application.

Finally, the interpretation of the MSD metric is useful to understand why the set of appropriate values for the parameter and is relatively small. The results on the simulated data are relevant with this interpretation and enable to find the best parameters for the design of the detector.

Number of detections

Figure 8 shows the detection rate of three selected detectors for all the measured samples. The scatter plot shows the detection rate increases when the sample concentration increases so the application of the MSD detector on real measurement is relevant. The results also show a high variance in the detection rates for the same detector applied to the three replicates of the same concentration. This can be explained by the experimental conditions. Indeed, hyperspectral measurements are representative of the material through a depth of some millimeters. However, the peanut concentration is a global characteristic of the sample volume. As it is not possible to make sure that the sample is homogenous, the global concentration and the apparent concentration at the surface may not be equal.

Consequently, sample replicates may exhibit real concentration variance at the surface.

Additionally, the results show a nonlinear behavior in the evolution of the detection rate according to the peanut concentration. This is visible on the left-graph of Figure 8. There are three phenomena to take into account to explain this behavior.

1) Let us assume a sample is perfectly homogeneous at the pixel-scale with a concentration of 20% and a detector with a sensitivity of 10% is applied on a hyperspectral data cube. Then, each pixel has a contribution of 20% of peanut and is detectable. As a consequence, every pixel is detected and the detection rate is 1.

Conversely, if the concentration is 5%, no pixel is detected and the detection rate is 0. Figure 7 shows the explained behavior and the resulting detection rate curve. This phenomenon may explain part of the nonlinearity of the results. This highly depends on the scale of scrutiny for which the sample is declared to be homogeneous: that is to say, the sampling size for which the homogeneity is guaranteed.

[Insert Figure 7]

Figure 7 : A virtual Matched Subspace Detector is considered to have a pixel-wise sensitivity of 10%. The sample is assumed to be perfectly homogeneous so that the situation of every pixel is identical: there is one target particle in the pixel field of view for the sample with 5 % of peanut, and there is four target particles in the pixel field of view for the sample with 20% of peanut.

2) More realistically, the fact a sample has a global concentration of 20% does not mean each pixel surface has the same concentration. In other words, it is reasonable to assume the sample is heterogeneous at the pixel-scale. Assuming there is no spatial relationship between pixels, an image of 100 000 pixels can be statistically considered as 100 000 independent experiments. Each one can be seen as a series of Bernoulli processes for which there are as many trials as the number of particles in the pixel. The probability of selecting a peanut particle corresponds to the global concentration of peanut. The pixel-wise concentration is thus given by a binomial distribution. Such a simulation provides a detection rate curve shown on the left graph of Figure 8. The result is obtained by simulating a sensitivity of 25%. This is relevant to the estimated sensitivity using the simulated data.

3) In practice, pixels are not independent, so several neighbor pixels are likely detected as targets. This phenomenon may happen for two main reasons. Firstly, because flour samples do contain several particle sizes. Some may be higher than 150 µm.

With such a size, some particles may overlap several pixels and make all of them detectable. Secondly, because the particle size study of flours shows that particles tend to agglomerate with each other. Despite the fact the median particle size is approximately 50 µm in wheat and peanut flours, some agglomerations may have a size of several millimeters which leads to the coverage of several pixels. The agglomeration occurs particularly often for peanut flour because of the remaining fatty acids.

[Insert Figure 8]

Figure 8: The detection rate for the three selected detectors for all samples (from 0.02% to 20% of peanut concentration). The graph on the right shows a focus in the concentration range from 0.02% to 2.5%. The dot line is the result of the statistical simulation performed with a virtual detector of sensitivity 25%.

These arguments show all detector designs provide relevant results regarding the real samples with different concentrations. However, the relationship between the number of detections and the peanut concentration of samples is complex. This is related to the scale for which the sample homogeneity is assumed as well as the sensitivity of the detector and the particle size.

Detection positions

Previous results show that detector 3 seems to provide the most sensitive results.

Figure 9 shows the detection maps of three different concentrations and three replicates.

These maps show that the number of detections is repeatable among the replicates as Figure 8 showed. They also show the detection locations are credible: most of them are made on neighbor pixels so that peanut agglomeration can be seen. Furthermore, the location of these agglomerations is randomly distributed across the sample. These results

show the MSD can be used to detect difficult targets as peanut flour in wheat and give their position efficiently.

[Insert Figure 9]

CONCLUSIONS

The purpose of this study was to tackle a difficult detection problem dealing with similar materials with high spectral variability and particle size involving subpixel detections. The development of a Matched Subspace Detector was proposed to overcome these difficulties.

The spectral variability was tackled using subspace modeling through PCA whereas the Linear Mixing Model was used to consider the subpixel detection issue. Moreover, data simulation of different peanut concentrations was proposed to provide an estimation of the performances of the detectors. This technique was used to choose the most appropriate design.

As a result, the data simulation method provides realistic data regarding the measurement variability. Detector designs giving a high detector rate for the low concentrated samples were conserved and applied on real measurements. Results show that the MSD and the data simulation were relevant to overcome the detection issue. Despite the lack of local reference values, the number and the position of the detections show that MSD provides reliable results.

Additional work could be provided for further improvements for this kind of detection situation. Firstly, the data simulation process could be improved by selecting a subset of loadings to simulate the data. This may provide more reliability on the simulation. Indeed, the expected sensitivity obtained on simulated data (20%) does not seems to be reached in practice. Then, even if no spatial a priori hypothesis can be done regarding the particle size, the detection results on real data show that most of the detections are made on neighbor pixels. This is because of particle agglomeration which is a phenomenon that applies on very small flour particles. Such an effect could be taken into account to improve the detection by adding some spatial dependence in the detector algorithm. Finally, the statistical simulation for the number of detected pixels according to the peanut concentration may be improved.

For example, the hypothesis that each pixel is independent of the other could be changed to get a more accurate approach.

Figure 1 :

 1 Figure1: Procedure for spectral simulation. For the first two steps (1-2), the procedure is only described for the wheat flour sample. The same procedure is applied on the peanut flour sample.

Figure 2

 2 Figure2shows the factorial plan of PCA performed on pure wheat and pure peanut

Figure 2 :

 2 Figure 2: Simulated data are projected on the score plot of the PCA performed on real measurements of pure samples. Real measurements are shown with empty square markers. Projected simulated data are shown for low concentrations between 5% and 20% of peanut. Only 400 representative data points among 130 000 are plotted for peanut and wheat flour.

Figure 3 :

 3 Figure 3: The detection rate according to the peanut concentration in simulated data for three detector design. For zero peanut concentration, real wheat flour measurement data are used to calculate the detection rate.

Figure 4

 4 Figure4shows the comparison of a detection map for the three detectors presented

Figure 6 :

 6 Figure 6: Geometrical interpretation of the matched subspace detector. A spectrum is represented by a vector x in a 3dimensional subspace (). The target and background subspaces are represented as 1-dimensional and the common subspace is the red plane defined by the union. For more clarity, the triangle OWA is rectangle in W and the triangle OAP is rectangle in P.

Figure 9 :

 9 Figure 9: Detection maps obtained after applying Detector 3 on real samples for three concentration: 20%, 5% and 0.2%. A black pixel means no detection, a white pixel means target is detected.


Christoph Böhm

Pascal Brault
email: pascal.brault@irisa.fr

Quentin Delamare
email: quentin.delamare@irisa.fr

Paolo Robuffo Giordano

Stephan Weiss

COP: Control & Observability-aware Planning

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

A. Motivation & Related Work

Motion planning algorithms are a crucial component of autonomous task execution for an UAV or robots in general.

1) Trajectory Generation: A common approach is to generate a series of way-points in discrete time intervalsa so-called trajectory. These trajectories are often 4D flat outputs (3D position and yaw orientation) of the UAV with their subsequent derivatives. These generation methods often let the UAV move from point A to point B with additional goals and constraints, e.g., reducing flight time or energy consumption. Examples of such generation methods based on continuous-space models can be found in [START_REF] Berscheid | Jerk-limited Real-time Trajectory Generation with Arbitrary Target States[END_REF]- [START_REF] Morbidi | Minimum-Energy Path Generation for a Quadrotor UAV[END_REF], acting on different levels of abstraction and even taking probabilities into account. [START_REF] Chen | Dynamic Obstacle Avoidance for UAVs Using a Fast Trajectory Planning Approach[END_REF]- [START_REF] Candido | Minimum Uncertainty Robot Path Planning using a POMDP Approach[END_REF] create a discrete map of the taskspace, a so-called sampling-based method, to apply search algorithms to find a path from A to B. With the increase of computational power, learning-based approaches attracted attention in recent years. [START_REF] Liu | Control-oriented UAV highly feasible trajectory planning: A deep learning method[END_REF]- [START_REF] Xi | Unmanned Aerial Vehicle Trajectory Planning via Staged Reinforcement Learning[END_REF] use machine learning for trajectory generation with reduced calculation times.

All previously mentioned approaches do not consider model and state uncertainties to ensure robust trajectory tracking or accurate state/parameter estimation. 2) Estimation-aware Trajectories: State estimation with proper system modeling [START_REF] Böhm | Filter-Based Online System-Parameter Estimation for Multicopter UAVs[END_REF] relies on the design of sufficiently informative system input for accurate and fast estimate convergence. The works in [START_REF] Hausman | Observability-Aware Trajectory Optimization for Self-Calibration With Application to UAVs[END_REF]- [START_REF] Ansari | Minimum Sensitivity Control for Planning with Parametric and Hybrid Uncertainty[END_REF] show that taking the estimation or parametric uncertainty into account drastically improves system parameter estimation results while allowing task execution. In this work we will use the Expanded Empirical Local Observability Gramian (E 2 LOG) [START_REF] Preiss | Simultaneous self-calibration and navigation using trajectory optimization[END_REF]. These informative trajectories might show poor robustness against uncertainties in the robot model for the tracking controller that executes them [START_REF] Lee | Geometric Tracking Control of a Quadrotor UAV on SE(3)[END_REF], [START_REF] Kamel | Model Predictive Control for Trajectory Tracking of Unmanned Aerial Vehicles Using Robot Operating System[END_REF].

3) Control-aware Trajectories: To address this issue, the notion of Closed-Loop State and Input Sensitivities (S/I-S) has been recently introduced in [START_REF] Brault | Robust Trajectory Planning with Parametric Uncertainties[END_REF], [START_REF] Giordano | Trajectory Generation for Minimum Closed-Loop State Sensitivity[END_REF] as a suitable metric to be optimized. Minimizing the norm of the S/I-S generates a trajectory whose tracking results are minimally sensitive to model uncertainties of the robot states and inputs. This is important as it increases the robustness and accuracy of the trajectory tracking and improves the repeatability of the control inputs when the model parameters vary. The S/I-S, however, needs knowledge of the actual robot state and nominal values of the model parameters, which may not be directly available but can be provided by state estimation.

4) Link & Combination: Therefore, a link between observability-related metrics and S/I-S metrics exists, since the evaluation of S/I-S needs a good knowledge of states and parameters, and the tracking of maximally observable trajectories benefits from increased robustness in the trajectory execution. That fact motivates our study on how to combine both methods in a unified trajectory optimization problem taking into account that these two objectives can conflict among themselves, as seen in Sec. IV.

B. Contributions

This paper proposes Control & Observability-aware Planning (COP) as a way to balance two possibly contradicting optimization problems, namely (i) generating trajectories whose execution is minimally sensitive to model uncertainties and (ii) generating trajectories that can be sufficiently informative for accurate state estimation. We present a method to address these problems by leveraging previous contributions to the topics of observability-aware and minimally-sensitive trajectory planning, combining them in the formulation and resolution of a Single-Objective Optimization Problem (SOOP) based on the Augmented Weighted Tchebycheff method. In a case study, considering the state estimation and robust trajectory tracking for a 3D quadrotor UAV, we discuss the statistical results of a realistic simulation campaign that shows the potential of the proposed contribution.

II. PRELIMINARIES

A. Quadrotor model

Let us consider a frame M as the quadrotor body frame, attached to its center of mass (CoM), with its z-axis z M aligned with the thrust of the four rotors. The state vector x of this system consists of its linear position r = (x, y, z) and velocity v = (v x , v y , v z ), both expressed in the world frame W. It also includes the body orientation expressed through the unit length quaternion q = (q w , q x , q y , q z ) (Tait-Bryan angle definition with yaw first (312-sequence)) as well as its angular velocity ω = (ω x , ω y , ω z ), expressed in the body frame M, therefore

x = [r T , v T , q T , ω T ] T ∈ R 13 .
As the quadrotor's orientation is expressed by the Hamiltonian quaternion q, we recall that z M = R(q)z W = q -1 ⊗ z W ⊗ q, where ⊗ is the quaternion product. R(q) is the rotation matrix as a function of the quaternion q.

It is possible to link the squared rotor speeds of the quadrotor (control inputs

) u = [ω 2 1 ω 2 2 ω 2 3 ω 2 4 ]
T to [f, τ T ] T , the total effective thrust and torque. These are related by the allocation matrix S, [f, τ T ] T = Su (e.g., [START_REF] Mahony | Multirotor Aerial Vehicles: Modeling, Estimation, and Control of Quadrotor[END_REF] Eq. ( 8)), which includes the rotor thrust force coefficient k f , the drag moment coefficient k m , and the arm length from the center of mass to each motor/rotor group.

With these definitions, compared to [START_REF] Brault | Robust Trajectory Planning with Parametric Uncertainties[END_REF], the quadrotor dynamical model is

ẋ =          ṙW = v W vW = -gz W + f m z M q = 1 2 q ⊗ ω M ωM = J -1 (τ -ω M × (J ω M )) (1)
where m and J are the quadrotor's mass and its inertia matrix, respectively. g is the Earth's gravitational pull. The allocation matrix S implies that the dynamics are also effected by the parameters k f , k m and . The quantities k f and k m are aerodynamic parameters that depend on the rotors characteristics and the aerodynamic interaction with the environment (e.g., presence of wind, ground effects). Therefore, an accurate value for these parameters can be hard to obtain, and we thus consider p = [k f k m ] T ∈ R 2 as the uncertain parameters of the dynamical model.

B. Tracking controller

With the dynamic model detailed in the previous subsection, we can now present the controller that has been used in this work. The chosen control task is to let the system output y(x) = [x y z ϕ]

T ∈ R 4 track a desired motion r d (a, t) ∈ R 4 , where ϕ is the yaw (or heading) angle of the quadrotor. The DFL (Dynamic Feedback Linearization) controller with an integral term used in previous works [START_REF] Brault | Robust Trajectory Planning with Parametric Uncertainties[END_REF], [START_REF] Giordano | Trajectory Generation for Minimum Closed-Loop State Sensitivity[END_REF] is not robust against parameter uncertainties and is not capable of considering input constraints. Therefore, we chose to use a different controller, namely the so-called Lee controller [START_REF] Lee | Geometric Tracking Control of a Quadrotor UAV on SE(3)[END_REF], which performs slightly worse in an ideal case compared to the DFL, but it is much less complex to implement and tune.

Although we implemented the same structure of the controller as in [START_REF] Lee | Geometric Tracking Control of a Quadrotor UAV on SE(3)[END_REF], we added some minor changes in order to match the dynamics of the quadrotor in 3D (especially the use of the quaternion q). In particular we consider attitude and angular velocity errors defined as

e q = 1 2 R T d R(q) -R(q) T R d ∨ and e ω = w. ( 2
)
The resulting control inputs are then

f = (-k r e r -k v e v -k i ξ + m(gz w + rd ))•R(q)z w , (3)
τ = -k q e q -k ω e ω , (4)
where e r and e v are the position and velocity errors, and ξ = [ξ x ξ y ξ z ] T is the position integrator, and k r , k v , k i , k q , k ω are suitable control gains. We then compute u via the inverse of the allocation matrix, u = S -1 [f, τ T ] T .

C. Curve representation

The controller is designed to let the quadrotor follow a reference trajectory r d (a, t), where a is the parameter vector for the chosen class of curve. In [START_REF] Giordano | Trajectory Generation for Minimum Closed-Loop State Sensitivity[END_REF] 'plain polynomials' were used, with the drawback of introducing possible numerical instability during the optimization. Due to this reason, we switched in [START_REF] Brault | Robust Trajectory Planning with Parametric Uncertainties[END_REF] to the use of Bézier curve representation, as they are more stable from a numerical point of view.

This work goes one step further by implementing piecewise Beziér curves for the trajectory representation to avoid the use of a single high degree Beziér curve as in [START_REF] Brault | Robust Trajectory Planning with Parametric Uncertainties[END_REF]. An aditional abstraction happens on the parameter vector a, which now contains way-points (with velocity, acceleration, and even higher order constraints) instead of the control points themselves. These way-points define the curve at the start, end, and in between curve pieces.

Let B i∈ 1, n , with n 1 number of pieces, be the Bézier curve of degree d 2 shaping the trajectory (C d-1 continuity). In total, there are n + 1 way-points (t i , P i ), where t i is the time associated to the point P i ∈ R n dim ×n jc , with n dim the number of dimensions of the trajectory (e.g., x, y, z, and yaw angle), and n jc the number of joining conditions (e.g., n jc = 1 for position only, which means every Bézier curve piece is a straight line segment). The degree d of the Bézier curve, the number of joining conditions n jc and the number of control points n c are linked by d = 2n jc -1 = n c -1. With these conditions and the way-points as constraints, one can formulate a linear system of equations that solves for the control points of each Bezièr curve piece.

III. CONTROL & OBSERVABILITY-AWARE PLANNING

A. Objectives for Trajectory Optimization

A trajectory r d (a, t) parameterized by the coefficient vector a, Sec. II-C, can be optimized for different goals by changing its shape. We represent this goal by the so-called utility function U (a) which, in our case, is a scalar cost to be minimized subject to constraints minimize a U (a),

subject to a ∈ A, equ. & inequ. constraints, (5)
with A as the feasible set for the parameter vector.

1) State Sensitivity Metric: The state sensitivity metric as minimization objective was introduced in [START_REF] Giordano | Trajectory Generation for Minimum Closed-Loop State Sensitivity[END_REF] and is based on a generic robot model ẋ = f (x, u, p), where x ∈ R nx is the state vector, u ∈ R nu is the control input vector, and p ∈ R np the vector of system parameters (which are assumed uncertain). This is combined with a tracking control law ξ = g(ξ, x, r d (a, t), p c ) and u = c(ξ, x, r d (a, t), p c ),

where ξ ∈ R n ξ are the internal controller states (e.g., an integrator), p c a nominal value for the parameters p, and r d (a, t) a desired trajectory. The state sensitivity matrix for the closed-loop system (i.e., considering both the robot and chosen control) is defined as

Π(t) = ∂x(t) ∂p p=pc ∈ R nx×np (6)
representing the variation of the states x in relation to variations in the parameter vector p, evaluated on the nominal value p c . We refer the reader to [START_REF] Giordano | Trajectory Generation for Minimum Closed-Loop State Sensitivity[END_REF] for further details. The integral of the matrix norm of the state sensitivity over the whole trajectory (duration T ) can be used to reduce the influence of parameter uncertainty on the states.

U (a) = F Π (a) = T 0 ||Π(t)|| dt (7)
2) Input Sensitivity Metric: As natural evolution, [START_REF] Brault | Robust Trajectory Planning with Parametric Uncertainties[END_REF] added to the state sensitivity metric the so-called input sensitivity metric which is defined as

Θ(t) = ∂u(t) ∂p p=pc ∈ R nu×np (8)
and, similarly to Eq. ( 6), maps how variations of the parameter vector p result in variations of the control inputs. The integral of its matrix norm, again, gives us the objective function Eq. ( 9) which can be used to reshape the trajectory to be less sensitive in its control inputs against parameter uncertainties. Details on the derivations can be found in [START_REF] Brault | Robust Trajectory Planning with Parametric Uncertainties[END_REF].

U (a) = F Θ (a) = T 0 ||Θ(t)|| dt (9)
Note that it is not possible to compute Π(t) in closedform. However, it is possible to obtain a closed-form expression of its dynamics, and other quantities, e.g., Θ(t), can be derived from it. The evolution of Π(t) (consequently Θ(t) too) over a time interval T of interest then can be obtained by numerically integrating over time [START_REF] Brault | Robust Trajectory Planning with Parametric Uncertainties[END_REF], [START_REF] Giordano | Trajectory Generation for Minimum Closed-Loop State Sensitivity[END_REF]. One of the main hypotheses/key aspects in [START_REF] Brault | Robust Trajectory Planning with Parametric Uncertainties[END_REF], [START_REF] Giordano | Trajectory Generation for Minimum Closed-Loop State Sensitivity[END_REF] is to consider that the whole state of the system is known during the tracking of the desired trajectory.

3) Observability Metric: Online state estimation is one way of making the state x and parameters p available at run-time [START_REF] Böhm | Filter-Based Online System-Parameter Estimation for Multicopter UAVs[END_REF], [START_REF] Weiss | Versatile Distributed Pose Estimation and Sensor Self-Calibration for an Autonomous MAV[END_REF], [START_REF] Hausman | Self-Calibrating Multi-Sensor Fusion with Probabilistic Measurement Validation for Seamless Sensor Switching on a UAV[END_REF].

How well such estimates perform depends on the accuracy of the system model and sensors used, but control inputs given to the robot are equally important. The E 2 LOG [START_REF] Hausman | Observability-Aware Trajectory Optimization for Self-Calibration With Application to UAVs[END_REF], [START_REF] Preiss | Simultaneous self-calibration and navigation using trajectory optimization[END_REF] works on the idea of the quality of observability, proposed in [START_REF] Hermann | Nonlinear controllability and observability[END_REF], [START_REF] Krener | Measures of Unobservability[END_REF], which evaluates it over a whole trajectory r d (a, t) with a duration T . It uses a n-th order Taylor expansion to approximate the Jacobian matrix which models the sensitivity of the measurements with respect to the control inputs, the state and its changes on a small time horizon H ( W t 0 ,H (a) in [START_REF] Preiss | Simultaneous self-calibration and navigation using trajectory optimization[END_REF]). Summing all these quality measure segments along the trajectory r d (a, t) gives us

W O (a) = N k=0 W k∆t,∆t (a) ∈ R nx×nx (10)
with ∆t = T N and N the number of trajectory segments. The objective is to improve the least sensitive state (or combination of states) through the smallest eigenvalue of W O (a). Adding a minus to the smallest eigenvalue makes it usable in a minimization problem as objective function

U (a) = F E 2 LOG (a) = -σ min ( W O (a)). ( 11
)
The result is a trajectory with optimum observability properties, improving the states' convergence by lowering the uncertainty and increasing the overall accuracy.

The estimator in [START_REF] Böhm | Filter-Based Online System-Parameter Estimation for Multicopter UAVs[END_REF] is the base for the derivation of the observability-aware trajectory optimization in our quadrotor UAV case study, Sec. IV.

B. Multi-Objective Optimization Problem

The problem presented in this work is an example of a Multi-Objective Optimization Problem (MOOP) trying to optimize for different objectives with constraints.

1) Pareto Optimality: Ideally, one would try to find the non-dominated set in the entire feasible set A for the parameter vector a, a so-called Pareto optimal set or Pareto front [START_REF] Steuer | An interactive weighted Tchebycheff procedure for multiple objective programming[END_REF]- [START_REF] Holzmann | Solving discrete multi-objective optimization problems using modified augmented weighted Tchebychev scalarizations[END_REF]. The term non-dominated or optimal means that the current set does not improve one objective while worsening another. We compute only one point in the Pareto optimal set due to the complexity of the objective functions and the resulting computation times, see Sec. IV-B. [START_REF] Morbidi | Minimum-Energy Path Generation for a Quadrotor UAV[END_REF].

In [START_REF] Brault | Robust Trajectory Planning with Parametric Uncertainties[END_REF], Linear Scalarization Problem (LSP) was chosen as a method because it achieved good results in balancing the S/I-S of a 2D quadrotor. LSP combines objective functions to a single cost through a linear weighting of each objective. This is only a feasible option if the Pareto front is convex. To be more specific, in the case of concave Pareto fronts, LSP tends to converge towards extrema solutions, an optimal solution for only one of the objectives. Evenly distributed weights do not produce an evenly distributed representation of the Pareto optimal set. The addition of the observabilityawareness through the E 2 LOG as a third objective is the natural evolution of the approach, however, non-convexity (concavity) is possible with this addition.

2) Augmented Weighted Tchebycheff Method: To solve this issue, we balance the objectives through their distance between the objective value F i (a) and an aspiration point F O i , compared to [START_REF] Brault | Robust Trajectory Planning with Parametric Uncertainties[END_REF]. This aspiration point

F O i
is the individual objective's minimal value from Eq. ( 5) -called utopia point. Another important point in the Pareto set is the so-called nadir point

F N i = max 1≤j≤k {F i (a O j )}.
It is the largest cost of all objectives with respect to the j-th utopia point.

U (a) = max i {λ i |F i (a) -F O i |} + ρ k j=1 |F j (a) -F O j | (12)
The utility function in Eq. ( 12) includes F(a) = {F Π (a), F Θ (a), F E 2 LOG (a)}, and the scale of each objective

λ i = w i |F N i -F O i | .
Every point on a Pareto front is a minimum of the Tchebycheff function for some λ i (convex or nonconvex) and achieves Pareto optimality of the solution. The weight w i is the user defined preference of one objective, which selects one solution from the possible set of Pareto optimal solutions (bias),

k i=1 w i = 1. |F N i -F O i | in the denominator of λ i normalizes the Tchebycheff function to the interval [0, 1].
According to [START_REF] Steuer | An interactive weighted Tchebycheff procedure for multiple objective programming[END_REF], ρ values should be selected between 0.0001 and 0.01. We refer the reader to [START_REF] Wierzbicki | Multiple Criteria Games -Theory and Applications[END_REF], [START_REF] Marler | Survey of multi-objective optimization methods for engineering[END_REF], [START_REF] Emmerich | A tutorial on multiobjective optimization: fundamentals and evolutionary methods[END_REF] for further reading. As can be seen, this method reduces the MOOP into a SOOP. We will discuss how to use Eq. [START_REF] Böhm | Filter-Based Online System-Parameter Estimation for Multicopter UAVs[END_REF] in Sec. III-C.3 and Sec. III-C.4.

C. Implementation

COP uses a multi-step approach to the trajectory optimization, see Fig. 1, as the utopia point F O i and nadir point F N i of each objective F i need to be known before combining objectives. The framework is implemented in Python and uses a local derivative-free optimization, namely Constrained Optimization BY Linear Approximations (COBYLA) of the open-source library for nonlinear optimization (NLOPT). The numerical integration method dopri5 of SciPy enables us to calculate the individual costs (Eq. ( 7), Eq. ( 9), and Eq. ( 11)) during each iteration of the optimization. 1) Preconditioning: This important step ensures that the initial trajectories are dynamically feasible and let the system reach the final target accurately. The preconditioning starts with a random trajectory with initial way-point P 0 (for a 3D quadrotor UAV position and yaw orientation), target way-point P T , duration T , and n number of Bézier curve pieces supplied. All way-points between the initial and target represent the decision variables a i of the optimization, and can be chosen freely within the admissible set A. Then, dynamical constraints are applied to the trajectory through a short optimization (e.g., min. and max. rotor speeds). Afterwards, we precondition the trajectory to take the controller tracking imperfections into account, ensuring that the system reaches the final target accurately with the nominal value p = p c , thus, without uncertainties in the model. The vector a INIT parameterizes the shape of the initial trajectory.

2) Optimizing Individual Objectives: As mentioned before, to allow the combination of multiple objectives, one needs to compute F O i and F N i of each objective first. We do the minimization of each objective F Π , F Θ , and F E 2 LOG , defined in Eq. ( 7), Eq. ( 9), and Eq. ( 11) respectively, along with Eq. ( 5) to get those needed extrema points. The results are the utopia point and coefficient vector of each objective,

F O Π /a O Π , F O Θ /a O Θ , and F O E 2 LOG /a O E 2 LOG
, respectively. These coefficient vectors are used to calculate the nadir points of the objectives F N Π , F N Θ , and F N E 2 LOG . 3) State and Input Sensitivities Optimization: A first example for Eq. ( 12) is the computation of the S/I-S. The balancing of state and input sensitivities similar to [START_REF] Brault | Robust Trajectory Planning with Parametric Uncertainties[END_REF] is possible by setting w S/I-S = 1 2 , 1 2 , 0 and using the utopia and nadir points of F Π and F Θ for normalization. The rationale behind w S/I-S is as follows: (i) Eq. ( 7) aims at minimizing the state sensitivity over the trajectory, so that the tracking accuracy of r d (a, t) is made most insensitive to uncertainties in the model parameters; (ii) Eq. ( 9) aims to minimize the input sensitivity during the whole trajectory, in order to obtain control inputs that are most insensitive against variations of the robot parameters. For this optimization we chose ρ = 0.0001 in Eq. ( 12) and the results of this step are F O S/I-S /a O S/I-S . Both these objectives provide the ability to generate control-aware trajectories.

4) Control & Observability-aware Optimization: The possible antagonistic nature of the S/I-S and E 2 LOG needs utility functions like Eq. ( 12) to be balanced successfully. We weight all three objectives equally, w COP = 1 3 , 1 3 , 1 3 , in Eq. ( 12) with ρ = 0.0001. This makes the combination of State and Input Sensitivity-based with Observability-aware trajectory planning possible and results in F O COP /a O COP . The even distribution of weights might not result in equally improved objectives in practice, due to the possible skewed normalization from the approximation of nadir and utopia points. To ensure we get proper trajectories that reduce both, we apply filtering based on the costs history available from each optimization run -a posterior preference

F S/I-S (a O COP ) < F S/I-S (a INIT ) and F E 2 LOG (a O COP ) < F E 2 LOG (a INIT ).

IV. RESULTS

To show the use of our proposed approach, we conduct a case study focusing on a 3D quadrotor UAV's rotors thrust force and drag moment coefficients

k f and k m (p = [k f k m ]
T ). These coefficients have a significant impact on the tracking performance of the control and are hard to estimate.

A. Setup & Method

The system model of Sec. II-A, control law of Sec. II-B, and optimization of Sec. III are implemented in Python. All system parameters are based on the quadrotor Hummingbird model of [START_REF] Böhm | Filter-Based Online System-Parameter Estimation for Multicopter UAVs[END_REF]. This previous work shows that it is possible to estimate the thrust force and drag moment coefficients, k f and k m . The method of using simulations allows for better repeatability of experiments and avoids the introduction of other artifacts due to uncertainties of the real system.

In this empirical evaluation we look at four types of trajectories: (i) the initial preconditioned trajectory (INIT); (ii) the S/I-S objective optimized one; (iii) the E 2 LOG objective optimized trajectory; (iv) the new COP objective optimized trajectory. The initial one serves as baseline for all other trajectories. Each trajectory has a duration of T = 20s, 5 Bézier curve pieces, and a random target way-point P T in 3D space with (U(2, 5), U(2, 5), U(-0.5, 1)) in meters.

In total, optimizations have been completed for 20 targets, each for the S/I-S, E 2 LOG, and COP objective. The individual cost functions evolutions, positional mean integral error norms, and estimation uncertainties are evaluated from this set of trajectories. For the positional mean integral error norm evaluation, we chose to randomly perturb the coefficients k f and k m in the ranges of ±1 % and ±5 % and fulfill 30 closed-loop flights, changing them for every trajectory. Greater perturbations are not considered as such a deviation from the nominal value might hint at problems at the parameter identification/estimation. The estimation uncertainty results are based on 10 runs of each trajectory with different randomly wrong initial guesses (±30 %) from ground truth (k f = 3.375 × 10 -4 N/s -2 and k m = 0.016 m). As E 2 LOG and S/I-S have different orders of magnitudes, a normalization with the initial cost value was performed. In addition, we chose to use the inverted values of E 2 LOG as they can grow unbound and allows for better comparisons. A decrease (<1) means an improvement, while an increase (>1) indicates a decline in the performance of the respective objective. As the optimization uses a gradient-free method, the average shows some spikes.

Fig. 2 presents the results of the 20 individual COP runs. The top left plot shows the average and 1σ standard deviation of the S/I-S-based optimization, Sec. III-C.3, with its cost in blue and the E 2 LOG's cost in orange. One can see the decrease in the sensitivity cost, meaning that the state and input sensitivities are minimized (as expected). This comes, however, at the expense of the overall quality of observability, indicated by the increase of the inverse E 2 LOG cost. Therefore, these results seem to indicate that S/I-S and E 2 LOG can be two conflicting objectives. On the top right are both costs, again blue S/I-S and orange E 2 LOG, depicted for the quality of observability optimization, Sec. III-C.2. As we chose the inverse here, a decrease is equal to an improvement of the E 2 LOG. Once again, we see the behavior of the left plot reflected in this optimization as well. From those two plots one can infer that if one improves the other might get worse. The two plots at the bottom of Fig. 2 are the results of the optimization runs based on the COP objective, Sec. III-C. [START_REF] Morbidi | Minimum-Energy Path Generation for a Quadrotor UAV[END_REF]. We can see that even distributed weights, as described in Sec. III-C.4, can still result in slightly skewed solutions caused by the approximation of the individual utopia and nadir points F O i and F N i . The solutions themselves are Pareto optimal, and the straightforward filtering ensures an overall decrease of all considered objectives. This indicates that COP can balance and decrease all objectives.

2) Control Tracking Error: The evaluation of the tracking performance of the system with its controller is based on the aforementioned set of trajectories, and is done by considering the positional mean integral error norm of position r(t) with respect to the desired position r d (t).

The results can be seen in Fig. 3 where we depict the quartile box plots from the statistical data. As mentioned before, each trajectory is flown in simulation 30 times with a changed set of k f and k m for each flight, normally distributed around ±1 % and ±5 % of their nominal values (which is used to evaluate the various sensitivity quantities).

The quartile boxplots of Fig. 3 show the average tracking performance of each trajectory with different amplitudes of perturbation represented by the mean of the integral of the error norm at each point on the trajectory. The median of the S/I-S optimized trajectory performs the best and the E 2 LOG one the worst, with COP performing between those two, which confirms Fig. 2 and our expectations. COP can not reach the same performance level as S/I-S because of the balancing in Eq. ( 12), however, we can improve estimation performance at the same time, Fig. 4. Looking at the evolution of those boxplots, one can see that the farther away p gets from p c the less difference is between each objective. This is also expected as the optimization evaluates the trajectories at nominal value. Therefor, one might start with estimation-aware trajectories to get p as close to p c and then switches to control-aware ones for improved tracking.

3) Estimation Error: We used IEKF implementation in Matlab of [START_REF] Böhm | Filter-Based Online System-Parameter Estimation for Multicopter UAVs[END_REF] to evaluate the influence of the trajectories on the estimation of k f and k m , which has proven to be good at estimating those parameters. The trajectory optimization is used to generate artificial position sensor and IMU measurements together with rotor speed input for each trajectory. Each of these recordings has been tested using the IEKF with initial guesses of nominal values k f and k m perturbed randomly by ±30 % 10 times. Note that both parameters are poorly observable. To visualize the influence, we once again use quartile boxplots for each individual optimization objective. Fig. 4 shows on the left plot the quartile boxplot of k f and on the right for k m . The estimation performance is evaluated by the reduction of uncertainty (represented by the standard deviation) at the end of the trajectory at T . One can see that the S/I-S optimized trajectories perform slightly better than the initial ones. This is because we already gain more motion and excitation from the optimization objective. As expected, based on the optimization data, the E 2 LOG optimized trajectories perform best with the new COP optimized one in between. Note that this is expected, due to the balancing of two objectives we are not able to perform as well as a single objective optimized trajectory. These comparable results were already indicated in Fig. 2. All the presented results support our claim that our proposed COP objective can balance two possibly opposing objectives and maintain good performance. 4) Computation Times: All trajectory optimizations were done on a PC with an AMD Ryzen 5 3600 CPU (6 cores/12 threads), 16GB Ram, and NVidia RTX 2070 (Super) GPU. Note that one instance of the COP implementation only uses one CPU thread, and therefore, parallelization is possible by starting multiple instances of COP. The duration of each optimization run and each objective were logged. The initial trajectory generation is done in under 120 s. Minimizing the state and input sensitivity metrics (F Π (a) and F Θ (a)) takes on average 38 min and 15 min, respectively. Improving the observability through the E 2 LOG (F E 2 LOG (a)) needs on average 26 min. The trajectory optimization with S/I-S (F S/I-S (a)) finishes in 36 min on average. The last objective of the COP approach, minimizing both E 2 LOG and S/I-S, runs for 27 min on average. One complete optimization run needs in total around 2.4 hours which is due to the numerical integration of the complex metrics over the whole trajectory during each iteration of the optimization.

V. CONCLUSION

This paper wanted to answer the question: How to combine Closed-Loop State and Input Sensitivity-based with Observability-aware trajectory planning? Our proposed Control & Observability-aware Planning (COP) framework and its statistical evaluation provide an answer to it.

Intuitively, taking state and input sensitivities into account during the trajectory generation might result in noninformative trajectories for state/parameter estimation. However, informative motions for the estimation process are often difficult to control and likely cause higher tracking errors. Our statistical case study of a 3D quadrotor UAV, focused on system parameters that have a significant impact on the tracking error and are poorly observable (thrust force coefficient k f and drag moment coefficient k m ), provided insights into the negative correlation between closed-loop state and input sensitivity-based and observability-aware trajectory planning. We were able to show that both objectives work against each other, meaning that while one can improve, the other will get worse. Applying the Augmented Weighted Tchebycheff Method in a multi-step approach to such Multi-objective Optimization Problem (MOOP) allows to balance both in a Single-Objective Optimization Problem (SOOP), improving trajectory tracking and, at the same time, state estimation.

To summarize, we have successfully shown that it is important to consider both objectives as they correlate with each other. The insights and results in this paper give a motivation to go forward with more sophisticated MOOP approaches. Further research should be conducted towards a real-world closed-loop flight using the estimation in the optimization as feedback.

Fig. 1 .

 1 Fig. 1. Graphic overview of the multi-step Single-Objective Optimization Problem (SOOP) solved with Control & Observability-aware Planning (COP), a i being the resulting trajectory coefficients after each operation.

Fig. 2 .

 2 Fig.2. Overview of the cost function evolution during the optimization depicted by averages and 1σ standard deviations over 20 optimization runs with different initial trajectories. All cost values are normalized with the initial cost prior to the mean calculation. S/I-S cost (blue), E 2 LOG cost (orange), COP cost (green). All optimizations minimize their respective cost function, therefore, a decrease below 1 (dotted line) is an improvement of the respective objective. Note that the inverse of the E 2 LOG is used.

Fig. 3 .

 3 Fig.3. Quartile box plots showing the positional mean integral error norm over the whole trajectory, average over 20 trajectories each with 30 simulated closed-loop flights. The two plots show the influence of different perturbation amplitudes on the parameters k f and km. As expected, the S/I-S optimized one performs best followed by the COP. The E 2 LOG based trajectories perform worse as they just improve estimation performance. COP as well as S/I-S are most effective with small perturbations, because the sensitivity is evaluated at p = pc.B. Discussion1) Cost Function Behavior: We recorded each objective function's cost value (F INIT , F S/I-S , F E 2 LOG , F COP ) at all iteration steps during each optimization run to evaluate the behavior of the costs by averaging all runs.As E 2 LOG and S/I-S have different orders of magnitudes, a normalization with the initial cost value was performed. In addition, we chose to use the inverted values of E 2 LOG as they can grow unbound and allows for better comparisons. A decrease (<1) means an improvement, while an increase (>1) indicates a decline in the performance of the respective objective. As the optimization uses a gradient-free method, the average shows some spikes.Fig.2presents the results of the 20 individual COP runs. The top left plot shows the average and 1σ standard deviation of the S/I-S-based optimization, Sec. III-C.3, with its cost in blue and the E 2 LOG's cost in orange. One can see the decrease in the sensitivity cost, meaning that the state and input sensitivities are minimized (as expected). This comes, however, at the expense of the overall quality of observability, indicated by the increase of the inverse E 2 LOG cost. Therefore, these results seem to indicate that S/I-S and E 2 LOG can be two conflicting objectives. On the top right are both costs, again blue S/I-S and orange E 2 LOG, depicted for the quality of observability optimization, Sec. III-C.2. As we chose the inverse here, a decrease is equal to an improvement of the E 2 LOG. Once again, we see the behavior of the left plot reflected in this optimization as well. From those two plots one can infer that if one improves the other might get worse. The two plots at the bottom of Fig.2are the results of the optimization runs based on the COP objective, Sec. III-C.[START_REF] Morbidi | Minimum-Energy Path Generation for a Quadrotor UAV[END_REF]. We can see that even distributed weights, as described in Sec. III-C.4, can still result in slightly skewed solutions caused by the approximation of the individual utopia and nadir points F O i and F N i . The solutions themselves are Pareto optimal, and the straightforward filtering ensures an overall decrease of all considered objectives. This indicates that COP can balance and decrease all objectives.2) Control Tracking Error: The evaluation of the tracking performance of the system with its controller is based on the aforementioned set of trajectories, and is done by considering the positional mean integral error norm of position r(t) with respect to the desired position r d (t).

×10 - 3 Fig. 4 .

 34 Fig. 4. Quartile box plots showing the IEKF's uncertainty based on the state's standard deviation at the end of the trajectory, average over 20 trajectories (initial, S/I-S optimized, E 2 LOG optimized, and E 2 LOG optimized each) closed-loop flights. Each trajectory is tested with 10 different initial guesses ±30 % uniformly distributed around ground truth. (left) is the thrust force coefficient k f and (right) the drag moment coefficient km depicted. As expected, the initial trajectory has the worst estimation performance and the E 2 LOG ones are the best with COP generated motions being comparable.

C. Böhm and S. Weiss are with the Control of Networked Systems Group at the University of Klagenfurt, Austria. email:

This work was partially supported by the project ANR-20-CE33-0003 "CAMP" and by the Austrian Ministry for Transport, Innovation and Technology (BMVIT) under the grant agreement 878661 (SCAMPI)