Andres Aravena
email: andres.aravena@dim.uchile.cl

Carito Guziolowski
email: carito.guziolowski@bioquant.uni-heidelberg.de

Anne Siegel
email: anne.siegel@irisa.fr

Alejandro Maass
email: amaass@dim.uchile.cl

Andrés Aravena

Using Mutual Information and Answer Set Programming to refine PWM based transcription regulation network

Keywords: Transcriptional regulatory network reconstruction, mutual information, logical programming

come   L'archive ouverte pluridisciplinaire

Introduction

Regulation plays a major role in the transition from cell genotypes to phenotypes. The transcription of each gene is controlled by a regulatory region upstream of the transcription start site, where transcription factors can bind [START_REF] Alberts | Molecular Biology of the Cell[END_REF]. The chain of interactions between genes, transcription factors (TFs) and binding sites (BS) is known as the transcription regulation network, usually represented as a signed oriented graph [START_REF] Davidson | Gene regulatory networks[END_REF]. Knowledge of this graph is the foundation of the qualitative or quantitative modeling of cell behavior, [START_REF] Barabási | Network biology: understanding the cell's functional organization[END_REF] thus constituting one of the main tools in system biology modeling of development as well as organism maintenance.

Experimental methods are used to identify the regions in the chromosome to which a protein binds, and for purifying and identifying a DNA-bound protein [START_REF] Bouveret | Bacterial interactomes: from interactions to networks[END_REF]. These experiments have not been widely applied because they are expensive and time consuming [START_REF] Baumbach | Reliable transfer of transcriptional gene regulatory networks between taxonomically related organisms[END_REF]. In addition, most of the results of these methods are focused on human and other eukaryotic organisms. However, some results for bacteria are organized in databases such as Prodoric [START_REF] Grote | PRODORIC (release 2009): a database and tool platform for the analysis of gene regulation in prokaryotes[END_REF], RegulonDB, RegTransBase and CoryneRegNet. Besides the direct wet-lab experimentation, there are two main strategies for regulatory network reconstruction: [START_REF] Alberts | Molecular Biology of the Cell[END_REF] directly recognizing regulatory elements in the genomic sequence, and (2) using expression data to determine gene regulatory relationships based on statistical correlation.

The first strategy for in silico network reconstruction uses the experimental sequences to extrapolate to other taxonomically related organisms [START_REF] Baumbach | Reliable transfer of transcriptional gene regulatory networks between taxonomically related organisms[END_REF]. Usual computational approaches focus either on: (1) mapping TFs and their target genes from one genome to another, or (2) using the binding site motif described in one organism to look for binding sites in a new one [START_REF] Venancio | Reconstructing prokaryotic transcriptional regulatory networks: lessons from actinobacteria[END_REF]. Neither of these approaches is perfect since orthology by itself cannot ensure that DNA-binding properties are conserved, or that both proteins respond to the same signals. On the other hand, significant variability has been observed among the binding site motifs for the same TF (e.g. see LexA in [START_REF] Mazón | LexA-binding sequences in Gram-positive and cyanobacteria are closely related[END_REF]) in the different bacterial phylogenetic groups. Finally, TFs can act either as activators or repressors, depending on the precise placement of their binding sites respect to the transcription start site [START_REF] Alberts | Molecular Biology of the Cell[END_REF]. In practice, this strategy requires two steps: (1) orthologous detection for TFs and target genes; as well as (2) new binding site prediction using position weight matrices (PWM) or other motif model. The first step can be accomplished using bidirectional Blast. For the second step there are several alternative tools such as MEME/MAST or PoSSuMsearch. Both utilize existing motif descriptions characterized as PWM and a probabilistic model to determine plausible binding sites. Unfortunately all these models produce a large number of false positives [START_REF] Guhathakurta | Computational identification of transcriptional regulatory elements in DNA sequence[END_REF]. One strategy to reduce the number of false positives is the use of comparative genomics to see if the binding site is conserved among phylogenetically related organisms [START_REF] Xie | Systematic discovery of regulatory motifs in human promoters and 3' utrs by comparison of several mammals[END_REF]. These strategies are known to limit the number of false positives at the cost of increasing the number of false negatives [START_REF] Lepoivre | TranscriptomeBrowser 3.0 : introducing a new compendium of molecular interactions and a new visualization tool for the study of gene regulatory networks[END_REF].

The second strategy uses a set of differential expression results, from the same organism in several stress conditions, to statistically determine the relationship between genes, from the point of view of expression regulation. Mutual information is a statistical concept describing a probabilistic relationship between variables, which can be used to identify and characterize relationships that are not detected by linear correlation [START_REF] Cellucci | Statistical validation of mutual information calculations: comparison of alternative numerical algorithms[END_REF]. In this case these variables are gene expression levels measured along several stress conditions. If two genes have a related behavior -for example if one goes up every time another goes down-then knowing the expression of one gives some information about the other. This is evaluated by mutual information. These relationships are not necessarily causal, for example if genes A and B are regulated by a third one C, mutual information will be significative between A and B as well as between A and C [START_REF] Hartemink | Reverse engineering gene regulatory networks[END_REF]. High mutual information is a necessary but not sufficient condition for causality. Once mutual information is evaluated, there are also several strategies to determine which are the mutual information values that correspond to significant gene relationships, and to prune the resulting graph. The general idea is to keep only the direct causal links and to drop the edges which can be better explained by a third gene. Among these strategies are Relevance Networks, ARACNe, CLR, MRNET and C3NET [START_REF] Altay | Inferring the conservative causal core of gene regulatory networks[END_REF], which is the method used in this work.

In this paper, we propose an hybrid strategy to reduce the size of the set of putative TFs and binding sites using mutual information relationships derived from differential expression experimental data coded as rules in a logical program. The result is an unoriented graph which may have a small intersection with the transcriptional regulation network, because transcription factors can act at very low concentrations, and thus may be below the detection threshold of differential expression experiments. To handle this we consider operons instead of genes as the network nodes. Thus, if a TF regulates a given gene, all the operon where the TF belongs appears to control all the operon of the target gene. This makes sense because, in bacteria, operons are the minimal transcriptional unit [START_REF] Brown | Genomes[END_REF]. This logical program is coded in Answer Set Programming (ASP), a declarative problem solving paradigm in logic programming and knowledge representation, which offers a rich modeling language along with highly efficient inference engines based on Boolean constraint solving technology. In ASP, a problem encoding is a set of logic programming rules which are first transformed into an equivalent propositional logic program and then processed by an answer set solver, which searches for specific solutions to the rules, called Answer Sets. ASP allows solving search problems of high complexity [START_REF]Knowledge Representation, Reasoning and Declarative Problem Solving[END_REF].

In summary we have two sources of data for regulation network: (1) the putative TFs and their possible binding sites, and (2) behavioral relationships stated by mutual information. Complementing these with operon prediction, we ask two questions:

-Can each of the mutual information relationships be explained by a shared regulatory element controlling both genes?

-If so, which is the minimal core of TFs and binding sites required to explain all mutual information relationships?

The plan of this work is the following. First, we reconstructed a genetic regulatory network for Acidithiobacillus ferrooxidans based in the publicly available sequences. Second, we estimate mutual information from microarray experiments in several stress conditions. Finally, we encode both information (based on regulation and experiments) in Answer Set Programming (ASP) to find the minimal set of TFs and binding site families that explains the observed correlations. This procedure optimizes over the complete space of possible configurations and defines a core regulatory network whose properties are further analyzed.

Materials

Genomic DNA sequences were downloaded from NCBI with accession number NC 011761. The functional annotation corresponds to [START_REF] Valdés | Acidithiobacillus ferrooxidans metabolism: from genome sequence to industrial applications[END_REF]. TFs where determined as the best Blast hit with e-value under 1E-5 against Prodoric database of proteins [START_REF] Grote | PRODORIC (release 2009): a database and tool platform for the analysis of gene regulation in prokaryotes[END_REF]. Binding sites where located using FIMO from the MEME/MAST suite [START_REF] Bailey | MEME suite: tools for motif discovery and searching[END_REF][START_REF] Grant | FIMO: scanning for occurrences of a given motif[END_REF] and motifs from Prodoric database, represented as position weight matrices. Gene grouping in operons is given by the prediction of the recent database ProOpDB [START_REF] Taboada | ProOpDB: Prokaryotic operon database[END_REF].

Microarray slides where printed with DNA segments from a native Acidithiobacillus ferrooxidans strain (Wenelen DSM16786). Genomic DNA was shotgun by sonication and 5 568 segments of 2Kbp nominal size where sequenced by both ends and printed in duplicate [START_REF] Levicán | Comparative genomic analysis of carbon and nitrogen assimilation mechanisms in three indigenous bioleaching bacteria: predictions and validations[END_REF]. Acidithiobacillus ferrooxidans grows naturally in ferric or sulfuric medium. The set of experiments evaluated expression in ferric medium in the green channel and compared it versus: (i) sulfur medium, (ii) shift to sulfur, that is, ferric medium with last minut addition of sulfur, (iii) shift to Chalcopyrite (CuFeS 2 ), (iv) shift to Pyrite (FeS 2 ), (v) shift to coveline (CuS), (vi) shift to raw mine ore, and (vii) shift to quartz (SiO 2 ) in the red channel. Cell culture, RNA extraction, hybridization and scanning were performed by BioSigma S.A. (Colina, Chile) as described in [START_REF] Levicán | Comparative genomic analysis of carbon and nitrogen assimilation mechanisms in three indigenous bioleaching bacteria: predictions and validations[END_REF].

Spot quality assessment was performed using several indices summarized in a Qcom value following [START_REF] Wang | Quantitative quality control in microarray image processing and data acquisition[END_REF], which is further used as spot weight. Expression was normalized intra-slide using Lowess regression and interslides using Gquantile, since green channel was the same condition for all 7 experiments considered [START_REF] Smyth | limma: Linear Models for Microarray Data User's Guide[END_REF]. Each DNA segment expression value was calculated as the weighted average of the spots representing it. Since each spot contains DNA segments which can include several genes, we developed a protocol to determine each gene contribution to total spot luminescence. Using Blast each spot sequence was mapped to ATCC23270 genome. To recover the individual gene expression we solve the minimization

g = arg min || c -M g|| 2 s.t. g ≥ 0. ( 1
)
where c is the vector of observed spots luminescence, g has components representing the equivalent luminescence (and thus concentration) of each gene, and M is a matrix whose element M i,j is the length of the intersection between spot i sequence and gene j in the Blast alignment.

This transformation is applied to each channel in each slide. The resulting experiment set was analyzed using the standard procedure implemented in Limma framework [START_REF] Smyth | Limma: linear models for microarray data[END_REF] for the R statistical package [START_REF] Irizarry | R: A Language and Environment for Statistical Computing[END_REF]. A linear model was fitted to each gene using lmFit method and eBayes was used to estimate statistical significance of observed values [START_REF] Smyth | Linear models and empirical bayes methods for assessing differential expression in microarray experiments[END_REF].

Differential expression data from 56 microarrays corresponding to 7 stress conditions was used to estimate the mutual information between genes using c3net library in R [START_REF] Altay | Inferring the conservative causal core of gene regulatory networks[END_REF]. ASP code was executed using the Potassco implementation [START_REF] Gebser | Potassco: The Potsdam answer set solving collection[END_REF].

Results

In this section we detail the results of network reconstruction based on sequence data and its further pruning using logical programming to constraint the network to experimental evidence from differential expression results. Our main result is a reduced network, consistent with the experimental data, with a lower number of spurious regulatory elements and with a size more amenable to be analyzed and simulated using classical tools. Transcription regulation networks can be represented as bipartite graphs, with genes and BS motifs as nodes. Our first result was obtained using MAST and motifs from Prodoric database to find 182 motifs represented in 27 435 putative binding sites upstream of 3 143 genes in A. ferrooxidans. Using reference protein sequences from the same database, we found by homology 66 genes whose product could encode TFs which can bind in 36 of the 182 motifs. This first network is summarized in the first column of TABLE 1.

Mutual information was calculated from 56 microarrays representing 7 stress conditions as described in Materials, and then filtered using the Conserved Causal Core Network protocol (C3NET, [START_REF] Altay | Inferring the conservative causal core of gene regulatory networks[END_REF]). This protocol discards the weakest mutual information links and keeps a spanning-tree with 1 245 edges between genes whose behavior is experimentally related and constitute restrictions that any network simplification must satisfy.

A second network, also summarized in TABLE 1, was obtained simplifying the first one by grouping the 3 143 genes into 1 542 predicted operons, as predicted in ProOpDB. The 27 435 edges connecting motifs to genes are replaced by 22 011 edges, while the 1 245 mutual information relations are reduced to 1 206.

Our third result uses our method to express the network and minimize the number of relations subject to the mutual information restrictions. We encode biological constraints as disjunctive rules that can be processed by ASP, that is as a finite set of rules of the form a 1 ; . . . ; a l ← a l+1 , . . . , a m , not a m+1 , . . . , not a n where a n are atoms. Intuitively, atoms can be viewed as facts and rules as deductions to determine new facts. Rules shall be read from right to left: at least one fact in the left part a 1 ; . . . ; a l (called "head") shall be true whenever all facts in the right part (called "body") a l+1 , . . . , a m , not a m+1 , . . . , not a n are satisfied. Consequently, the rule with empty head ← a. means that the fact a is always false. The answers set of a logical program is a set of atoms that satisfy all the logical rules, together with minimality and stability properties, ensuring that every atom appears in at least one rule.

A strong feature that we use in the following is that atoms can be stated as predicates describing relationships between variables. For instance, upstream(M, G) means that the motif M is present in the regulatory region of the gene G; canBind(G, M ) means that gene G codes for a transcription factor which can bind to motif M ; and miRel(G, G ′ ) represent that genes G and G ′ are related by mutual information, that is, this edge is conserved by C3NET protocol. The predicate inOperon(G, O) is used to group a set of genes G into an operon O. Several genes can be part of the same operon, but no gene can be in two operons.

The first ASP program takes each gene pair related by mutual information and looks for shared motifs either in their upstream region or in the respective promoting region for operons which can regulate indirectly. Since This graph represents the connected component of the union of all optimal solutions, being its intersection represented by yellow nodes and dashed edges. Different operons having a common predecessor in the graph contain genes which are correlated in mutual information outputs. To facilitate visualization we show only the operons involved in at least one of the optimal solutions and omit "leaf" operons, that is, those that not control other operons. mutual information is an index of a related behavior, we expect that each mutual-information-related gene pair should share at least one motif in the sense previously described. We found only one case: AFE 0560 and AFE 2588 are two hypothetical proteins, each in a monocystronic operon, related by mutual information but no common motif was found. This will be further examined, but in a first instance we ignore them.

Leaving this particular case aside we can restrict the model to forbid mutual information relationships which can not be explained. Since several TFs are represented by many genes, we explore the combinatorial space looking for the smallest subsets of TFs satisfying the constraints given by mutual information. This optimization procedure, made feasible by ASP code, also reduces the number of binding site motif used. For each gene G that can bind we generate an eventual bound instance

{bound(G, M )} ← canBind(G, M ).
In ASP syntax this means that for each gene G which can bind to motif M , there will be at most one instance of bound(G, M ), or none. This is the key part of the program, as this generates all the combinatorial possibilities.

Then a motif M controls an operon O either if it is directly upstream it or if there is a transcript factor G mediating

controles(M, O) ← upstream(M, G), inOperon(O, G). controles(M, O) ← controles(M, O ′ ), inOperon(O ′ , G ′ ), bound(G ′ , M ′ ), upstream(M ′ , G), inOperon(O, G).
Now a mutual information link is explained when there is a motif M that can control both

explained(O, O ′ ) ← controles(M, O), controles(M, O ′ ), inOperon(O, A), inOperon(O ′ , B), miRel(A, B).
So the main restriction is that all mutual information links must be explained, or equivalently, no mutual information link must remain unexplained

← miRel(A, B), inOperon(O, A), inOperon(O ′ , B), not explained(O, O ′ ).
Under this ruleset we minimize the number of bound(G, M ) predicates, that is, we look for minimization of the number of TFs used, which can explain the observed mutual information relationships, in particular when several genes could encode TFs binding to the same motif. The minimal number of TFs achieved is 7, which can be realized by 97 different configurations. In the union of all 97 optimal solutions there are 14 motifs (see TABLE 2) and 27 TFs (see TABLE 3), forming a connected subgraph of the initial regulation network, shown in FIG. 1. All operons can be regulated by one of these TFs through binding sites matching one of these motifs. This intersection of all optimal solutions is represented by dashed edges.

As seen in TABLE 3, in all 97 solutions the TF PyrR binds to motif MX000061 and regulates the operon op0404, while FleQ binds always to motif MX000104 regulating three operons. Finally, NnrR binds always to motif MX000192 which regulates gene AFE 2696 (op1250) and operons op1040 and op1543 (which contains gene dnaA). At the same time, the motif MX000140 (LexA) was reported in 55 of the 97 instances to regulate gene lexA (AFE 1868), and in the other 42 instances regulates gene umuD (AFE 1750). Both genes encode putative LexA TFs.

The presence in the solution set of several sigma dependent TFs is consistent with the literature, which describes them as part of the basic transcription modification mechanisms. It is worth noticing that the housekeeping sigma-70 appears in only one solution, meaning that it is not involved in the regulations triggered by the experimental conditions. On the other side sigma-54, associated with nitrogen regulation, appears twice, as encoded by gene AFE 2696 in all solutions, and alternatively by genes AFE 0447, AFE 0957, AFE 174, AFE 2955 and AFE 3025. There are also other TFs related to nitrogen assimilation, as AFE 0024 and AFE 1527. This suggests that sigma-54 plays an important role in the adaptation from an pure iron energy source to a medium including sulfur as an energy source.

An enrichment analysis (data not shown) of GO categories associated to all genes in the operons of FIG. 1 shows significative presence of transcription factors, iron-sulfur cluster assembly and lypopolysaccharide transport, which has been related to the bioleaching activity of A. ferrooxidans.

Conclusions

We have presented an hybrid method for regulatory network reconstruction, which combines probabilistic tools for identification of regulatory elements and experimental relationships, with ASP logical programming to identify the core regulation elements consistent with experimental evidence considering the whole space of possible configurations between regulators and correlated genes. This last point is essential when comparing our approach with other integrative approaches, such as DISTILLER [START_REF] Lemmens | DISTILLER: a data integration framework to reveal condition dependency of complex regulons in Escherichia coli[END_REF] which cannot afford exploring the whole space of configurations and therefore introduces additional criterias of minimum number of correlated genes and regulators in the motifs that could generate biased results. 3. Genes encoding the transcription factors used in any of the optimal solutions. These are the elements which can explain all mutual information correlations in experimental data. A Gene putatively encoding a TF is contained in an Operon (as predicted in ProOpDB), its Name and Description are taken from NCBI's annotation. Column Count shows the number of optimal solutions in which each transcription factor is used, the ones marked in boldface are used in all solutions. Note that several sigma-54 dependent transcriptional regulators are included, which suggest that this TF plays an important role in adaptation to different energy sources.

We applied this method to A. ferrooxidans, an industrial relevant bacteria whose characteristics difficult classical experimental methods. The proposed network drastically reduces the number of nodes necessary to explain co-expressions in differential expression results, and thus can be a significant contribution to the understanding of the biology of this microorganism.

Our perspective is to improve this network by incorporating signs to the edges, either from experimental data or by homology to other organisms, to get a model which can be used to predict behavior under new experimental conditions. Transcriptional units other than operons can also be considered in an improved version of this model. Those predictions will be eventually validated experimentally.

Figure 1 .

 1 Figure1. Core regulatory network of A. ferrooxidans. Nodes stand for operons, while colored edges stand for different transcription factors transcribed by the operons. This graph represents the connected component of the union of all optimal solutions, being its intersection represented by yellow nodes and dashed edges. Different operons having a common predecessor in the graph contain genes which are correlated in mutual information outputs. To facilitate visualization we show only the operons involved in at least one of the optimal solutions and omit "leaf" operons, that is, those that not control other operons.

Table 1 .

 1 Summary of results for each stage of network reconstruction. Regulatory networks are bipartite graphs, connecting Nodes 1 (genes or operons) to Nodes 2 (binding site motifs). The first network is the full reconstruction linking transcription factor (determined by homology to known transcription factor sequences in Prodoric database) to motifs (found with MEME/FIMO using PWM profiles in Prodoric database) and back to genes (from official annotation in NCBI). Edges 1→2 derive from literature relating transcription factor and motifs, while Edges 2→1 exists when a motif has a instance in the promoter region of a gene. Restrictions were determined from mutual information calculated on experimental data as described in the text. Network 2 consolidates genes into operons (as predicted in ProOpDB), reducing network size based on biological hypothesis. The final network is the result of the method we present in the text. The size of the final network, supporting the same experimental evidence, is much smaller than the original one, TF are reduced by one third, BS are reduced to 8% of the initial number, as marked in boldface.

	Element	Network 1	Network 2	Network 3
		(genes)	(operons)	(ASP optimization)
	Nodes 1	3143	1542	1542
	Nodes 2	182	182	14
	Edges 1→2	66	66	27
	Edges 2→1	27435	22011	2428
	MI Restrictions	1245	1206	1206

Table 2 .

 2 Binding Site Motifs used in any of the optimal solutions. All mutual information relationships can be explained using just 7 TFs at a time, in what we call an optimal solution. There are 97 different optimal solutions, which use only these motifs. The four motifs marked in boldface are used in all solutions, suggesting they play a key role in the transcriptional regulation. Motif and protein identifiers were taken from Prodoric Database.

	Motif	Name	Protein Description

Table

	Gene	Operon	Name	Description	Count
	AFE 0024	op0008	ntrX	nitrogen assimilation regulatory protein NtrX, putative	11
	AFE 0119	op0055	-	transcriptional regulator, Crp/Fnr family	97
	AFE 0185	op0087	-	type IV fimbriae expression regulatory protein PilR, putative	14
	AFE 0447	op0170	-	sigma-54 dependent DNA-binding response regulator	11
	AFE 0470	op0180	rpoS	RNA polymerase sigma-38 factor	1
	AFE 0672	op0282	iscR	iron-sulfur cluster assembly transcription factor IscR	20
	AFE 0857	op0367	-	transcriptional regulator, Crp/Fnr family	12
	AFE 0957	op0404	-	sigma-54 dependent transcriptional regulator	11
	AFE 1434	op0640	phoB	DNA-binding response regulator PhoB	9
	AFE 1527	op0668	nifA	Nif-specific regulatory protein	11
	AFE 1658	op0734	-	transcriptional regulator, Crp/Fnr family	65
	AFE 1741	op0780	-	sigma-54 dependent transcriptional regulator	11
	AFE 1750	op0787	umuD	UmuD protein	42
	AFE 1868	op0848	lexA	LexA repressor (EC:3.4.21.88)	55
	AFE 1901	op0864	-	hypothetical protein	12
	AFE 1990	op0900	-	transcriptional regulator, AsnC family	5
	AFE 2271	op1040	-	transcriptional regulator, putative	9
	AFE 2342	op1074	-	RNA polymerase sigm-70 factor family	1
	AFE 2696	op1250	-	sigma-54 dependent transcriptional regulator	97
	AFE 2750	op1276	rpoH	RNA polymerase sigma-32 factor	3
	AFE 2798	op1307	ogt-1	methylated-DNA-protein-cysteine methyltransferase (EC:2.1.1.63) 16
	AFE 2934	op1383	-	transcriptional regulator, putative	3
	AFE 2955	op1394	-	sigma-54 dependent DNA-binding response regulator, putative	8
	AFE 3025	op1420	rpoN	RNA polymerase sigma-54 factor	2
	AFE 3060	op1427	pyrR	PyrR bifunctional protein (EC:2.4.2.9)	97
	AFE 3137	op1469	-	DNA-binding response regulator	27
	AFE 3309	op1543	dnaA	chromosomal replication initiator protein DnaA	29

Acknowledgements

This work was funded by Center for Genome Regulation (Fondap 15090007), Universidad de Chile; Laboratory of Bioinformatics and Mathematics of Genome, Center for Mathematical Modeling, UMI 2807 CNRS-Universidad de Chile; INRIA-U. de Chile IntegrativeBioChile Associate Team; a mobility grant from the International College Doctoral (IDC) of Université Européenne de Bretagne (UEB), INRIA-Conicyt 2010 mobility grant 2010-55 and "Estadías cortas de Investigación para Estudiantes de Doctorado de la Universidad de Chile" grant.


Adam Houenou
email: adam.houenou@hds.utc.fr

Philippe Bonnifait

Veronique Cherfaoui

Jean-François Boissou

A Track-To-Track Association Method for Automotive Perception Systems

Keywords: Multisensor track-to-track association, track level fusion, target tracking, processing architecture

Recent and future driver assistance systems use more and more sensors, that have individual tracking modules. For target tracking, it becomes necessary to find techniques to manage as simply as possible the use of a great number of independent and heterogeneous sensors, at the different stages of the process. This paper presents a modular highlevel track-fusion architecture for a multisensor environment. This architecture allows the variation of the number and the types of the used sensors with no major change in the tracking algorithm. The paper also tackles the multisensor track-to-track association issue with a new algorithm based on a particular track-to-track distance computation. An example of target tracking method is shown to make use of the proposed architecture and the track-to-track association algorithm.

I. INTRODUCTION

Target tracking is a basic task for a number of Advanced Driver Assistance Systems (ADAS). It allows the system to get a representation of the vehicle's environment, containing potential obstacles with their kinematic parameters. This dynamic representation can then be used by different decision modules, depending on the considered functionality. Recent ADAS use more and more sensors (Fig 1) in order to increase confidence in detections, to enlarge their field of view (f.o.v.) and to get more complementary information. Many sensors are now equipped with an individual processing module that performs a local tracking task and reports a list of targets identified through time. Target tracking for a system that uses such sensors consists in fusing the received sensor-tracks in order to get upper level tracks, called system-tracks in this paper. It is a track-level data fusion. A discussion of the advantages and disadvantages of this architecture is presented in chapter 9 of [START_REF] Blackman | Design and analysis of modern tracking systems[END_REF]. The sensors work independently and the high level of abstraction of the data they provide allow the implementation of a tracking process that is not affected by the use of heterogeneous sensors.

In a monosensor target tracking process, the data association step is usually the first one; It consists in matching the new reported observations to the current tracks. In a multisensor environment, one must match the sensor reports that correspond to the same target in order to estimate the target's state by merging the different sensor reports. There Figure 1. Vehicles are equipped with more and more sensors are many data association algorithms such as Nearest Neighbor (NN), Probabilistic Data Association, Multi-Hypothesis Tracking (MHT) and their variants (see [START_REF] Blackman | Design and analysis of modern tracking systems[END_REF]) but most of them are only known for the two-data-classe (current tracks and reported observations) association cases. The cases with more data classes considered simultaneously are less frequently encountered. Data association errors can generate a number of problems such as ghost track, target loss and track ID permutation. Errors usually occur when the targets are close to each other (causing ambiguities) or because of sensors' lack of accuracy.

This paper aims at presenting a track-to-track association algorithm that can be used for as many sensors as wanted. This algorithm is used in a multisensor tracking system presented thereafter. The processing architecture that has been used makes it possible to vary the number and the types of the sensors with little modification in the structure of the tracking system. In section II, we present the tracking architecture and the internal tasks of the tracking system. In section III, we discuss the data synchronization issue in a multisensor environment. In section IV, the data association algorithm is described in details and in section V, we show some simulation results.

II. GLOBAL ARCHITECTURE OF THE TRACKING SYSTEM

This system is intended to be used in vehicles that may have different types and different number of sensor devices (lidars, radars, cameras). So it has been designed to be as sensor-independent as possible. The global processing architecture (Fig. 2) is inspired from the one used in [START_REF] Darms | Obstacle detection and tracking for the urban challenge[END_REF] for the vehicle that won the 2007 Darpa Urban Challenge. It separates the tasks into two main layers, a Sensor Layer (SL) and a Fusion Layer (FL). The common idea is to perform the sensor-specific tasks in the SL and to provide the FL with sensor-independent data so that there would be no modification in the tracking algorithm even if the number and the types of the sensors were modified (because of sensor failure, addition or removal). In our case, there are less interactions between the layers and the internal tasks of the layers are different.

The SL has the role of collecting sensor-specific data from each physical sensor and converting them into a common sensor-independent data type. These data are then sent to the FL that performs the target tracking.

A. Sensor Layer

The perception of the environment is done by means of several sensors that work independently and provide their reports at different rates. Since each sensor has its own internal tracking module the reports are lists of target tracks identified through time. We assume that a track is defined by a timestamp, an identification number, a target state estimate and its error covariance matrix. The target state is composed of different parameters depending on the sensor specifications. In this paper, we will focus only on the tracking of the position and velocity which are usually the minimal reported parameters. A sensor report is relative to the sensor's individual mobile coordinate system and time base. In order to make a global tracking on the reunion of the individual f.o.v., it is necessary to convert the collected sensor data to a common mobile coordinate system and a common time base. These conversions are respectively the spatial alignment and the temporal alignment. Sensors that have the same data types are managed by the same Sensor Processing Module that can be run in parallel instances as new sensor reports are being received (see Fig2). If a sensor with a new data type is added, it will be necessary to add a new Sensor Processing Module. The Sensor Processing Modules perform the following tasks:

-Feature extraction: Sensors reports are usually received in many packages that first need be reassembled in order to form utilizable data structure. The reassembling rule is specific to each sensor. One must check the sensor's documentation. -Spatial alignment: For automotive applications, it is generally sufficient to have a 2D positioning since one can assume that the road is locally plane all around the vehicle and that there is no flying or underground object. Therfore, the sensors' measures are supposed to be made in a plane considered to be the road plane. The common coordinate system is also defined in that plane. A rigid planar transformation can then be used to convert measures from an individual sensor coordinate system to the common coordinate system. Indeed, the relative position of the common coordinate system's origin and orientation in each sensor coordinate system is supposed to be known (sensor configuration and physical mounting point). -Temporal alignment: The timestamp given to a report by a given sensor originates from the sensor's internal clock. In practice, the sensor clocks have different times (due to unknown sensor booting time) but it is possible to estimate the relative offset between them by using the global time (coming from the FL's clock). The individual timestamps can then be converted into global time. Because of clock drifts, the temporal alignment is periodically repeated. . -Data validation: Invalid data are not transmitted to the FL. The validity can be checked by known sensor properties (validity flag or check-sum value). -Sensor meta-data: In complement of the sensors estimates, other information are useful at the FL level. For instance, the identification of the sensors is used for data association and data synchronization, the f.o.v. are used to check overlapping areas for data association. Here is an example of sensor-independent data for a given track: {Timestamp, Sensor meta-data, Track ID, State estimate (position and velocity), Estimation error covariance matrix}.

B. Fusion Layer

The target tracking task itself is performed in the FL. The sensor tracks are asynchronously received from the SL and fused to form system tracks. The fusion process works as follows: First, a software synchronization of the received data is implemented and a fusion rate is defined so that at each fusion time, the system has a report from each sensor. We discuss this choice in section III. The fusion operation consists in associating the reported tracks so that each cluster is composed of tracks originating from the same target. For that, we use a novel Track-To-Track Association (TTTA) algorithm. Often, data association algorithms are used to match track state predictions with sensor observations. In [START_REF] Aziz | Fuzzy track-to-track association and track fusion approach in distributed multisensor multitarget multiple-attribute environment[END_REF] a multisensor TTTA technique based on the computation of degrees of membership was described but the number of clusters (equivalent to the number of targets) has to be known in advance. The TTTA algorithm here, is based on a specific Track-To-Track Distance (TTTD) computation, using buffered past track states. We present it in details in section IV. After the TTTA stage, each cluster stands for one detected target. Depending on the overlapping of the f.o.v. of the sensors, and on their detections, there is one or more reported tracks per cluster. The tracks in each cluster are simultaneously merged to form a meta-track. The merging is done by forming pairs of tracks that are hierarchically merged as explained in Figure 3. Two tracks (X 1 , P 1 ) and (X 2 , P 2 ) are merged using equations ( 1) and ( 2).

X = P 2 • (P 1 + P 2 ) -1 • X 1 + P 1 • (P 1 + P 2 ) -1 • X 2 (1)
P = P 2 • (P 1 + P 2 ) -1 P 1 (2)
This formula (derived from the Kalman filter gain in [START_REF] Smith | On the representation and estimation of spatial uncertainty[END_REF]) is optimal at minimizing the trace of the resulting covariance, under the assumption of uncorrelated measurement noise. Since it is also a commutative and associative fusion operator, fusion can be done in any order. As explained by Bar-Shalom in [START_REF] Bar-Shalom | On the track-to-track correlation problem[END_REF], there may be an inter-sensor correlation due to the temporal correlation of the filters. In our problem, it is difficult to determine this inter-sensor correlation because in practice the characteristics of the filtering performed by each sensor is generally unknown. In [START_REF] Bar-Shalom | Multitarget multisensor tracking: principles and techniques[END_REF] a formula was proposed to approximate it, but only in the case of Kalman filters.

The meta-tracks are then associated to the current systemtracks state predictions using the Global Nearest Neighbor algorithm [START_REF] Konstantinova | A study of target tracking algorithm using global nearest neighbor approach[END_REF]. The prediction is done by using a linear and constant velocity prediction model (described in [START_REF] Bar-Shalom | Tutorial on multisensor management and fusion algorithms for target tracking[END_REF]). The non-associated meta-tracks are used to initialize new systemtracks. The non-associated system-tracks are predicted at each iteration until they match with a meta-track or until they are deleted after a lap of time where the correspondent target is supposed to be definitely lost. The system-tracks that match with meta-tracks are updated by the meta-track's state with no extra filtering.

III. DATA SYNCHRONIZATION AND BUFFER MANAGEMENT

In order to perform a correct data fusion, it is crucial not only to keep coherence on the data chronology but also to fuse data only if they are relative to the exact same instant. Incoherence on data chronology is usually caused by sensors' latency; indeed, because of internal sensor processing time and data transfer time, it may happen that a data is received by the system after others that are more recent yet. This case is called out-of-sequence. Many solutions have been suggested to answer the out-of-sequence issue ( [START_REF] Blackman | Design and analysis of modern tracking systems[END_REF], [START_REF] Kaempchen | Data synchronization strategies for multi-sensor fusion[END_REF], [START_REF] Tessier | Système de localisation basé sur une strategie de perception cognitive appliqué a la navigation autonome d'un robot mobile[END_REF]). Most of them recommend the use of buffers. The one suggested in [START_REF] Tessier | Système de localisation basé sur une strategie de perception cognitive appliqué a la navigation autonome d'un robot mobile[END_REF] consists in fusing the data as they are The fusion of many data at a given instant is non sense if these data are not relative to the same time. Thus, in a multisensor environment where data are received at different rates, there are two possibilities: either an update of the system-tracks is performed every time a new sensor data is received (Fig 4 .a) either, a fusion rate is defined, then, at each fusion time an estimate of each sensor's data is computed and a system update is performed with all the estimates (Fig 4 .b). In the first case, the quality of the update (amount of information gained) depends on the accuracy and the type of information provided by the sensor whose data has just been used. To explain this, let's say that one of the sensors used does not make an estimate of the lateral velocity of targets; after the fusion of its reports, there won't be an update of the tracks lateral velocity. We opted for the second case, where we have steadier updates at steady time. Our solution works as follows: each sensor has a dedicated Sensor Buffer where only its last report is saved, meaning that a new report overwrites the previous one. The fusion module works at the same rate as the fastest sensor, meaning that a fusion task is performed every time a report from this sensor is received. When the fastest sensor sends a new report, a state prediction is made for the other sensors' tracks, using their (buffered) last report. The prediction is performed assuming that the targets have a linear trajectory with a constant velocity between two sensor reports. So, we obtain an apparent synchronization of the sensors. The state predictions and the fastest sensor report are saved in a History Buffer at each fusion time. When an out-of-sequence report is detected, the History Buffer is used to insert that report at "the right moment in the past". Then successive fusion operations are made to catch up the present time, using the following saved reports in the buffer. In order not to use an excessive amount of memory, the History Buffer's size is limited.

IV. TRACK-TO-TRACK ASSOCIATION

TTTA consists in clustering all the currents sensor tracks so that the tracks of a same cluster should be state estimates of the same target. For our multisensor TTTA algorithm, two important and logic assumptions have to be made:

1) in each cluster there is at most one track of each sensor 2) a given track can not be in two different clusters. This means that we assume a sensor won't generate more than one report per target and that two or more targets won't be reported as a single one; which is the ideal case. In practice, that may not be true, due to filtering errors but this algorithm will precisely allow the detection of ghost reports. Indeed, in areas covered by more than one sensor, the ghost reports will be isolated in their cluster. The algorithm is a generalization of the Nearest-Neighbor algorithm to more than two data classes. However it is based on computations of two-by-two distances between tracks.

A. Track-To-Track Distance computation

A track corresponds to the state estimate of a single target identified through time. We can then introduce the notion of track history. Figure 5 is the representation of three tracks from instant k -n to instant k. The state estimate here is limited to the Cartesian position (x, y) (for the sake of representation). At instant k Track C is closer to Track A than Track B. But considering all the known history (past states) of the three tracks, it is clear that Track B is more similar to Track A than Track C. The idea is to use the track histories in order to make better association. We use a similar formula (Eq 3) as the one shown in [START_REF] Müller | A generic video and radar data fusion system for improved target selection[END_REF] to compute the "distance" between two tracks, using their past estimates.

D (a,b) k = 1 n n-1 i=0 d (a,b) k-i (3)
with n being the history size (number of past reports considered) and

d (a,b) k = X a k -X b k T • P a k + P b k -1 • X a k -X b k + ln P a k + P b k
where a and b are two tracks and X k and P k are respectively the state estimate and the error covariance matrix of the corresponding state estimate at instant k. The value assigned to n is equivalent to a lap of time that must not be greater than the more recent track's age.

With this formula, the resulting distance between Track C and Track A will be bigger than the distance between Track B and Track A and the association error will be avoided at instant k.

B. Multisensor Track-To-Track Association algorithm

The multisensor TTTA process is described in algorithm 1.

Iteratively finding the minimal distance guaranties that every new association made is the best possible at each iteration and that the first associations are better than the followings. That is why, at stage 4, if both tracks already have a cluster, nothing is changed. The value given to the threshold (stage 3.e) depends on the sensors' accuracy. When the sensors are accurate, their reports are close to each other because they are close to the actual target state. In that case, the threshold can be defined as small as possible. If they are inaccurate, reports about a same target may be far from each other, so the threshold needs to be enlarge. The size of the history in the Track-To-Track Distance (TTTD) formula (Eq. 3) is the second parameter of the algorithm. When it is defined big, there are less chances for association errors because more past track states are taken into account but the computational requirement becomes bigger.

If at stage 4 there is a multiple minimum value (which should be rare), each hypothesis, in accordance with the first considered cell, must be evaluated; then if the results are different, the chosen hypothesis is the one where the sum of the TTTD of tracks in the same clusters is the smallest.

An example of execution of the algorithm is shown in Appendix.

V. SIMULATION RESULTS

In order to evaluate the performance of this tracking method in various situations, we designed a simulator. Real automotive scenarios can be simulated and we obtain the position and the velocity of the vehicles (targets and host) during the simulation, with a chosen sample rate, in a fixed global Cartesian coordinate system. These are considered as reference data. Sensor measurements are also elaborated. It is possible to simulate several kinds of sensors by adding to the reference data a uniform white noise with a variable amplitude w.r.t. the sensor's basic accuracy and the target's range. Target visibility is simulated by checking sensors' f.o.v. and alignment of objects w.r.t. the sensor's position. Then, the tracking module of each sensor is simulated, by performing a Kalman filter with a linear and constant velocity prediction model (see [START_REF] Bar-Shalom | Tutorial on multisensor management and fusion algorithms for target tracking[END_REF]). A virtual sensor's report at time t is a timestamped list of its tracks. For simplification purpose, the global tracking is done in a fixed Cartesian coordinate system.

We created an urban driving context with four target cars, containing straight road parts, a roundabout and an intersection. The host vehicle is equipped with three sensors, S1, S2 and S3, all in front view sensing. The sensors' f.o.v. are made quite redundant in order to evaluate the performance of the TTTA algorithm in various situations and to compare the global tracking to the individual sensor trackings. The simulation has been run several times with different sensor accuracies as reported in table I. The values in the table are the maximal measurement error rates, expressed in percentage, for each state parameter, for a target at 100m range. When the target is closer, the maximal measurement error rate is smaller. The actual measurement error is a random value in the correspondent interval. One can see that from setting 1 to setting 3 sensor S1 does not change while sensors S2 and S3 are more and more degraded. These settings are used to evaluate the TTTA algorithm. Setting 4 is used to evaluate the result of the global tracking.

A. TTTA Algorithm

The algorithm has two parameters that need to be tuned: the gating threshold (see stage 3.e of the algorithm) and the size of the history considered for the TTTD computation (see eq. 3). We study here the impact of each parameter on the result in terms of erroneous associations. The simulation has been run 5 times per sensor setting per parameter and per chosen parameter value. For a complete execution of the scenario, the TTTA algorithm is run about 400 times (at each sample time). The association errors were of two types: 1) Some tracks permuted their clusters: This can happen when the correspondent targets evolve very closely from each other for a lap of time greater than the defined history size (see Fig. 6a). The advantage of TTTD is to help distinguishing between tracks that are close only for a moment and the ones that have always been quite close because they probably originate from the same target. If the track history is not long enough, target that move closely for a long time may be mistaken due to measurement noise. A solution is to enlarge the history size but a trade-off has to be made between this issue and the computational requirement. Indeed adding one unit to the history size (equivalent to the simulation period) increases the execution time by approximately 1% in our case (three sensors and four targets). It gets worse if the numbers of sensors and targets are greater. 2) There were more clusters than actual targets: it happens when a sensor reports ghost target or when we have erroneous target state estimates that fall out of the association gate (see Fig. 6b). The solution is to enlarge the threshold but that increases the computational requirement and can lead to associating tracks that are not similar at all. A third possible type of errors, that did not occur during the tests is to have less clusters than real targets. If the gating threshold is defined too large and in case of miss detections, reports originating from 2 different targets may be clustered 6c). However, this case should be rare because it needs each sensor to have a miss detection; indeed if at least one does not miss a detection, the number of cluster will be at least equal to the number of targets.

Table II.a reports the percentage of erroneous results w.r.t. the threshold value. The history size is fixed to 10 here. One can notice that the more inaccurate sensors are (as shown by settings 2 and 3) the larger the threshold should be chosen to have a good result. An expression of the relationship between the sensor accuracies and the required threshold cannot be clearly defined because the statistical distribution followed by the TTTD (Eq. 3) is not known. The right threshold should be chosen experimentally.

Table II.b reports the percentage of erroneous results as a function of the history size. The gating threshold is fixed to 30. One can notice that taking into account more past sensor reports helps to solve ambiguous situations caused by temporary closeness of some targets and leads to less erroneous associations.

Whatever the sensor accuracies, it is possible to parametrize the algorithm so that it performs well.

B. Fusion result

The efficiency of the global tracking has also been evaluated. Here we use setting 4 described in Table I, for the sensor accuracies. S1 is accurate for target velocity estimation, S3 for target localization and S2 is less accurate for both. The gating threshold is set to 30 and the the history size equals 10. Figure 7 shows the longitudinal position and speed errors for the tracking of a chosen target. The colored curves are the result for each sensor tracking and the black curve is for the global tracking. One can notice that for each state parameter the global tracking result is always close to the estimate provided by the best sensor even if it is not always better than it. This is because in the merging formula (equations 1 and 2), sensors with a large error covariance matrix have a low (but not null) weight. Nevertheless, globally speaking (for all state parameters), the global tracking is better than each sensor's tracking.

VI. CONCLUSION AND FUTURE WORK

In this article, track-level multisensor target tracking has been studied for automotive applications. We have proposed a processing architecture that is robust to the addition, removal or failure of sensors. The use of many sensors has a lot of advantages but one can take benefit on it only if correct associations between the reported sensor tracks is made. For this we have proposed a multisensor TTTA method that uses track histories. This approach has several advantages like the ability to reduce ambiguities and to partially cover the lack of accuracy of sensors. We have also reminded a buffer management that solves the out-of-sequence and the data synchronization issues.

We focused only on the position and the velocity tracking but future work include other attributes such as object dimensions and classification. We are also working on global track existence probability management.

Figure 2 .

 2 Figure 2. Processing architecture

Figure 3 .

 3 Figure 3. Track merging architecture.

Figure 4 .

 4 Figure 4. Data synchronization. a) Asynchronous fusion. b) Software synchronizationbeing received (regardless of the fusion technique used) and keeping them in a buffer for a limited lap of time. When an out-of-sequence case is detected, the buffer is used to go back into the past and redo the fusion, with the right chronology, up to present time.The fusion of many data at a given instant is non sense if these data are not relative to the same time. Thus, in a multisensor environment where data are received at different rates, there are two possibilities: either an update of the system-tracks is performed every time a new sensor data is received (Fig 4.a) either, a fusion rate is defined, then, at each fusion time an estimate of each sensor's data is computed and a system update is performed with all the estimates (Fig 4.b). In the first case, the quality of the update (amount of information gained) depends on the accuracy and the type of information provided by the sensor whose data has just been used. To explain this, let's say that one of the sensors used does not make an estimate of the lateral velocity of targets; after the fusion of its reports, there won't be an update of the tracks lateral velocity. We opted for the second case, where we have steadier updates at steady time. Our solution works as follows: each sensor has a dedicated Sensor Buffer where only its last report is saved, meaning that a new report overwrites the previous one. The fusion module works at the same rate as the fastest sensor, meaning that a fusion task is performed every time a report from this sensor is received. When the fastest sensor sends a new report, a state prediction is made for the other sensors' tracks, using their (buffered) last report. The prediction is performed assuming that the targets have a linear trajectory with a constant velocity between two sensor reports. So, we obtain an apparent synchronization of the sensors. The state predictions and the fastest sensor report are saved in a History Buffer at each fusion time. When an out-of-sequence report is detected, the History Buffer is used to insert that report at "the right moment in the past". Then successive fusion operations are made to catch up the present time, using the following saved reports in the buffer. In order not to use an excessive amount of memory, the History Buffer's size is limited.

Figure 5 .

 5 Figure 5. The use of track history for a better TTTA

Algorithm 1

 1 Multisensor Track-To-Track Association 1) Collect the tracks of all the sensors 2) Assign a number from 1 to N to each track, N being the total number of tracks. 3) Create a N ×N array for the TTTDs between the tracks a) Set cells over the diagonal to a defined maximal value (M axV al) in order not to compute twice the distance between two same tracks. b) Set cells corresponding to two tracks of the same sensor to M axV al, in order not to associate two tracks of a same sensor. c) Set the remaining cells to the distance between the corresponding two tracks. d) Set cells where the distance is greater than a defined threshold to M axV al. The threshold symbolizes a gate out of which we assume that the two tracks cannot originate from the same target. 4) Loop: Determine the minimal value (M inV al) of the array and its position (lin, col) in the array. While M inV al is smaller than M axV al a) If none of the corresponding two tracks has not been inserted in a cluster yet, then put both of them in a new cluster. b) If only one has already been inserted in a cluster, then add the second to that cluster. c) If both have already been inserted in a cluster (the same or not), then do nothing. d) Set to M axV al cells in the line lin and the ones in the column col, that correspond to tracks reported by the two concerned sensors. 5) Each track that has not been inserted in a cluster forms a new cluster (singletons).

Figure 6 .

 6 Figure 6. Association errors. Targets are represented in black and sensor reports in color. a) Cluster permutation b) More clusters than targets c) Less clusters than targets

Figure 7 .

 7 Figure 7. Global tracking errors versus sensor tracking errors

Figure 8 .

 8 Figure 8. Example of TTTA. T ij are sensor reported positions. Red dots are the actual positions of the targets. Red ellipses are correct track clusters

Table I SENSOR

 I ACCURACY SETTINGS USED DURING THE SIMULATION

Table II PERCENTAGE

 II OF ERRONEOUS ASSOCIATIONS W.R.T. SENSOR ACCURACIES AND A) GATING THRESHOLD B)HISTORY SIZE together (see Fig.

APPENDIX

Let's consider a host vehicle with four sensors, tracking four target vehicles in their overlapped f.o.v. Figure 8 depicts an example of sensor reports where T ij are the positions


Fabrizio Dabbene

Didier Henrion

Minimum volume semialgebraic sets for robust estimation

Motivated by problems of uncertainty propagation and robust estimation we are interested in computing a polynomial sublevel set of fixed degree and minimum volume that contains a given semialgebraic set K. At this level of generality this problem is not tractable, even though it becomes convex e.g. when restricted to nonnegative homogeneous polynomials. Our contribution is to describe and justify a tractable L 1 -norm or trace heuristic for this problem, relying upon hierarchies of linear matrix inequality (LMI) relaxations when K is semialgebraic, and simplifying to linear programming (LP) when K is a collection of samples, a discrete union of points.

Introduction

In this paper, we consider the problem of computing reliable approximations of a given set K ⊂ R n . The set K is assumed to have a complicated shape (e.g. nonconvex, nonconnected), expressed in terms of semialgebraic conditions, and we seek for approximations which should i) be easy computable and ii) have a simple description.

The problem of deriving reliable approximations of overly complicated sets by means of simpler geometrical shapes has a long history, and it arises in many research fields related to optimization, system identification and control.

In particular, outer bounding sets, i.e. approximations that are guaranteed to contain the set K, are widespread in the technical literature, and they find several applications in robust control and filtering. For instance, set-theoretic state estimators for uncertain discrete-time nonlinear dynamic systems have been proposed in [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF][START_REF] Durieu | Multi-input multi-output ellipsoidal state bounding[END_REF][START_REF] Ghaoui | Robust filtering for discrete-time systems with bounded noise and parametric uncertainty[END_REF][START_REF] Shamma | Approximate set-valued observers for nonlinear systems[END_REF]. These strategies adopt a set-membership approach [START_REF] Garulli | Robustness in Identification and Control[END_REF][START_REF] Schweppe | Uncertain Dynamical Systems[END_REF], and construct compact sets that are guaranteed to bound the systems states which are consistent with the measured output and the norm-bounded uncertainty. Outer approximation also arise in the context of robust fault detection problems (e.g., see [START_REF] Ingimundarson | Robust fault detection using zonotope-based set-membership consistency test[END_REF]) and of reachability analysis of nonlinear and/or hybrid systems [START_REF] Hwang | Applications of polytopic approximations of reachable sets to linear dynamic games and a class of nonlinear systems[END_REF][START_REF] Kurzhanski | Ellipsoidal techniques for reachability analysis[END_REF].

Similarly, inner approximations are employed in nonlinear programming [START_REF] Nesterov | Interior Point Polynomial Algorithms in Convex Programming[END_REF], in the solution of design centering problems [START_REF] Wojciechowski | Ellipsoidal method for design centering and yield estimation[END_REF] and for fixed-order controller design [START_REF] Henrion | Convex inner approximations of nonconvex semialgebraic sets applied to fixed-order controller design[END_REF].

Recently, the authors of [START_REF] Dabbene | On the complexity of randomized approximations of nonconvex sets[END_REF] have proposed an approach based on randomization, which constructs convex approximations of generic nonconvex sets which are neither inner nor outer, but they enjoy some specific probabilistic properties. In this context, an approximation is considered to be reliable if it contains "most" of the points in K with prescribed high probability. The key tool in this framework is the generation of random samples inside K, and the construction of a convex set containing these samples.

In all the approaches listed above, several geometric figures have been adopted as approximating sets. The application of ellipsoidal sets to the state estimation problem has been introduced in the pioneering work [START_REF] Schweppe | Uncertain Dynamical Systems[END_REF] and used by many different authors from then on; see, for example, [START_REF] Durieu | Multi-input multi-output ellipsoidal state bounding[END_REF][START_REF] Ghaoui | Robust filtering for discrete-time systems with bounded noise and parametric uncertainty[END_REF]. The use of polyhedrons was proposed in [START_REF] Kuntsevich | Synthesis of optimal and adaptative control systems: The game approach[END_REF] to obtain an increased estimation accuracy, while zonotopes have been also recently studied in [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF][START_REF] Guibas | Zonotopes as bounding volumes[END_REF].

More recent works, like for instance [START_REF] Bertsimas | Optimal data fitting: a moment approach[END_REF][START_REF] Henrion | Convex inner approximations of nonconvex semialgebraic sets applied to fixed-order controller design[END_REF][START_REF] Magnani | Tractable fitting with convex polynomials via sum-of-squares[END_REF], employ sets defined by semialgebraic conditions. In particular, in [START_REF] Magnani | Tractable fitting with convex polynomials via sum-of-squares[END_REF] the authors use polynomial sum-of-squares (SOS) programming, a particular class of SDP, to address the problem of fitting given data with a convex polynomial, seen as a natural extension of quadratic polynomials and ellipsoids. Convexity of the polynomial is ensured by enforcing that its Hessian is matrix SOS, and volume minimisation is indirectly enforced by increasing the curvature of the polynomial. In [START_REF] Bertsimas | Optimal data fitting: a moment approach[END_REF] the authors propose moment relaxations for the separation and covering problems with semialgebraic sets, thereby also extending the classical ellipsoidal sets used in data fitting problems.

Our contribution is to extend further these works to cope with volume minimization of arbitrary (e.g. non-convex, non-connected, higher degree) semialgebraic sets containing a given semialgebraic set (e.g. a union of points). Since there is no analytic formula for the volume of a semialgebraic set, in terms of the coefficients of the polynomials defining the set, we have no hope to solve this optimization problem globally. Instead, we describe and justify analytically and geometrically a computationally tractable heuristic based on L 1 -norm or trace minimization.

Problem statement

Given a compact basic semialgebraic set

K := {x ∈ R n : g i (x) ≥ 0, i = 1, 2, . . . , m}
where g i (x) are real multivariate polynomials, we want to compute a polynomial sublevel set V(q) := {x ∈ R n : q(x) ≤ 1} ⊃ K of minimum volume that contains K. Set V(q) is modeled by a polynomial q belonging to P d , the vector space of multivariate real polynomials of degree less than or equal to d.

In other words, we want to solve the following optimization problem:

inf q∈P d vol V(q) s.t. K ⊂ V(q). (1)
In the above problem vol V(q) :=

V(q) dx = R n I V(q) (x)dx
is the volume or Lebesgue measure of set V(q) and I X (x) is the indicator function, equal to one when x ∈ X and zero otherwise.

Note that, typically, the set K has a complicated description in terms of polynomials g i (e.g. coming from physical measurements and/or estimations) and set V(q) has a simple description (in the sense that the degree of q is small, say less than 10). Minimization of the volume of V(q) means that we want V(q) to capture most of the geometric features of K.

If K is convex and q is quadratic, then the infimum of problem ( 1) is attained, and there is a unique (convex) ellipsoid V(q) of minimum volume that contains K, called Löwner-John ellipsoid. It can be computed by convex optimization, see e.g. [3, §4.9].

In general, without convexity assumptions on K, a solution to problem (1) is not unique.

There is also no guarantee that the computed set V(q) is convex. Optimization problem ( 1) is nonlinear and semi-infinite, in the sense that we optimize over the finite-dimensional vector space P d but subject to an infinite number of constraints (to cope with set inclusions).

If we denote by π d (x) a (column vector) basis of monomials of degree up to d, we can write q(x) = π T d (x)q where q is a vector of coefficients of given size. In vector space P d , the set Q := {q : π T d (x)q ≤ 1} is by definition the polar of the bounded set {π d (x) : x ≤ R} where R > 0 is a constant chosen sufficiently large so that all vectors x ∈ V(q) have norm less than R. As the polar of a compact set whose interior contains the origin, set Q is compact, see e.g. [START_REF] Lay | Convex Sets and Their Applications[END_REF]. It follows that the feasible set of problem (1) is compact, and since the objective function is continuous, the infimum in problem (1) is attained.

Convex conic formulation

In this section it is assumed that q is a nonnegative homogeneous polynomial, or form, of even degree d = 2δ in n variables. Under this restriction, in [START_REF] Lasserre | evel sets and non gaussian integrals of positively homogeneous functions[END_REF]Lemma 2.4] it is proved that the volume function q → vol V(q)

is convex in q. The proof of this statement relies on the striking observation [START_REF] Morosov | New and old results in resultant theory[END_REF] that

vol V(q) = C d R n e -q(x) dx
where C d is a constant depending only on d. Note also that boundedness of V(q) implies that q is nonnegative, since if there is a point x 0 ∈ R n such that q(x 0 ) < 0, and hence

x 0 ∈ V(q), then by homogeneity of q it follows that q(λx 0 ) = λ 2δ q(x 0 ) < 0 for all λ and hence λx 0 ∈ V(q) for all λ which contradicts boundedness of V(q). This implies that problem (1), once restricted to nonnegative forms, is a convex optimization problem.

In [START_REF] Lasserre | evel sets and non gaussian integrals of positively homogeneous functions[END_REF]Lemma 2.4] explicit expressions are given for the first and second order derivatives of the volume function, in terms of the moments R n

x α e -q(x) dx

for α ∈ N n , |α| ≤ 2d. In an iterative algorithm solving convex problem [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF], one should then be able to compute repeatedly and quickly integrals of this kind, arguably a difficult task.

When q is not homogeneous, we do not know under which conditions on q the volume function V(q) is convex in q.

Motivated by these considerations, in the remainder of the paper we propose a simpler approach to solving problem (1) which is not restricted to forms, and which does not require the potentially intricate numerical computation of moments (2). Our approach is however only a heuristic, in the sense that we do not provide guarantees of solving problem (1) globally.

L 1 -norm minimization

Let us write V(q) as a polynomial superlevel set

U(p) := V(q) = {x ∈ R n : p(x) := 2 -q(x) ≥ 1}. with polynomial p(x) = π T δ (x)P π δ (x) (3)
expressed as a quadratic form in a given (column vector) basis π δ (x) of monomials of degree up to δ := d 2 , with symmetric Gram matrix P . Then, optimization problem (1)

reads v * d := min p∈P d vol U(p) s.t. K ⊂ U(p). (4)
Note that in problem (4) we can indifferently optimize over coefficients of p or coefficients of matrix P , since they are related linearly via constraint (3).

Since K is compact by assumption and U(p) is compact for problem (4) to have a finite minimum, we suppose that we are given a compact semialgebraic set

B := {x ∈ R n : b i (x) ≥ 0, i = 1, 2, . . . , m b }
such that U(p) ⊂ B and hence

U(p) = {x ∈ B : p(x) ≥ 1}.
The particular choice of polynomials b i will be specified later on. Now, observe that by definition p(x) ≥ I U(p) (x) on R n and hence, integrating both sides we get

B p(x)dx ≥ B I U(p) (x)dx = vol U(p),
an inequality known as Chebyshev's inequality, widely used in probability, see e.g. [2, §2.4.9]. If polynomial p is nonnegative on B then the above left-hand side is the L 1 -norm of p, and the inequality becomes

p 1 ≥ vol U(p). (5)
Now consider the following L 1 -norm minimization problem: Proof: The graph of polynomial p lies above I K , the indicator function of set K, while being nonnegative on B, so minimizing the L 1 -norm of p on B yields an upper approximation of I K . Monotonicity of the sequence w * d follows immediately since polynomials of degree d + 1 include polynomials of degree d. When its degree increases, p converges in L 1 -norm to I K , hence p 1 converges to vol K. The convergence is pointwise almost everywhere, and almost uniform, but not uniform since I K is discontinuous on B.

w * d := min p∈P d p 1 s.t. p ≥ 0 on B p ≥ 1 on K. (6
Note that this L 1 -norm minimization approach was originally proposed in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] to compute numerically the volume and moments of a semialgebraic set.

Trace minimization

In this section we give a geometric interpretation of problem [START_REF] Durieu | Trace versus determinant in ellipsoidal outerbounding with application to state estimation[END_REF]. First note that the objective function reads

p 1 = B p(x)dx = B π T δ (x)P π δ (x)dx = trace P B π δ (x)π T δ (x)dx = trace P M where M := B π δ (x)π T δ (x)dx
is the matrix of moments of the Lebesgue measure on B in basis π δ (x). In equation ( 3) if the basis is chosen such that its entries are orthonormal w.r.t. the (scalar product induced by the) Lebesgue measure on B, then M is the identity matrix and inequality (5) becomes trace P ≥ vol U(p) [START_REF] Durieu | Multi-input multi-output ellipsoidal state bounding[END_REF] which indicates that, under the above constraints, minimizing the trace of the Gram matrix P entails minimizing the volume of U(p).

The choice of polynomials b i in the definition of the bounding set B should be such that the objective function in problem ( 6) is easy to compute. If

p(x) = π T d (x)p = α p α [π d (x)] α then B p(x)dx = α p α B [π d (x)] α dx = α p α y α
and we should be able to compute easily the moments

y α := B [π d (x)] α dx
of the Lebesgue measure on B w.r.t. basis π d (x). This is the case e.g. if B is a box.

Remark 1 (Minimum trace heuristic for ellipsoids) Note that, in the case of quadratic polynomials, i.e. d = 2, we retrieve the classical trace heuristic used for volume minimization, see e.g. [START_REF] Durieu | Trace versus determinant in ellipsoidal outerbounding with application to state estimation[END_REF].

If B = [-1, 1] n then the basis π 1 (x) = √ 6 2
x is orthonormal w.r.t. the Lebesgue measure on B and p 1 = 3 2 trace P . The constraint that p is nonnegative on B implies that the curvature of the boundary of U(p) is nonnegative, hence that U(p) is convex.

In [START_REF] Magnani | Tractable fitting with convex polynomials via sum-of-squares[END_REF], the authors restricted the search to convex polynomial sublevels U(p) = V(q) = {x : q(x) ≤ 1} by enforcing positive semidefiniteness of the Gram matrix of the quadratic Hessian form of q. They proposed to maximize (the logarithm of) the determinant of the Gram matrix, justifying this choice by explaining that this increases the curvature of the polynomial sublevel set along all directions, and hence minimizes the volume. Supported by Lemma 1 and the above discussion, we came to the consistent conclusion that minimizing the trace of the Gram matrix of p, that is, maximizing the trace of the Gram matrix of q, is a relevant heuristic for volume minimization. Note however that in our approach we do not enforce convexity of U(p).

Handling constraints

If K = {x : g i (x) ≥ 0, i = 1, 2, . . . , m} is a general semialgebraic set then we must ensure that polynomial p -1 is nonegative on K, and for this we can use Putinar's Positivstellensatz and a hierarchy of finite-dimensional convex LMI relaxations which are linear in the coefficients of p. More specifically, we write p -1 = r 0 + i r i g i where r 0 , r 1 , . . . , r m are polynomial sum-of-squares of given degree, to be found. For each fixed degree, the problem of finding such polynomials is an LMI, see e.g. [START_REF] Lasserre | evel sets and non gaussian integrals of positively homogeneous functions[END_REF]Section 3.2]. The constraint that p is nonnegative on B can be ensured similarly.

A particularly interesting case is when K is a discrete set, i.e. a union of points x i ∈ R n , i = 1, . . . , N . Indeed, in this case the inclusion constraint K ⊂ U(p) is equivalent to a finite number of inequalities p(x i ) ≥ 1, i = 1, . . . , N which are linear in the coefficients of p. Similarly, the constraint that p is nonnegative on B can be handled by linear inequalities p(x j ) ≥ 0 enforced at a dense grid of points x j ∈ B, j = 1, . . . , M , for M sufficiently large. Note that, with this pure linear programming (LP) approach, it is not guaranteed that p is nonnegative on B, but what matters primarily is that K ⊂ U(p), which is indeed guaranteed. This purely LP formulation allows to deal with problems with rather large N .

Examples

We illustrate the proposed approach for the case when K is a discrete set, a union of points of R n . For our numerical examples, we have used the YALMIP interface for Matlab to model the LMI optimization problem [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF], and the SDP solver SeDuMi to solve numerically the problem. Since the degrees of the semialgebraic sets we compute are typically low (say less than 20), we did not attempt to use appropriate polynomial bases (e.g. Chebyshev polynomials) to improve the quality and resolution of the optimization problems, see however [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] for a discussion on these numerical matters.

Line (n = 1)

To illustrate the behavior of the proposed optimization procedure, we first consider B = [-1, 1] and K = {-1 2 , 0, 1 4 }. On Figure 1 we represent the solutions p of degrees 2, 7, 17 and 26 of minimization problem [START_REF] Durieu | Trace versus determinant in ellipsoidal outerbounding with application to state estimation[END_REF]. We observe that the superlevel set U(p) = {x ∈ B : p(x) ≥ 1} is simply connected for degree 2, doubly connected for degree 7, and triply connected for degrees 17 and 26. Note that, as the degree increases of p, the length of the intervals for which p(x) > 1 tends to zero. This is consistent with the fact that the volume of a finite set is zero. We also observe that polynomial p shows increasing oscillations for increasing degrees, a typical feature of such discontinuous function approximations.

Plane (n = 2)

In B = [-1, 1] 2 we consider two clouds of 50 points each, i.e. N = 100. On Figure 2 we represent the solutions p of degrees 2, 5, 9 and 14 of minimization problem [START_REF] Durieu | Trace versus determinant in ellipsoidal outerbounding with application to state estimation[END_REF]. Here too we observe that increasing the degree of p allows to disconnect set U(p). The side effects near the border of B on the top right figure can be removed by enlarging the bounding set B.

Space (n = 3)

In B = [-1, 1] 3 we consider N = 10 points. The solutions p of degrees 4, 10, and 14 of minimization problem ( 6) is depicted in Figure 3. Here too we observe that increasing the degree of p allows to capture point clusters in distinct connected components.

It should be pointed out that all the semialgebraic sets in the previous examples were computed in a few seconds of CPU time on a standard PC.

Conclusion

In this paper, we proposed a simple technique for approximating a given set of complicated shape by means of the polynomial sub level set of minimal size that contains it. The proposed approximation has been shown to be convex in the case the polynomial is assumed to be homogeneous. Then, a tractable relaxation based on Chebychev's inequality has been introduced. Interestingly, this relaxation reduces to the classical trace minimization heuristic in the case of quadratic polynomials, thus indirectly providing an intuitive explanation to this widely used criterion.

Ongoing research is devoted to utilize the proposed approximation to construct new classes of set-theoretic filters, in the spirit of the works [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF][START_REF] Ghaoui | Robust filtering for discrete-time systems with bounded noise and parametric uncertainty[END_REF].

Figure 1 :

 1 Figure 1: Minimum L 1 -norm polynomials p (red) such that p ≥ 0 on [-1, 1] and p ≥ 1 at the 3 points (black), for degree 2 (upper left), 7 (upper right), 17 (lower left) and 26 (lower right).

Figure 2 :

 2 Figure 2: Minimum L 1 -norm polynomials p (red) and low surface semialgebraic sets (yellow) such that p ≥ 0 on [-1, 1] 2 and p ≥ 1 at 100 points (black), for degree 2 (upper left), 5 (upper right), 9 (lower left) and 14 (lower right).

Figure 3 :

 3 Figure 3: Including the same 10 space points (black) in low volume semialgebraic sets (yellow) of degree 4 (left), 10 (center), 14 (right).

CNR-IEIIT; Corso Duca degli Abruzzi 123, Torino; Italy. fabrizio.dabbene@polito.it

CNRS, LAAS, 7 avenue du colonel Roche, F-31400 Toulouse; France. henrion@laas.fr

Université de Toulouse, LAAS, F-31400 Toulouse; France

Faculty of Electrical Engineering, Czech Technical University in Prague, Technická 2, CZ-16626 Prague, Czech Republic

Acknowlegdments

This work was partially supported by a Bilateral Project CNRS-CNR.