Julien Stoehr

Pierre Pudlo

Lionel Cucala

Adaptive ABC model choice and geometric summary statistics for hidden Gibbs random fields

Keywords: Approximate Bayesian Computation, model choice, hidden Gibbs random fields, summary statistics, misclassification rate, k-nearest neighbors

) who exhibited sufficient statistics on directly observed Gibbs random fields. But when the random field is latent, the sufficiency falls and we complement the set with geometric summary statistics. The general approach to construct these intuitive statistics relies on a clustering analysis of the sites based on the observed colors and plausible latent graphs. The efficiency of ABC model choice based on these statistics is evaluated via a local error rate which may be of independent interest. As a byproduct we derived an ABC algorithm that adapts the dimension of the summary statistics to the dataset without distorting the model selection.

Introduction

Gibbs random fields are polymorphous statistical models, that are useful to analyse different types of spatially correlated data, with a wide range of applications, including image analysis [START_REF] Hurn | A Tutorial on Image Analysis[END_REF], disease mapping [START_REF] Green | Hidden Markov Models and Disease Mapping[END_REF], genetic analysis [START_REF] François | Bayesian Clustering Using Hidden Markov Random Fields in Spatial Population Genetics[END_REF] among others. The autobinomial model [START_REF] Besag | Spatial interaction and the statistical analysis of lattice systems (with Discussion)[END_REF] which encompasses the Potts model, is used to describe the spatial dependency of discrete random variables (e.g., shades * I3M -UMR CNRS 5149, Université Montpellier II, France of grey or colors) on the vertices of an undirected graph (e.g., a regular grid of pixels). See for example [START_REF] Alfò | A finite mixture model for image segmentation[END_REF] and [START_REF] Moores | Segmentation of cone-beam CT using a hidden Markov random field with informative priors[END_REF] who performed image segmentation with the help of the above modeling. Despite their popularity, these models present major difficulties from the point of view of either parameter estimation [START_REF] Friel | Bayesian Inference in Hidden Markov Random Fields for Binary Data Defined on Large Lattices[END_REF][START_REF] Friel | Bayesian inference for Gibbs random fields using composite likelihoods[END_REF][START_REF] Everitt | Bayesian Parameter Estimation for Latent Markov Random Fields and Social Networks[END_REF] or model choice [START_REF] Grelaud | ABC likelihood-free methods for model choice in Gibbs random fields[END_REF][START_REF] Friel | Evidence and Bayes Factor Estimation for Gibbs Random Fields[END_REF][START_REF] Cucala | Bayesian Inference on a Mixture Model With Spatial Dependence[END_REF], due to an intractable normalizing constant. Remark the exception of small latices on which we can apply the recursive algorithm of [START_REF] Reeves | Efficient recursions for general factorisable models[END_REF], [START_REF] Friel | Recursive computing and simulation-free inference for general factorizable models[END_REF] and obtain a reliable approximation of the normalizing constant. However, the complexity in time of the above algorithm grows exponentially and is thus helpless on large lattices.

The present paper deals with the challenging problem of selecting a dependency structure of an hidden Potts model in the Bayesian paradigm and explores the opportunity of approximate Bayesian computation (ABC, [START_REF] Tavaré | Inferring Coalescence Times From DNA Sequence Data[END_REF][START_REF] Pritchard | Population growth of human Y chromosomes: a study of Y chromosome microsatellites[END_REF], Marin et al., 2012[START_REF] Baragatti | An overview on Approximate Bayesian Computation[END_REF] to answer the question. Up to our knowledge, this important question has not yet been addressed in the Bayesian literature. Alternatively we could have tried to set up a reversible jump Markov chain Monte Carlo, but follows an important work for the statistician to adapt the general scheme, as shown by Caimo andFriel (2011, 2013) in the context of exponential random graph models where the observed data is a graph. [START_REF] Cucala | Bayesian Inference on a Mixture Model With Spatial Dependence[END_REF] addressed the question of inferring the number of latent colors with an ICL criterion but their complex algorithm cannot be extended easily to choose the dependency structure. Other approximate methods have also been tackled in the literature such as pseudo-likelihoods [START_REF] Besag | Statistical Analysis of Non-Lattice Data[END_REF], mean field approximations [START_REF] Forbes | Hidden Markov random field model selection criteria based on mean fieldlike approximations[END_REF] but lacks theoretical support.

Approximate Bayesian computation (ABC) is a simulation based approach that can addresses the model choice issue in the Bayesian paradigm. The algorithm compares the observed data y obs with numerous simulations y through summary statistics S(y) in order to supply a Monte Carlo approximation of the posterior probabilities of each model. The choice of such summary statistics presents major difficulties that have been especially highlighted for model choice [START_REF] Robert | Lack of confidence in approximate Bayesian computation model choice[END_REF][START_REF] Didelot | Likelihood-free estimation of model evidence[END_REF]. Beyond the seldom situations where sufficient statistics exist and are explicitly known (Gibbs random fields are surprising examples, see [START_REF] Grelaud | ABC likelihood-free methods for model choice in Gibbs random fields[END_REF], [START_REF] Marin | Relevant statistics for Bayesian model choice[END_REF] provide conditions which ensure the consistency of ABC model choice. The present work has thus to answer the absence of available sufficient statistics for hidden Potts fields as well as the difficulty (if not the impossibility) to check the above theoretical conditions in practice.

Recent articles have proposed automatic schemes to construct theses statistics (rarely from scratch but based on a large set of candidates) for Bayesian parameter inference and are meticulously reviewed by [START_REF] Blum | A Comparative Review of Dimension Reduction Methods in Approximate Bayesian Computation[END_REF] who compare their performances in concrete examples. But very few has been accomplished in the context of ABC model choice apart from the work of [START_REF] Prangle | Semi-automatic selection of summary statistics for ABC model choice[END_REF]. The statistics S(y) reconstructed by [START_REF] Prangle | Semi-automatic selection of summary statistics for ABC model choice[END_REF] have good theoritical properties (those are the posterior probabilities of the models in competition) but are poorly approximated with a pilot ABC run [START_REF] Robert | Lack of confidence in approximate Bayesian computation model choice[END_REF], which is also time consuming.

The paper is organized as follows: Section 2 presents ABC model choice as a k-nearest neighbor classifier, and defines a local error rate which is the first contribution of the paper. We also provide an adaptive ABC algorithm based on the local error to select automatically the dimension of the summary statistics. The second contribution is the introduction of a general and intuitive approach to produce geometric summary statistics for hidden Potts model in Section 3. We end the paper with numerical results in that framework.

Local error rates and adpative ABC model choice

When dealing with models whose likelihood cannot be computed analytically, Bayesian model choice becomes challenging since the evidence of each model writes as the integral of the likelihood over the prior distribution of the model parameter. ABC provides a method to escape from the intractability problem and relies on many simulated datasets from each model either to learn the model that fits the observed data y obs or to approximate the posterior probabilities.

We refer the reader to reviews on ABC (Marin et al., 2012, Baragatti and[START_REF] Baragatti | An overview on Approximate Bayesian Computation[END_REF] to get a wider presentation and will focus here on the model choice procedure. .

Background on Approximate Bayesian computation for model choice

(1)

When the goal of the Bayesian analysis is the selection of the model that best fits the observed data y obs , it is performed through the maximum a posteriori (MAP) defined by m MAP (y obs ) = arg max m π(m|y obs ).

(2)

The latter can be seen as a classification problem predicting the model number given the observation of y.

From this standpoint, m MAP is the Bayes classifier, well known to minimize the 0-1 loss [START_REF] Devroye | A probabilistic theory of pattern recognition, volume 31 of Applications of Mathematics[END_REF]. One might argue that m MAP is an estimator defined as the mode of the posterior probabilities which form the density of the posterior with respect to the counting measure. But the counting measure, namely δ 1 +• • •+δ M , is a canonical reference measure, since it is invariant to any permutation of {1, . . . , M } whereas no such canonical reference measure (invariant to one-to-one transformation) exists on compact subset of the real line. Thus (2) does not suffer from the drawbacks of posterior mode estimators [START_REF] Druilhet | Invariant HPD credible sets and MAP estimators[END_REF].

To approximate m MAP , ABC starts by simulating numerous triplets (m, θ m , y) from the joint Bayesian model. Afterwards, the algorithm mimics the Bayes classifier (2): it approximates the posterior probabilities by the frequency of each model number associated with simulated y's in a neighborhood of y obs . If required, we can eventually predict the best model with the most frequent model in the neighborhood, or, in other words, take the final decision by plugging in (2) the approximations of the posterior probabilities.

If directly applied, this first, naive algorithm faces the curse of dimensionality, as simulated datasets y can be complex objects and lie in a space of high dimension (e.g., numerical images). Indeed, finding a simulated dataset in the vicinity of y obs is almost impossible when the ambient dimension is high. The ABC algorithm performs therefore a (non linear) projection of the observed and simulated datasets onto some Euclidean space of reasonable dimension via a function S, composed of summary statistics. Moreover, due to obvious reasons regarding computer memory, instead of keeping track of the whole simulated datasets, one commonly saves only the simulated vectors of summary statistics, which leads to a table composed of iid replicates (m, θ m , S(y)), often called the reference save (m j , θ j , S(y j )) ← (m, θ, S(y)); end return the table of (m j , θ j , S(y j )), j = 1, . . . , n REF From the standpoint of machine learning, the reference table serves as a training database composed of iid replicates drawn from the distribution of interest, namely the joint Bayesian model. The regression problem of estimating the posterior probabilities or the classification problem of predicting a model number are both solved with nonparametric methods. The neighborhood of y obs is thus defined as simulations whose distances to the observation measured in terms of summary statistics, i.e., ρ(S(y), S(y obs )), fall below a threshold ε commonly named the tolerance level. The calibration of ε is delicate, but had been partly neglected in the papers dealing with ABC that first focused on decreasing the total number of simulations via the recourse to Markov chain Monte Carlo [START_REF] Marjoram | Markov chain Monte Carlo without likelihoods[END_REF] or sequential Monte Carlo methods [START_REF] Beaumont | Adaptive approximate Bayesian computation[END_REF][START_REF] Del Moral | An adaptive sequential monte carlo method for approximate bayesian computation[END_REF] whose common target is the joint Bayesian distribution conditioned by ρ(S(y), S(y obs )) ≤ ε for a given ε. By contrast, the simple setting we adopt here reveals the calibration question. In accepting the machine learning viewpoint, we can consider the ABC algorithm as a knearest neighbor (knn) method, see [START_REF] Biau | New insights into Approximate Bayesian Computation[END_REF]; the calibration of ε is thus transformed into the calibration of k. Unfortunately, when the summary statistics are not sufficient for the model choice problem, [START_REF] Didelot | Likelihood-free estimation of model evidence[END_REF] and [START_REF] Robert | Lack of confidence in approximate Bayesian computation model choice[END_REF] found that the above probability can greatly differ from the genuine π(m|y obs ). Afterwards [START_REF] Marin | Relevant statistics for Bayesian model choice[END_REF] provide necessary and sufficient conditions on S(•) for the consistency of the MAP based on π(m|S(y obs )) when the information included in the dataset y obs increases, i.e. when the dimension of y obs tends to infinity. Consequently, the problem that ABC addresses with reliability is classification, and the mentioned theoretical results requires a shift from the approximation of posterior probabilities. Practically the frequencies returned by Algorithm 2 should solely be used to order the models with respect to their fits to y obs and construct a knn classifier m that predicts the model number.

It becomes therefore obvious that the calibration of k should be done by minimizing the misclassification error rate of the resulting classifier m. This indicator is the expected value of the 0-1 loss function, namely 1{ m(y) = m}, over a random (m, y) distributed according to the marginal (integrated in θ m ) of the joint Bayesian distribution, whose density in (m, y) writes

π(m) f m (y|θ m )π m (θ m )dθ m . (3)
Ingenious solutions have been already proposed and are now well established to fullfil this minimization goal and bypass the overfitting problem, based on cross-validation on the learning database. But, for the sake of clarity, particularly in the following sections, we decided to take advantage of the fact that ABC aims at learning on simulated databases, and we will use a validation reference table, simulated also with Algorithm 1, but independently of the training reference table, to evaluate the misclassification rate with the averaged number of differences between the true model numbers m j and the predicted values m(y j ) by knn (i.e. by ABC) on the validation reference table.

Local error rates

The misclassification rate τ of the knn classifier m at the core of Algorithm 2 provides consistant evidence of its global accuracy. It supplies indeed a well-known support to calibrate k in Algorithm 2. The purpose of ABC model choice methods though is the analyse of an observed dataset y obs and this first indicator is irrelevant to assess the accuracy of the classifier at this precise point of the data space, since it is by nature a prior gauge. We propose here to disintegrate this indicator, and to rely on conditional expected value of the misclassification loss 1{ m(y) = m} knowing y as an evaluation of the efficiency of the classifier at y. We recall the following proposition whose proof is easy, but might help clarifying matters when applied to the joint distribution (3).

Proposition 1. Consider a classifier m that aims at predicting m given y on data drawn from the joint distribution f (m, y). Let τ be the misclassification rate of m, defined by P( m(Y ) = M ), where (M , Y ) is a random pair with distribution f under the probability measure P. Then, (i) the expectation of the loss function is

τ = m y 1{ m(y) = m}f (m, y) dy.
Additionally, (ii), the conditional expectation knowing y, namely τ (y

) = P m(Y ) = M Y = y , is τ (y) = m 1{ m(y) = m}f (m|y) (4)
and τ = y f (y)τ (y) dy, where f (y) denotes the marginal distribution of f (integrated over m) and f (m|y) = f (m, y)/f (y) the conditional probability of m given y. Furthermore, we have

τ (y) = 1 -f ( m(y) | y). (5)
The last result (5) suggests that a conditional expected value of the misclassification loss is a valuable indicator of the error at y since it is admitted that the posterior probability of the predicted model reveals the accuracy of the decision at y. Nevertheless, the whole simulated datasets are not saved into the ABC reference table but solely some numerical summaries S(y) per simulated dataset y, as explained above. Thus the disintegration process of τ is practically limited to the conditional expectation of the loss knowing some non one-to-one function of y. Its definition becomes thereforemuch more subtle than the basic (4). Actually, the ABC classifier can be trained on a subset S 1 (y) of the summaries S(y) saved in the training reference table, or on some deterministic function (we still write S 1 (y)) of S(y) that reduces the dimension, such as the projection on the LDA axes proposed by [START_REF] Estoup | Estimation of demo-genetic model probabilities with Approximate Bayesian Computation using linear discriminant analysis on summary statistics[END_REF]. To highlight this fact, the ABC classifier is denoted by m(S 1 (y)) in what follows. It is worth noting here that the above setting encompasses any dimension reduction technique presented in the review of [START_REF] Blum | A Comparative Review of Dimension Reduction Methods in Approximate Bayesian Computation[END_REF], though the review in oriented on parameter inference. Furthermore we might want to disintegrate the misclassification rate with respect to another projection S 2 (y) of the simulated data that can or cannot be related to the summaries S 1 (y) used to train the ABC classifier, albeit S 2 (y) is also limited to be a deterministic function of S(y). This yields the following.

Definition 2. The local error rate of the m(S 1 (y)) classifier with respect to S 2 (y) is

τ S1 (S 2 (y)) := P m(S 1 (Y )) = M S 2 (Y ) = S 2 (y) ,
where (M , Y ) is a random variable with distribution given in (3).

The purpose of the local misclassification rate in the present paper is twofold and requires to play with the distinction between S 1 and S 2 , as the last part will show on numerical examples. The first goal is the construction of a prospective tool that aims at checking whether a new statistic S ′ (y) carries additional information regarding the model choice, beyond a first set of statistics S 1 (y). In the latter case, it can be useful to localize the misclassification error of m(S 1 (y)) with respect to the concatenated vector S 2 (y) = (S 1 (y), S ′ (y)). Indeed, this local error rate can reveal concentrated areas of the data space, characterized in terms of S 2 (y), in which the local error rate rises above (M -1)/M , the averaged (local) amount of errors of the random classifier among M models, so as to approach 1. The interpretation of the phenomenon is as follows: errors committed by m(S 1 (y)), that are mostly spread on the S 1 (y)-space, might gather in particular areas of subspaces of the support of S 2 (y) = (S 1 (y), S ′ (y)). This peculiarity is due to the dimension reduction of the summary statistics in ABC before the training of the classifier and represents a concrete proof of the difficulty of ABC model choice already raised by [START_REF] Didelot | Likelihood-free estimation of model evidence[END_REF] and [START_REF] Robert | Lack of confidence in approximate Bayesian computation model choice[END_REF].

The second goal of the local error rate given in Definition 2 is the evaluation of the confidence we may concede in the model predicted at y obs by m(S 1 (y)), in which case we set S 2 (y) = S 1 (y). And, when both sets of summaries agree, the results of Proposition 1 extend to

τ S1 (S 1 (y)) = m π(m|S 1 (y))1{ m(S 1 (y)) = m} = 1 -π( m(S 1 (y)) | S 1 (y)). ( 6
)
Besides the local error rate we propose in Definition 2 is an upper bound of the Bayes classifier if we admit the loss of information committed by replacing y with the summaries.

Proposition 3. Consider any classifier m(S 1 (y)).

The local error rate of this classifier satisfies

τ S1 (s 2 ) = P ( m(S 1 (Y )) = M | S 2 (Y ) = s 2 ) ≥ P ( m MAP (Y ) = M | S 2 (Y ) = s 2 ) , (7)
where m MAP is the Bayes classifier defined in (2) and s 2 any value in the support of S 2 (Y ). Consequently,

P ( m(S 1 (Y )) = M ) ≥ P ( m MAP (Y ) = M ) . ( 8
)
Proof. Proposition 1, in particular ( 5), implies that m MAP (y) is the ideal classifier that minimizes the conditional 0-1 loss knowing y. Hence, we have

P ( m(S 1 (Y )) = M | Y = y) ≥ P ( m MAP (Y ) = M | Y = y) .
Integrating the above with respect to the distribution of Y knowing S 2 (Y ) leads to (7), and a last integral to (8).

Proposition 3 shows that the introduction of new summary statistics cannot distort the model selection insofar as the risk of the resulting classifier cannot decrease below the risk of the Bayes classifier m MAP . We give here a last flavor of the results of [START_REF] Marin | Relevant statistics for Bayesian model choice[END_REF] and mention that, if S 1 (y) = S 2 (y) = S(y) and if the classifiers are perfect (i.e., trained on infinite reference tables), we can rephrase part of their results as providing mild conditions on S under which the local error τ S (S(y)) tends to 0 when the size of the dataset y tends to infinity.

Estimation algorithm of the local error rates

The numerical estimation of the local error rate τ S1 (S 2 (y)), as a surface depending on S 2 (y), is therefore paramount to assess the difficulty of the classification problem at any s 2 = S S (y), and the local accuracy of the classifier. Naturally, when S 1 (y) = S 2 (y) for all y, the local error can be evaluated at S 2 (y obs ) by plugging in (6) the ABC estimates of the posterior probabilities (the relative frequencies of each model among the particles returned by Algorithm 2) as substitute for π(m|S(y obs )). This estimation procedure is restricted to the above mentioned case where the set of statistics used to localize the error rate agrees with the set of statistics used to train the classifier.

Moreover, the approximation of the posterior probabilities returned by Algorithm 2, i.e., a knn method, might not be trustworthy: the calibration of k performed by minimizing the prior error rate τ does not provide any certainty on the estimated posterior probabilities beyond a ranking of these probabilities that yields the best classifier in terms of misclassification. In other words, the knn method calibrated to answer the classification problem of discriminating among models does not produce a reliable answer to the regression problem of estimating posterior probabilities. Certainly, the value of k must be increased to face this second kind of issue, at the price of a larger bias that might even swap the model ranking (otherwise, the empirical prior error rate would not depend on k, see the numerical result section). For all these reasons, we propose here an alternative estimate of the local error. The core idea of our proposal is the recourse to a nonparametric method to estimate conditional expected values based on the calls to the classifier m on a validation reference table, already simulated to estimate the global error rate τ . Nadaraya-Watson kernel estimators of the conditional expected values

τ S1 (S 2 (y)) = E (1{ m(S 1 (Y )) = M } | S 2 (Y ) = S 2 (y)) (9)
rely explicitly on the regularity of this indicator, as a function of s 2 = S 2 (y), which contrasts with the ABC plug-in estimate described above. We thus hope improvements in the accuracy of error estimate and a more reliable approximation of the whole function τ S1 (S 2 (y)). Additionally, we are not limited anymore to the special case where S 1 (y) = S 2 (y) for all y. It is worth stressing here that the bandwidth of the kernels must be calibrated by minimizing the L 2 -loss, since the target is a conditional expected value.

Practically, this leads to Algorithm 3 which requires a validation or test reference table independent of the training database that constitutes the ABC reference table. We can bypass the requirement by resorting to cross validation methods, as for the computation of the global prior misclassification rate τ . But the ensued algorithm is complex and it induces more calls to the classifier (consider, e.g,. a ten-fold cross validation algorithm computed on more than one random grouping of the reference table) than the basic Algorithm 3, whereas the training database can always be supplemented by a validation database since ABC, by its very nature, is a learning problem on simulated databases. Moreover, to display the whole We performed numerical experiments (not detailed here) concluding that the resort to a Kriging algorithm provides results comparable to the evaluation of Nadaraya-Watson estimator at any point of the support of S 2 (y), and can reduce computation times.

Adaptive ABC

The local error rate can also represent a valuable way to adjust the summary statistics to the data point y and to build an adaptive ABC algorithm achieving a local trade off that increases the dimension of the summary statistics at y only when the additional coordinates add information regarding the classification problem. Assume that we have at our disposal a collection of ABC classifiers, m λ (y) := m λ (S λ (y)), λ = 1, . . . , Λ, trained on various projections of y, namely the S λ (y)'s, and that all these vectors, sorted with respect to their dimension, depend only on the summary statistics registered in the reference tables. Sometimes low dimensional statistics may suffice for the classification (of models) at y, whereas other times we may need to examine statistics of larger dimension. The local adaptation of the classifier is accomplished through the disintegration of the misclassification rates of the initial classifiers with respect to a common statistic S 0 (y). Denoting τ λ (S 0 (y)) the lo-cal error rate of m λ (y) knowing S 0 (y), this reasoning yields the adaptive classifier defined by m(S(y)) := m λ(y) (y), where λ(y) := arg min λ=1,...,Λ τ λ (S 0 (y)). (10) This last classifier attempts to avoid bearing the cost of the potential curse of dimensionality from which all knn classifiers suffer and can help reduce the error of the initial classifiers, although the error of the ideal classifier (2) remains an absolute lower bound, see Proposition 3. From a different perspective, (10) represents a way to tune the similarity ρ(S(y), S(y obs )) of Algorithm 2 that locally includes or excludes components of S(y) to assess the proximity between S(y) and S(y obs ). Practically, we rely on the following algorithm to produce the adaptive classifier, that requires a validation reference table independent of the reference table used to fit the initial classifiers. The local error surface estimated within the loop of Algorithm 4 must contrast the errors of the collection of classifiers. Our advice is thus to build a projection S 0 (y) of the summaries S(y) registered in the reference tables as follow. Add to the validation reference table a qualitative trait which groups the replicates of the table according to their differences between the predicted numbers by the initial classifiers and the model numbers m j registered in the database. For instance, when the collection is composed of Λ = 2 classifiers, the qualitative trait takes three values: value 0 when both classifiers m λ (y j ) agree (whatever the value of m j ), value 1 when the first classifier only returns the correct number, i.e., m 1 (y j ) = m j = m 2 (y j ), and value 2 when the second classifier only returns the correct number, i.e., m 1 (y j ) = m j = m 2 (y j ). The axes of the linear discriminant analysis (LDA) predicting the qualitative trait knowing S(y) provide a projection S 0 (y) which contrasts the errors of the initial collection of classifiers.

Finally it is important to note that the local error rates are evaluated in Algorithm 4 with the help of a validation reference table. Therefore, a reliable estimation of the accuracy of the adaptive classifier cannot be based on the same validation database because of the optimism bias of the training error. Evaluating the accuracy requires the simulation of a test reference table independently of the two first databases used to train and adapt the predictor, as is usually performed in the machine learning community.

Hidden random fields

Our primary intent with the ABC methodology exposed in Section 2 was the study of new summary statistics to discriminate between hidden random fields models. The following materials numerically illustrate how ABC can choose the dependency structure of latent Potts models among two possible neighborhood systems, both described with undirected graphs, whilst highlighting the generality of the approach.

Hidden Potts model

This numerical part of the paper focuses on hidden Potts models, that are representative of the general level of difficulty while at the same time being widely used in practice (see for example [START_REF] Hurn | A Tutorial on Image Analysis[END_REF][START_REF] Alfò | A finite mixture model for image segmentation[END_REF][START_REF] François | Bayesian Clustering Using Hidden Markov Random Fields in Spatial Population Genetics[END_REF], Moores et al., 2014). we recall that the latent random field x is a family of random variables x i indexed by a finite set S , whose elements are called sites, and taking values in a finite state space X := {0, . . . , K -1}, interpreted as colors. When modeling a digital image, the sites are lying on a regular 2D-grid of pixels, and their dependency is given by an undirected graph G which defines an adjacency relationship on the set of sites S : by definition, both sites i and j are adjacent if and only if the graph G includes an edge that links directly i and j. A Potts model sets a probability distribution on x, parametrized by a scalar β that adjusts the level of dependency between adjacent sites. The latter class of models differs from the auto-models of [START_REF] Besag | Spatial interaction and the statistical analysis of lattice systems (with Discussion)[END_REF], that allow variations on the level of dependencies between edges and introduce potential anisotropy on the graph. But the difficulty of all these models arises from the intractable normalizing

π(x|G , β) = 1 Z(G , β) exp   β i G ∼j 1{x i = x j }    .
The above sum i G ∼ j ranges the set of edges of the graph G and the normalizing constant Z(G , β) writes as

Z(G , β) = x∈X exp   β i G ∼j 1{x i = x j }    , (11)
namely a summation over the numerous possible realizations of the random field x, that cannot be computed directly (except for small grids and small number of colors K). In the statistical physic literature, β is interpreted as the inverse of a temperature, and when the temperature drops below a fixed threshold, values x i of a typical realization of the field are almost all equal (due to important dependency between all sites). These peculiarities of Potts models are called phase transitions.

In hidden Markov random fields, the latent process is observed indirectly through another field; this permits the modeling of a noise that may be encountered in many concrete situations. Precisely, given the realization x of the latent field, the observation y is a family of random variables indexed by the set of sites, and taking values in a set Y , i.e., y = (y i ; i ∈ S ), and are commonly assumed as independent draws that form a noisy version of the hidden fields. Consequently, we set the conditional distribution of y knowing x as the product π(y|x, α) = i∈S P (y i |x i , α), where P is the marginal noise distribution parametrized by some scalar α. Hence the likelihood of the hidden Potts model with parameter β on the graph G and noise distribution P α , denoted HPM(G , α, β), is given by

f (y|α, β, G ) = x∈X π (x|G , β) π α (y|x)
and faces a double intractable issue as neither the likelihood of the latent field, nor the above sum can be computed directly: the cardinality of the range of the sum is of combinatorial complexity. The following numerical experiments are based on two classes of noises, producing either observations in {0, 1, . . . , K-1}, the set of latent colors, or continuous observations that take values in R.

The common point of our examples is to select the hidden Gibbs model that better fits a given y obs composed of N = 100 × 100 pixels within different neighborhood systems represented as undirected graphs G . We considered the two widely used adjacency structures in our simulations, namely the graph G 4 (respectively G 8 ) in which the neighborhood of a site is composed of the four (respectively eight) closest sites on the two-dimensional lattice, except on the boundaries of the lattice, see Fig. 1. The prior probabilities of both models were set to 1/2 in all experiments. The Bayesian analysis of the model choice question adds another integral beyond the two above mentioned sums that cannot be calculated explicitly or numerically either and the problem we illustrate are said triple intractable. Up to our knowledge the choice of the latent neighborhood structure has never been seriously tackled in the Bayesian literature. We mentioned here the mean field approximation of [START_REF] Forbes | Hidden Markov random field model selection criteria based on mean fieldlike approximations[END_REF] whose software can estimate paramaters of such models, and compare models fitness via a BIC criterion. But the first results we try to obtain with this tool were worse than with ABC, except for very low values of β. This means that either we did not manage to run the software properly, or that the mean field approximation is not appropriate to discriminate between neighborhood structures. The detailed settings of our three experiments are as follows.

First experiment. We considered Potts models with K = 2 colors and a noise process that switches each pixel independently with probability exp(-α)/(exp(α) + exp(-α)), following the proposal of [START_REF] Everitt | Bayesian Parameter Estimation for Latent Markov Random Fields and Social Networks[END_REF]. The prior on α was uniform over (0.42; 2.3), where the bounds of the interval were determined to switch a pixel with a probability less than 30%. Regarding the dependency parameter β, we set prior distributions below the phase transition which occurs at different levels depending on the neighborhood structure. Precisely we used a uniform distribution over (0; 1) when the adjacency is given by G 4 and a uniform distribution over (0; 0.35) with G 8 .

Second experiment. We increased the number of colors in the Potts models and set K = 16. Likewise, we set a noise that changes the color of each pixel with a given probability parametrized by α, and conditionally on a change at site i, we rely on the least favorable distribution, which is a uniform draw within all colors except the latent one. To extend the parametrization of [START_REF] Everitt | Bayesian Parameter Estimation for Latent Markov Random Fields and Social Networks[END_REF], the marginal distribution of the noise is defined by

P α (y i |x i ) = exp α 21{x i = y i } -1 exp(α) + (K -1) exp(-α)
and a uniform prior on α over the interval (1.78; 4.8) ensures that the probability of changing a pixel with the noise process is at most 30%. The uniform prior on the Potts parameter β was also tuned to stay below the phase transition. Hence β ranges the interval (0; 2.4) with a G 4 structure and the interval (0; 1) with a G 8 structure.

Third experiment. We introduced a homoscedastic Gaussian noise whose marginal distribution is characterized by

y i | x i = c ∼ N (c, σ 2 ) c ∈ {0; 1}
over bicolor Potts models. And both prior distributions on parameter β are similar to the ones on the latent fields of the first experiment. The standard deviation σ = 0.39 was set so that the probability of a wrong prediction of the latent color with a marginal MAP rule on the Gaussian model is about 15%.

Geometric summary statistics

Performing a Bayesian model choice via ABC algorithms requires summary statistics that capture the relevant information from the observation y obs to discriminate among the competing models. When the observation is noise-free, [START_REF] Grelaud | ABC likelihood-free methods for model choice in Gibbs random fields[END_REF] noted that the joint distribution resulting from the Bayesian modeling falls into the exponential family, and they obtained consecutively a small set of summary statistics, depending on the collection of considered models, that were sufficient. In front of noise, the situation differs substantially as the joint distribution lies now outside the exponential family, and the above mentioned statistics are not sufficient anymore, whence the urge to bring forward other concrete and workable statistics. The general approach we developed reveals geometric features of a discrete field y via the recourse to colored graphs attached to y and their connected components. Consider an undirected graph G whose set of vertices coincides with S , the set of sites of y.

Definition 4. The graph induced by G on the field y, denoted Γ(G , y), is the undirected graph whose set of edges gathers the edges of G between sites of y that share the same color, i.e., i Γ(G ,y)

∼ j ⇐⇒ i G ∼ j and y i = y j .

We believe that the connected components of such induced graphs capture major parts of the geometry of y. Recall that a connected component of an undirected graph Γ is a subgraph of Γ in which any two vertices are connected to each other by a path, and which is connected to no other vertices of Γ. And the connected components form a partition of the vertices. Since ABC relies on the computation of the summary statistics on many simulated datasets, it is also worth noting that the connected components can be computed efficiently with the help of famous graph algorithms in linear time based on a breadthfirst search or depth-first search over the graph. The empirical distribution of the sizes of the connected components represents an important source of geometric informations, but cannot be used as a statistic in ABC because of the curse of dimensionality. The definition of a low dimensional summary statistic derived from these connect components should be guided by the intuition on the model choice we face.

Our numerical experiments discriminate between a G 4 -and a G 8 -neighborhood structure and we considered two induced graphs on each simulated y, namely Γ(G 4 , y) and Γ(G 8 , y). Remark that the twodimensional statistics proposed by [START_REF] Grelaud | ABC likelihood-free methods for model choice in Gibbs random fields[END_REF], which are sufficient in the noise-free context, are the total numbers of edges in both induced graphs. After very few trials without success, we fixed ourselves on four additional summary statistics, namely the size of the largest component of each induced graph, as well as the total number of connect components in each graph. See Fig. 2 for an example on a bicolor picture y. To fix the notations, for any induced graph Γ(G , y), we define

Γ(G 4 , y) Γ(G 8 , y)
• R(G , y) as the total number of edges in Γ(G , y),

• T (G , y) as the number of connected components in Γ(G , y) and

• U (G , y) as the size of the largest connected component of Γ(G , y).

And to sum up the above, the set of summary statistics that where registered in the reference tables for each simulated field y is

S(y) = R(G 4 , y); R(G 8 , y); T (G 4 , y); T (G 8 , y); U (G 4 , y); U (G 8 , y)
in the first and second experiments.

In the third experiment, the observed field y takes values in R and we cannot apply directly the approach based on induced graphes because no two pixels share the same color. All of the above statistics are meaningless, including the statistics R(G , y) used by [START_REF] Grelaud | ABC likelihood-free methods for model choice in Gibbs random fields[END_REF] in the noise-free case. We rely on a quantization preprocessing performed via a kmeans algorithm on the observed colors that forgets the spatial structure of the field. The algorithm was tuned to uncover the same number of groups of colors as the number of latent colors, namely K = 2. If q 2 (y) denotes the resulting field, the set of summary statistics becomes S(y) = R G 4 , q 2 (y) ; R G 8 , q 2 (y) ; T G 4 , q 2 (y) ; T G 8 , q 2 (y) ; U G 4 , q 2 (y) ; U G 8 , q 2 (y) .

We have assumed here that the number of latent colors is known to keep the same purpose of selecting the correct neighborhood structure. Indeed [START_REF] Cucala | Bayesian Inference on a Mixture Model With Spatial Dependence[END_REF] have already proposed a (complex) Bayesian method to infer the appropriate number of hidden colors. But more generally, we can add statistics based on various quantizations q k (y) of y with k groups.

Numerical results

In all three experiments, we compare three nested sets of summary statistics S 2D (y), S 4D (y) and S 6D (y) of dimension 2, 4 and 6 respectively. They are defined as the projection onto the first two (respectively four and six) axes of S(y) described in the previous section. We stress here that S 2D (y), which is composed of the summaries given by [START_REF] Grelaud | ABC likelihood-free methods for model choice in Gibbs random fields[END_REF], are used beyond the noise-free setting where they are sufficient for model choice. In order to study the information carried by the connected components, we add progressively our geometric summary statistics to the first set, beginning by the T (G , y)-type of statistics in S 4D (y). Finally, remark that, before evaluating the Euclidean distance in ABC algorithms, we normalize the statistics in each reference tables with respect to an estimation of their standard deviation since all these summaries take values on axis of different scales. Simulated images have been drawn thanks to the [START_REF] Swendsen | Nonuniversal critical dynamics in Monte Carlo simulations[END_REF] algorithm. In the least favorable experiment, simulations of one hundred pictures (on pixel grid of size 100 × 100) via 20, 000 iterations of this Markovian algorithm when parameters drawn from our prior requires about one hour of computation on a single CPU with our optimized C++ code. Hence the amount of time required by ABC is dominated by the simulations of y via the Swedsen-Wang algorithm. This motivated [START_REF] Moores | Preprocessing for approximate Bayesian computation in image analysis[END_REF] to propose a cut down on the cost of running an ABC experiment by removing the simulation of an image from hidden Potts model, and replacing it by an approximate simulation of the summary statistics. Another alternative is the clever sampler of [START_REF] Mira | Perfect slice samplers[END_REF] that provides exact simulations of Ising models and can be extended to Potts models.

First experiment. Fig. 3(a) illustrates the calibration of the number of nearest neighbors (parameter k of Algorithm 2) by showing the evolution of the prior error rates (evaluated on a validation reference table including 20, 000 simulations) when k increases. We compared the errors of six classifiers to inspect the differences between the three sets of summary statistics (in yellow, green and magenta) and the impact of the size of the training reference table (100, 000 simulations in solid lines; 50, 000 simulations in dashed lines). The numerical results exhibit that a good calibration of k can reduce the prior misclassification error. Thus, without really degrading the performance of the classifiers, we can reduce the amount of simulations required in the training reference table, whose computation cost (in time) represents the main obstacle of ABC methods, see also Table 1. Moreover, as can be guessed from Fig. 3(a), the sizes of the largest connected components of induced graphs (in- One can argue that the curse of dimensionality does not occur with such low dimensional statistics and sizes of the training set, but this intuition is wrong, as shown in Fig. 3(b). The latter plot shows indeed the prior misclassification rate as a function of k when we replace the last four summaries by ancillary statistics drawn independently of m and y. We can conclude that, although the three sets of summary statistics carry then the same information in this artificial setting, the prior error rates increase substantially with the dimension (classifiers are not trained on infinite reference tables!). This conclusion shed new light on the results of Fig. 3(a): the U (G , y)-type summaries, based on the size of the largest component, are not concretely able to help discriminate among models, but are either highly correlated with the first four statistics; or the resolution (in terms of size of the training reference table) does not permit the exploitation of the possible information they add.

Fig. 3(c) displays the local error rate with respect to a projection of the image space on a plan. We have taken here S 1 (y) = S 2D (y) in Definition 2. And S 2 (y) ranges a plan given by a projection of the full set of summaries that has been tuned empirically in order to gather the errors committed by calls of m(S 2D (y)) on the validation reference table. The most striking fact is that the local error rises above 0.9 in the oval, reddish area of Fig. 3(c). Other reddish areas of Fig. 3(c) in the bottom of the plot correspond to parts of the space with very low probability, and may be a dubious extrapolation of the Kriging algorithm. We can thus conclude that the information of the new geometric summaries depends highly on the position of y in the image space and have confidence in the interest of Algorithm 4 (adaptive ABC) in this frame-work. As exhibited in Table 1(d), this last classifier does not decrease dramatically the prior misclassification rates. But the errors of the non-adaptive classifiers are already low and the error of any classifier is bounded from below, as explained in Proposition 3. Interestingly though, the adaptive classifier relies on m(S 2D (y)) (instead of the most informative m(S 6D (y))) to take the final decision at about 60% of the images of our test reference table of size 30, 000.

Second experiment. The framework was designed here to study the limitations of our approach based on the connected components of induced graphs. The number of latent colors is indeed relatively high and the noise process do not rely on any ordering of the colors to perturbate the pixels. Table 2 indicates the difficulty of capturing relevant information with the geometric summaries we propose. Only the sharpness introduced by a training reference table composed of 100, 000 simulations distinguishes m(S 4D (y)) and m(S 6D (y)) from the basic classifier m(S 2D (y)). This conclusion is reinforced by the low value of number of neighbors after the calibration process, namely k = 16, 5 and 5 for m(S 2D (y)), m(S 4D (y)) and m(S 6D (y)) respectively. Hence we do not display in the paper other diagnosis plots based on the prior error rates or the conditional error rates, which led us to the same conclusion. The adaptive ABC algorithm did not improve any of these results.

Third experiment. The framework here includes a continuous noise process as described at the end of Section 3.1. We reproduced the entire diagnosis process performed in the first experiment and we obtained the results given in Fig. 4 and Table 3. The most noticeable difference is the extra information carried by the U (G , y)-statistics, representing the size of the largest connected component, and the adaptive ABC relie on the simplest m(S 2D (y)) in about 30% of the data space (measured with the prior marginal

Conclusion and perspective

In the present article, we considered ABC model choice as a classification problem in the framework of the Bayesian paradigm (Section 2.1) and provided a local error in order to assess the accuracy of the classifier at y obs (Sections 2.2 and 2.3). We derived then an adaptive classifier (Section 2.4) which is an attempt to fight against the curse of dimensionality locally around y obs . This method contrasts with most projection methods which are focused on parameter estimation [START_REF] Blum | A Comparative Review of Dimension Reduction Methods in Approximate Bayesian Computation[END_REF]. Additionally, most of them perform a global trade off between the dimension and the information of the summary statistics over the whole prior domain, while our proposal adapts the dimension with respect to y obs (see also the discussion about the posterior loss approach in [START_REF] Blum | A Comparative Review of Dimension Reduction Methods in Approximate Bayesian Computation[END_REF]. Besides the inequalities of Proposition 3 complement modestly the analysis of [START_REF] Marin | Relevant statistics for Bayesian model choice[END_REF] on ABC model choice. Principles of our proposal are well founded by avoiding the well-known optimism of the training error rates and by resorting to validation and test reference tables in order to evaluate the error practically. And, finally, the machine learning viewpoint gives an efficient way to calibrate the threshold of ABC (Section 2.1).

Regarding latent Markov random fields, the proposed method of constructing summary statistics based on the induced graphs (Section 3.2) yields a promising route to construct relevant summary statistics in this framework. This approach is very intuitive and can be reproduced in other settings. For instance, if the goal of the Bayesian analysis is to select between isotropic latent Gibbs models and anisotropic models, the averaged ratio between the width and the length of the connect components or the ratio of the width and the length of the largest connected components can be relevant numerical summaries. We have also explained how to adapt the method to a continuous noise by performing a quantization of the observed values at each site of the fields (Section 3.2). And the detailed analysis of the numerical results demonstrates that the approach is promising. However the results on the 16 color example with a completely disordered noise indicate the limitation of the induced graph approach. We believe that there exists a road we did not explore above with an induced graph that add weights on the edges of the The numerical results (Section 3.3) highlighted that the calibration of the number of neighbors in ABC provides better results (in terms of misclassification) than a threshold set as a fixed quantile of the distances between the simulated and the observed datasets (as proposed in [START_REF] Marin | Approximate Bayesian Computational methods[END_REF]. Consequently, we can reduce significantly the number of simulations in the reference table without increasing the misclassification error rates. This represents an important conclusion since the simulation of a latent Markov random field requires a non-negligible amount of time. The gain in misclassification rates of the new summaries is real but not spectacular and the adaptive ABC algorithm was able to select the most performant classifier.

  Assume we are given M stochastic models with respective paramater spaces embedded into Euclidean spaces of various dimensions. The joint Bayesian distribution sets (i) a prior on the model space, π(1), . . . , π(M ), (ii) for each model, a prior on its parameter space, whose density with respect to a reference measure (often the Lebesgue measure of the Euclidean space) is π m (θ m ) and (iii) the likelihood of the data y within each model, namely f m (y|θ m ). The evidence of model m is then defined as e(m, y) := f m (y|θ m )π m (θ m )dθ m and the posterior probability of model m as π(m|y) = π(m)e(m, y) m ′ π(m ′ )e(m ′ , y)

Algorithm 4 :

 4 Adaptive ABC model choice Input: A collection of classifiers m λ (y), λ = 1, . . . , Λ and a validation reference table Output: An adaptive classifier m(y) for each λ ∈ {1, . . . , Λ} do estimate the local error of m λ (y) knowing S 0 (y) with the help of Algorithm 3; end return the adaptive classifier m as a function computing (10);

Figure 1 :

 1 Figure 1: Neighborhood graphs G of hidden Potts model. (a) The four closest neighbour graph G 4 defining model HPM(G 4 , α, β). (b) The eight closest neighbour graph G 8 defining model HPM(G 8 , α, β)

Figure 2 :

 2 Figure 2: The induced graph Γ(G 4 , y) and Γ(G 8 , y) on a given bicolor image y of size 5 × 5. The six summary statistics on y are thus R(G 4 , y) = 22, T (G 4 , y) = 7, U (G 4 , y) = 12, R(G 8 , y) = 39, T (G 8 , y) = 4 and U (G 8 , y) = 16

Figure 3 :

 3 Figure 3: First experiment results. (a) Prior error rates (vertical axis) of ABC with respect to the number of nearest neighbors (horizontal axis) trained on a reference table of size 100, 000 (solid lines) or 50, 000 (dashed lines), based on the 2D, 4D and 6D summary statistics. (b) Prior error rates of ABC based on the 2D summary statistic compared with 4D and 6D summary statistics including additional ancillary statistics. (c) Evaluation of the local error on a 2D surface

Figure 4 :

 4 Figure 4: Third experiment results. (a) Prior error rates (vertical axis) of ABC with respect to the number of nearest neighbors (horizontal axis) trained on a reference table of size 100, 000 (solid lines) or 50, 000 (dashed lines), based on the 2D, 4D and 6D summary statistics. (b) Prior error rates of ABC based on the 2D summary statistics compared with 4D and 6D summary statistics including additional ancillary statistics. (c) Evaluation of the local error on a 2D surface

  table in the ABC literature, see Algorithm 1. Simulation of the ABC reference table Output: A reference table of size n REF

	Algorithm 1:

for j ← 1 to n REF do draw m from the prior π; draw θ from the prior π m ; draw y from the likelihood f m (•|θ); compute S(y);

  The Algorithm we have to calibrate is given in Algorithm 2.

	Algorithm 2: Uncalibrated ABC model choice
	Output: A sample of size k distributed
	according to the ABC approximation
	of the posterior
	simulate the reference table T according to
	Algorithm 1;
	sort the replicates of T according to
	ρ(S(y j ), S(y obs ));
	keep the k first replicates;
	return the relative frequencies of each model
	among the k first replicates and the most
	frequent model;
	Before entering into the tuning of k, we highlight
	that the projection via the summary statistics gen-
	erates a difference with the standard knn methods.
	Under mild conditions, knn are consistent nonpara-
	metric methods. Consequently, as the size of the ref-
	erence table tends to infinity, the relative frequency
	of model m returned by Algorithm 2 converges to
	π(m|S(y obs )).

table Output :

 Output Algorithm 3: Estimation of τ S1 (S 2 (y)) given an classifier m(S 1 (y)) on a validation or test reference table Input: A validation or test reference table and a classifier m(S 1 (y)) fitted with a first reference Estimations of (9) at each point of the second reference table for each (m j , y j ) in the test table do compute δ j = { m(S 1 (y j )) = m j }; end calibrate the bandwidth h of the Nadaraya-Watson estimator predicting δ j knowing S 2 (y j ) via cross-validation on the test table; for each (m j , y j ) in the test table do evaluate the Nadaraya-Watson estimator with bandwidth h at S 2 (y j ); end surface τ S1 (S 2 (y)), we can interpolate values of the local error between points S 2 (y) of the second reference table with the help of a Kriging algorithm.

Table 1 :

 1 Evaluation of the prior error rate on a test reference table of size 30, 000 in the first experiment.

	Prior error rates
	Train size	5, 000 100, 000
	2D statistics	8.8%	7.9%
	4D statistics	6.5%	6.1%
	6D statistics	7.1%	7.1%
	Adaptive ABC	6.2%	5.5%

Table 2 :

 2 Evaluation of the prior error rate on a test reference table of size 20, 000 in the second experiment.

	Prior error rates
	Train size 50, 000 100, 000
	2D statistics	4.5%	4.4%
	4D statistics	4.6%	4.1%
	6D statistics	4.6%	4.3%
	cluded only in S 6D (y)) do not carry additional infor-
	mation regarding the model choice and Table 1 con-
	firms this results through evaluations of the errors on
	a test reference table of 30, 000 simulations drawn in-
	dependently of both training and validation reference
	tables.

Table 3 :

 3 Evaluation of the prior error rate on a test reference table of size 30, 000 in the third experiment. Likewise, the gain in misclassification errors is not spectacular, albeit positive.

	Prior error rates
	Train size	5, 000 100, 000
	2D statistics	14.2%	13.8%
	4D statistics	10.8%	9.8%
	6D statistics	8.6%	6.9%
	Adaptive ABC	8.2%	6.7%

Acknowledgments

The three author were financially supported by the Labex NUMEV. We are grateful to Jean-Michel Marin for his constant feedback and support. Some part of the present work was presented at MCMSki 4 in January 2014 and benefited greatly from discussions with the participants during the poster session. We would like to thank the anonymous referrees and the Editors whose valuable comments and insightful suggestions led to an improved version of the paper.


Motion Planning

Problem Given x, y ∈ M, find an admissible trajectory steering the system from x to y , possibly under some constraints.

Possible constraints:

1 Avoiding some obstacles 2 Rendez-vous problem, i.e., being near certain places at certain times

Assumption

A metric with balls B(q, ε) is fixed on M.

Method

Different approaches are possible. We consider the following method:

x y

x y

Γ

x y

Γ

We focus on quantifying the difficulty of the second step.

Motion planning problem

Definitions of complexity

Asymptotic estimates in affine-control systems

Complexity

Let J : U → [0, +∞) be a cost function.

Definition (Complexity)

A measure of the cost of approximation of a given curve/path with a certain precision

In general:

1 we fix a set Adm(Γ, ε) of admissible controls for precision ε 2 we define complexity as

σ(γ, ε) = inf u∈Adm(Γ,ε) cost of u cost of an "ε piece" of u = 1 ε inf u∈Adm(Γ,ε) J(u, T ).

Obstacle-avoidance problem

Let Γ ⊂ M be a curve, Tube(Γ, ε) = q∈γ B(q, ε), and

A(Γ, ε) = u ∈ L 1 ([0, T ], R m )) | T > 0, q u (T ) = y , q u (•) ⊂ Tube(Γ, ε) .
With this set we define the tubular approximation complexity

Σ a (Γ, ε) = 1 ε inf u∈A(Γ,ε) J(u, T ).
x

Γ y

Tube(Γ, ε)

q u (•)

Rendez-vous problem

Let γ : [0, T ] → M be a path and

N (γ, ε) = u ∈ L 1 ([0, T ], R m )) | q u (T ) = y and q u (t) ∈ B(γ(t), ε) for any t ∈ [0, T ] .
This set defines the neighborhood approximation complexity

σ n (γ, ε) = 1 ε inf u∈N (γ,ε) J(u, T ). x y q u (t) γ(t) B(γ(t), ε)
Tubular approximation complexity ⇐⇒ Neighborhood approximation complexity

Nonholonomic complexities

Introduced by Gromov (1996) in a different context.

Weak equivalence:

σ(Γ, ε) ≍ g (ε) ⇐⇒ C 1 ≤ σ(Γ, ε) g (ε) ≤ C 2 for ε ↓ 0.
Complete results (Jean 2003).

Strong equivalence:

σ(Γ, ε) ≃ g (ε) ⇐⇒ lim ε↓0 σ(Γ, ε) g (ε) = 1.
Results in particular cases (Gauthier, Zakalyukin, et al., 2004-2013) General case

Recall the general form of a control-affine system

q(t) = f 0 (q(t)) + m i=1 u i (t)f i (q(t)).
We will consider: strong Hörmander condition: Lie q {f 1 , . . . , f m } = T q M for any q ∈ M.

The set of controls is U =

T ∈(0,T ]

L1 ([0, T ], R m ).

The cost J is the L 1 -norm of u. Consequences:

Different approaches are possible. We consider the following method:1 Find an (non-admissible) curve Γ ⊂ M or a path γ : [0, T ] → M solving the problem.

Find an (non-admissible) curve Γ ⊂ M or a path γ : [0, T ] → M solving the problem.

Track Γ or γ with an admissible trajectory.

Find an (non-admissible) curve Γ ⊂ M or a path γ : [0, T ] → M solving the problem.→ global topology

Track Γ or γ with an admissible trajectory. → local behavior of the control system

The value function associated to this system w.r.t. the L 1 cost is a distance, called sub-Riemannian distance.

Due to the linearity of the system, we can always reparametrize trajectories without changing their L 1 cost. Hence,

Small time local controllability.

The associated driftless system (f 0 = 0) is a nonholonomic system.

Particular case: nonholonomic control systems Nonholonomic control system = control-affine system without drift q(t) = m i=1 u i (t)f i (q(t)), that satisfies the Hörmander condition, i.e., such that Lie q {f 1 , . . . , f m } = T q M, for any q ∈ M.

Complexities for control-affine systems

We will use the sub-Riemannian metric to define the complexities.

Since the system is not linear, we cannot reparametrize the trajectories, and hence

Tubular approximation complexity

✘ ✘ ✘

⇐⇒

Neighborhood approximation complexity

For any q ∈ M, s ∈ N, let

Hypothesis Equiregularity: for any s ∈ N, dim ∆ s does not depend on the point q ∈ M.

The complexity of curves is not sensible to the drift.

The complexity of paths depends on the drift. In particular, when f 0 ⊂ ∆ r \ ∆ r -1 where r is such that ∆ r = T q M, the complexity is always maximal, i.e., σ n (γ, ε) ≍ ε -r .

Techniques and Remarks

Weak estimates of the value function near a point (generalization of the sub-Riemannian Ball-Box theorem).

Example f 1 and f 2 control vector fields on R 3 satisfying the Hörmander condition, Drift s.t.

Techniques and Remarks (continued)

Estimates obtained by reducing the control system with drift to a driftless but time-dependent system.

For this system we can define a generalization of the nilpotent approximation, that yields the estimates.

Final remarks

We studied also two other notions of complexity, where we track the curve/path by interpolation, and no metric is assumed.

We studied also another cost

Thank you for your attention.

Ball-Box

Let {∂ z i } n i=1 be the canonical basis of R n and R f0 (q, ε) the reachable set from q with cost ≤ ε. We define

Theorem

Let z = (z 1 , . . . , z n ) a privileged coordinate system at q for {f 1 , . . . , f m }, rectifying f 0 as the k-th coordinate vector field ∂ z ℓ , for some 1 ≤ ℓ ≤ n. Them, there exist C , ε 0 , T 0 s.t., if T < T 0 , it holds


A second order finite-volume scheme on cartesian grids for Euler equations Y. Gorsse a , A. Iollo b , L. Weynans b a. CEA Cadarache b. Univ. Bordeaux, IMB, UMR5251, F-33400, Talence, France CNRS, IMB, UMR5251, F-33400, Talence, France INRIA, F-33400, Talence, France

Résumé :

Nous présentons un schéma de type volumes-finis d'ordre deux pour la résolution des équations d'Euler compressible dans des géométries complexes. La discrétisation est faite sur une grille cartésienne, qui en général ne coincide pas avec la géométrie considérée. Dans le domaine fluide, nous utilisons une méthode de volumes-finis classique. Sur les cellules situées près de la frontière avec le solide, nous résolvons un problème de Riemann ad hoc qui tient compte des conditions aux limites à imposer. Pour éviter les oscillations de pression près du solide, nous effectuons une pondération de la condition aux limites avec une extrapolation des valeurs du fluide, selon l'angle entre la normale au solide et celle à la cellule de la grille cartésienne. Le schéma est simple à implémenter et précis à l'ordre 2. Nous présentons une étude de la convergence de l'erreur, des comparaisons avec d'autres méthodes de la littérature, et des exemples d'écoulements en une, deux et trois dimensions.

Abstract :

We present a second-order finite-volume scheme for the resolution of Euler compressible equations in complex geometries. The discretization is performed on a Cartesian grid, which in general does fit with the considered geometry. In the fluid domain we use a classical finite-volume method. On the cells located near the boundary with the solid, we solve an ah hoc Riemann problem taking into account the boundary conditions that we want to impose. To avoid pressure oscillations near the solid we balance the boundary condition with an extrapolation of the fluid values, as a function of the angle between the normal to the solid and the normal to the cell. The scheme is straightforward to implement and secondorder accurate. We present a convergence study, comparisons with other methods of the literature, and examples of flows simulations in one, two and three dimensions.

Mots clefs : Euler equations ; complex geometries ; cartesian grids

Introduction

The computation of flows in complex unsteady geometries is a crucial issue to perform realistic simulations of physical or biological applications like for instance biolocomotion (fish swimming or insect flight), turbomachines, windmills... To this end several class of methods exist. Here we are concerned with immersed boundary methods, i.e., integration schemes where the grid does not fit the geometry. These methods have been widely developed in the last 15 years, though the first methods were designed earlier (see for example [START_REF] Peskin | The fluid dynamics of heart valves : experimental, theoretical and computational methods[END_REF], [START_REF] Berger | An adaptive cartesian mesh algorithm for the Euler equations in arbitrary geometries[END_REF], [START_REF] Berger | Stable boundary conditions for cartesian grid calculations Computing systems in[END_REF]). The general idea behind immersed boundary methods is to take into account the boundary conditions by a modification of the equations to solve, either at the continuous level or at the discrete one, rather than by the use of an adapted mesh. The main advantages of using these approaches, compared to methods using body-conforming grids, are that they are easily parallelizable and allow the use of powerful line-iterative techniques. They also avoid to deal with grid generation and grid adaptation, a prohibitive task when the boundaries are moving. A recent through review of immersed boundary methods is provided by Mittal and Iaccarino [START_REF] Mittal | Immersed Boundary Methods[END_REF].

In this paper we present a simple globally second order scheme inspired by ghost cell approaches to solve compressible inviscid flows [START_REF] Gorsse | A simple second order cartesian scheme for compressible Euler flows[END_REF]. In the fluid domain, away from the boundary, we use a classical finite-volume method based on an approximate Riemann solver for the convective fluxes. At the cells located on the boundary, we solve an ad hoc Riemann problem taking into account the relevant boundary condition for the convective fluxes by an appropriate definition of the contact discontinuity speed.

To avoid pressure oscillations near the solid we balance the boundary condition with an extrapolation of the fluid values, as a function of the angle between the normal to the solid and the normal to the cell. Our objective is to device a method that can easily be implemented in existing codes and that is suitable for massive parallelization.

In section 2 we describe the finite volume scheme used in the fluid domain, away from the interface. In section 3 we introduce our method for solving the compressible Euler equations near the solid. Finally, in section 4 we present numerical tests to validate the expected order of convergence and to discuss performance compared to others immersed boundary or body fitted methods.

Resolution in the fluid domain 2.1 Governing equations

The compressible Euler equations are :

∂ρ ∂t + ∇ • ρu = 0 (1)
∂ρu ∂t + ∇ • (ρu ⊗ u + pn) = 0 (2) ∂E ∂t + ∇ • ((E + p)u) = 0 (3) (4)
where E denotes the total energy per unit volume. For a perfect gas

E = p γ -1 + 1 2 ρu 2 and p = ρRT (5)

Discretization

We focus on a two-dimensional setting. Let i and j be integers and consider the rectangular lattice generated by i and j, with spacing h x and h y in the x and y direction, respectively.

Let W be the conservative variables, F x (W ), F y (W )) the convective flux vectors in the x and y directions, respectively. By averaging the governing equations over any cell of the rectangular lattice we have

dW i j dt + 1 h x (F x i+1/2 j -F x i-1/2 j ) + 1 h y (F y i j+1/2 -F y i j-1/2 ) = 0 (6)
where W i j is the average value of the conservative variables on the cell considered, F x i+1/2 j the average flux in the x direction taken on the right cell side, and similarly for the other sides. The average convective fluxes at cell interfaces are approximated using the Osher numerical flux function. A second order Runge-Kutta scheme is used for the time integration.

A second order impermeability condition

For Euler equations, the boundary condition on the interface is the impermeability assumption, i.e., given normal velocity to the boundary (zero for a static wall, but non-zero for a moving solid). We are concerned with recovering second order accuracy on the impermeability condition.

Level set method

In order to improve accuracy at the solid walls crossing the grid cells we need additional geometric information. This information, mainly the distance from the wall and the wall normal, is provided by the distance function. The level set method, introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed : Algorithms based on HamiltonÐJacobi formulations[END_REF], is used to implicitly represent the interface of solid in the computational domain. The zero isoline of the level set function represents the boundary Σ of the immersed body. The level set function is defined here by :

ϕ(x) = dist Σ (x) outside of the solid -dist Σ (x) inside of the solid (7)
A useful property of this level set function is :

n(x) = ∇ϕ(x) (8)
where n(x) is the outward normal vector of the isoline of φ passing on x. In particular, this allows to compute the values of the normal to the interface, represented by the isoline ϕ = 0.

The impermeability condition in one dimension

The typical situation in one dimension for a grid that does not fit the body is shown in Fig. 1. For a fixed body, we want to impose u b = u(x b ) = 0 at the boundary point x b where ϕ(x b ) = 0.

u * φ=0 i i+1 i-1 u i-1
Figure 1 -Mesh near the solid. Red line is ϕ = 0. The interface lies between the center of cell i (fluid) and the center of cell i + 1 (solid). The flux in i + 1/2 has to be modified in order to account for the boundary conditions. Let u * be the contact discontinuity speed resulting from the solution of the Riemann problem defined at the interface between cell i and cell i + 1. We define a fictitious fluid state to the right of i + 1/2 such that u * takes into account, at the desired degree of accuracy, the boundary condition u b in x b . In particular, the value of u * at x i+1/2 is determined by the following formula :

u * = u b + 1 2 -d s b , (9)
where d = ϕ i ∆x , u b is the velocity of the obstacle (= 0 for a steady body) and s b is a slope defined by :

s b = d s 1 + (1 -d)s 2 , (10)
where

s 1 = u b -u i d and s 2 = u b -u i-1
1 + d are the slopes represented in Fig. 1. This formulation has been chosen to avoid stability issues in the slope computation when x b is close to x i since

s b = u b -u i + 1 -d 1 + d (u b -u i-1 ) . ( 11
)
If a slope limiter is needed, the limited slope is defined by s l b = minmod(s b , s 3 ), where s 3 = u i -u i-1 . An higher polynomial reconstruction to retrieve s b does not significantly modify the convergence rates shown in the following. At x i+1/2 the left fluid state of the Riemann problem U -= (u -, p -, c -) is computed as usual with the MUSCL reconstruction. The right state is U + = (-u -+ 2u * , p -, c -), so that u * is the contact discontinuity speed of the resulting Riemann problem. The left and right states of the variables p and c are identical in order to impose the correct wave reflection.

The scheme is non conservative at the numerical interface x i+1/2 . However, the loss of conservativity at this point is negligible compared to all the other points and the shocks are correctly resolved as shown in the following.

The impermeability condition in two dimensions

In two dimensions the flow equations are solved by computing independently the flux in each direction, so we want to apply in each direction the same kind of ideas as in one dimension. The interface points are the intersections between the interface (ϕ = 0) and the segment connecting the two cell centers concerned by the sign change (for example the points A, B and C on Fig. 2(a)). For the flux computation, a fictitious state is created for instance between the cells (i, j) and (i + 1, j) on Fig. 2(a). The boundary condition that we have to impose now is u A .n A = 0, where u A is the speed of the fluid at the boundary, and n A the outward normal vector of the body. When n A is parallel to the cell side normal n cell we recover the one-dimensional case. However, when the scalar product between the normal to the physical boundary and the normal to the mesh side is close to zero, we assume that the boundary condition only weakly affects the numerical flux, as it would happen for a body fitted mesh. The state corresponding to a cell in the solid is therefore modified in order to take into account this requirement thanks to a convex interpolation.

B A (i,j) (i+1,j) φ=0 n A n cell (a) Geometric configuration (u .n )n (u . ) u u i-1 j n * i-1 j i-1 j n T φ=0 (i, j) (i+1, j) (i-1, j)
We compute as before the left state primitive variables u -, p -, c -relative to the Riemann problem at the concerned numerical interface by a standard MUSCL reconstruction. According to Fig. 2(b) and considering the (i + 1/2, j)-flux, the right state will be U + = (-u -+ 2u w , v w , p w , c w ), where :

-

    u w v w p w c w     = α     u * v * p - c -     + (1 -α)     u f v f p f c f     ; -α = n A .n cell ; -(u f , v f , p f , c f ) is an extrapolated fluid state detailed hereafter. and U -= (u -, v w , p -, c -).
The speeds u * and v * : The normal component of the contact discontinuity speed u * is calculated as in the one-dimensional case. We determine the value of the contact discontinuity speed u * , relative to a Riemann problem defined in the direction normal to the cell side through x i+1/2,j , consistent at second order accuracy with u A • n A = 0 in A. The vector u * is determined as follows :

u * .n A = u * n = u A .n A + 1 2 -d s n A u * .τ A = u * τ = u -.τ A      ⇒ u * = u * n n x + u * τ τ x u * n n y + u * τ τ y
Where u A is the velocity of the obstacle (= 0 for a steady body), n A = (n x , n y ) t and τ A = (τ x , τ y ) t are respectively the normal and tangential vectors to the boundary at point A and the slope s n A is defined as in 1D :

s n A = u A .n A -u i .n A + 1 -d 1 + d (u A .n A -u i-1 .n A ) . (12)
If a slope limiter is needed, the limited slope is defined by s n,l A = minmod(s n A , s n 3 ), where

s n 3 = u i .n A - u i-1 .n A .
Bordeaux, 26 au 30 août 2013 Extrapolated fluid state : (u f , v f , p f , c f ) is computed as a linear extrapolation of the fluid variables to the point x i+1/2 j using the variables and the slopes of the closest upstream fluid cell.

Numerical illustrations

The Ringleb flow The Ringleb flow refers to an exact solution of Euler equations. The solution is obtained with the hodograph method. In our test case, the computational domain is [-0.5; -0.1] × [0; 0.6] and we numerically solve the flow between the streamlines Ψ 1 = 0.8 and Ψ 2 = 0.9. The inlet and outlet boundary condition are supersonic for y = 0 and y = 0.6 respectively.

The convergence orders are calculated in L 2 norm on four different grids 32 × 48, 64 × 96, 128 × 192 and 256 × 384. The results for the L 2 norm of our method are compared to a simple symmetry technique, the ghost-cell CCST method [START_REF] Dadone | Ghost-cell method for analysis of inviscid three dimensionnal flows on cartesian grids[END_REF], that relies on a local isoentropic flow model at the wall, and a standard finite-volume scheme case with a body-fitted meshs in Fig. 3. The overall results show that the classical symmetry scheme is first order accurate in the L 2 norm. The other schemes have overall comparable accuracy, although the amplitude of the error is lower for the present scheme compared to the ghost-cell CCST method. For the same test case, Coirier and Powell [START_REF] Coirier | An accuracy assessment of cartesian mesh approaches for the Euler equations[END_REF] observed also a convergence order between one and two in the case of their own cartesian method.

Mach 10 shock over three spheres The computation of a planar shock reflecting over three spheres is performed. The spheres are located at (0,0,0), (-0.75, 1.4,-1.4), (0,-2,0) with radii 1, 0.4, 0.3, respectively. The size of the domain is [-2.5; 2.5] 3 . The numerical computation is performed on a 256 3 mesh. Four snapshots are shown on Fig. 4. The interactions of the bow shocks give rise to complicated flow structures in the wake of the spheres.

Conclusion

We have presented a new cartesian method to solve compressible flows in complex domains. This method is based on a classical finite volume approach, but the values used to compute the fluxes at the cell interfaces near the solid boundary are determined so to satisfy the boundary conditions with a second order accuracy. Several numerical validations have been presented to assess its accuracy. This method is particularly simple to implement, as it doesn't require any special cell reconstruction at the solid-wall interface.

Figure 2 -

 2 Figure 2 -Fig. 2(a) : Example of geometric configuration at the interface. A is the interface point located between (i, j) and (i + 1, j). The flux on cell interface (i + 1/2, j) is modified to enforce the boundary condition on A. Fig. 2(b) : Graphical illustration of the construction of the u * vector.

Figure 3 -

 3 Figure 3 -Comparison of the L 2 accuracy of the present scheme with several methods. The convergence orders are detailed in the legend.

Figure 4 -Mach 10 planar shock reflecting on a 3D sphere. 20 isopressure surfaces.

Références


F Wurtz

B Delinchant

Xavier Brunotte

Julien Pouget

Dinh Van Binh

Les enjeux de la conception en phase d'esquisse pour les systèmes du génie électrique : illustration sur le cas des systèmes énergétiques pour les bâtiments

Keywords: conception, dimensionnement optimisation, esquisse énergétique, approche système, stockage, système hybrides, conception intégrée taille des systèmes & stratégie de supervision

RESUME

-Ce papier propose de donner un statut, et un outil, aux phases de conception initiale, dans lesquelles les caractéristiques principales du système énergétique utilisant l'énergie électrique doivent être esquissées. Nous proposons ainsi, d'utiliser un type de modèle adapté (modèle macroscopique et statique basé sur des bilans énergétiques et financiers). Ces modèles sont réalisés à l'échelle des composants du système et ils comportent des équations et des contraintes, qui sont ensuite composées à l'échelle du système, via des ports énergétiques et financiers. Il en résulte un modèle qui peut dès lors, être utilisé par des techniques d'optimisation. Cela permet de fournir un dimensionnement global et simultané des caractéristiques du système (taille et dimension des systèmes énergétiques) tout en esquissant une stratégie de gestion optimale sur des cycles de fonctionnement type. La fonction objectif est un coût global du système intégrant les coûts d'investissement, de maintenance et d'exploitation. Cette approche est illustrée sur un exemple de gare, représentatif à la fois des enjeux de la filière bâtiment, mais aussi représentatif de cette problématique de conception générique des systèmes hybrides comportant des sources, des charges et des moyens de stockage, auquel le génie électrique est de plus en plus confronté via des application dans les transports terrestres, aéronautiques, le smart-building, le smart-grid, … Mots-clés

-des approches d'optimisation comme aide à la décision -des modèles d'esquisse énergétique, qui doivent être d'une nature (hypothèses et niveaux de modélisations) et qui doivent porter des informations (comme des contraintes) compatibles avec l'objectif d'esquisse.

On illustrera ainsi comment ce type d'outil et d'approche permet d'« esquisser » simultanément les caractéristiques principales du système énergétique et la stratégie optimale de gestion énergétique. Ce concept et ce type d'outil ont été imaginés et appliqués ici pour des gares à énergie positive, mais les concepts et outils proposés s'avèrent être génériques et utilisables pour système ayant vocation à utiliser de l'énergie électrique.

L'ESQUISSE ENERGETIQUE : UN CONCEPT GENERIQUE MIS EN RELIEF DANS LE CONTEXTE DE LA CONCEPTION DES SYSTEMES BATIMENTS

Si le concept d'esquisse énergétique est général aux systèmes hybrides que l'on peut rencontrer en génie électrique, nous nous proposons néanmoins d'en montrer l'importance dans ce papier en s'appuyant sur les enjeux et la conception des systèmes bâtiments.

Ce n'est sans doute pas un hasard non plus, si cette notion d'esquisse énergétique a été formalisée dans notre groupe de recherche dans la thématique bâtiment, tant il est vrai que dans cette filière, la notion d'esquisse est valorisée, en faisant en particulier la noblesse du métier d'architecte. Or nous comptons montrer que cette mise en valeur doit dépasser la simple filière bâtiment et concerner tous les domaines de la conception en général, et celle des systèmes énergétiques en particulier.

IMPORTANCE, PLACE ET ROLE DE LA PHASE D'ESQUISSE DANS LE PROCESSUS DE CONCEPTION

Les phases d'esquisses sont les phases du processus de conception qui font partie des phases préliminaires. Les études théoriques sur l'activité de conception [VIS-04] montrent que, si ces phases préliminaires ne représentent qu'un faible coût (5% du coût des projets), elles sont pour autant fondamentales, car on y prend les décisions qui figent 75% des coûts totaux du projet. Il est donc important d'y d'explorer le maximum de possibilités, avec des outils de modélisation et d'aide à la décision pertinents, afin de fixer les bons choix dès les phases initiales, et en intégrant au plus tôt des aspects habituellement décidés très tard dans le cycle de conception des produits, comme la stratégie de gestion optimale. Par ailleurs, la montée en efficacité dans cette phase du processus de conception nécessité de franchir les verrous, d'une prise en compte dès les phases initiales du processus de conception :

-d'un maximum de paramètres et de contraintes, -d'un maximum de composants, -avec une vue « système et couplée » qui s'illustre à titre d'exemple dans la filière bâtiment par la volonté de concevoir au plus tôt simultanément l'enveloppe + les systèmes énergétiques + la stratégie de pilotage. Pour cette phase d'esquisse énergétique, nous proposons d'introduire pour chaque composant du système, un type de modèle dédié à l'optimisation énergétique qui comporte une dimension physique, avec pour chaque composant des équations reliant les principales dimensions physiques aux performances énergétiques et flux énergétiques (pour l'énergie électrique, l'énergie thermique, …). Ceci est réalisé par des modèles analytiques et statiques. Il n'est pas utile et pertinent, dans cette phase d'esquisse, de lancer des outils de simulation fine car : le système et les composants ne sont pas encore connus finement, il existe un besoin d'estimation rapide pour explorer le maximum de combinaisons et de possibilités. Il s'agit aussi d'introduire une dimension économique en estimant le coût économique sur la durée de vie prévue pour le composant sera la somme des coûts d'investissement et des coûts d'exploitation (ce qui est une manière d'esquisser un calcul d'Analyse de Cycle de Vie, même si le démantèlement n'est pas intégré, et doit plutôt faire l'objet d'ACV ultérieures [START_REF]L'éco-conception dans le génie électrique : application au cas du transformateur[END_REF] sans doute dans une phase plus avancée du processus de conception).

Dans ce but, chaque composant possible du système doit fournir des équations pour calculer : les coûts d'investissement, les coûts d'exploitation et des contraintes qui peuvent porter sur : la taille des équipements (cf. la taille du panneaux solaire qui doit rester inférieure à une surface max qui est typiquement la surface du toit), le fonctionnement des équipements (cf. l'état de charge de la batterie qui doit rester inférieur à la capacité max de la batterie), des ports qui doivent permettre la composition à l'échelle du système au niveau des flux énergétiques et économiques.

Exemple d'un modèle de batterie dédié à l'esquisse énergétique

La figure 4 illustre la structure d'un tel modèle pour l'exemple d'une batterie électrique. Nous allons détailler les équations et contraintes définies, avec les hypothèses associées, de sorte à avoir un modèle permettant le dimensionnement de cette batterie en phase d'esquisse, et la détermination de la stratégie optimale de gestion de cette batterie sur 1 journée de 24h discrétisée par pas de 1H (ce qu'on appellera par la suite un cycle). Équations du modèle d'esquisse de la batterie Calcul de l'état de charge de la batteriebilan énergétique du composant batterieconnexion au port énergétique du système On détermine à chaque heure l'état de charge de la batterie (en W) en fonction de l'état de charge à l'heure précédente et du bilan de puissance entrée-sortie de la batterie modulé par son rendement. Il s'agit typiquement d'une équation de bilan d'énergie typique de celle que l'on peut trouver dans les modèles d'esquisses que nous avons développés :

[eqbat1] : Pour i Є[0;23] : Bat_EnergieStockee[i+1]=Bat_EnergieStockee[i] +Bat_EnergieEntrante[i]*Rendement_convertisseur Bat_EnergieSortante[i]/Rendement_convertisseur C'est

Contraintes du modèle d'esquisse de la batterie

Le modèle est destiné à être utilisé avec un algorithme d'optimisation, qui aura pour objectif de chercher simultanément un dimensionnement de la batterie et une stratégie de gestion optimale sur cycle de fonctionnement (stratégie de pilotage). Il est en résulte que le modèle doit comporter, en plus des équations précédentes de bilan énergétique et financier, des contraintes liées à la recherche de la stratégie de pilotage optimale et au dimensionnement optimal. Contraintes liées à la recherche de la stratégie de gestion optimale Contrainte à annuler pour trouver lʹétat de charge initial optimal de la batterie

On définit en fait ici un ensemble d'une équation et d'une contrainte permettant de faire le lien entre les états de charge de la batterie aux heures 0 et 23 :

[eqbat4] : Bat_EnergieStockee_Annuler=Bat_EnergieStockee[0] -Bat_EnergieStockee[23] +Bat_EnergieEntrante[23]*Rendement_convertisseur -Bat_EnergieSortante[23]/Rendement_convertisseur)
On associera donc à cette équation la contrainte :

[cbat1] : Bat_EnergieStockee_Annuler=0 Cet ensemble équation+contrainte impose ainsi à l'optimiseur de trouver une stratégie de pilotage de la batterie sur un jour de fonctionnement avec une condition de cyclicité, permettant typiquement de reproduire le cycle de fonctionnement sur un grand nombre de jours. On pourra ainsi supposer qu'il peut y avoir dans l'année N cycles de fonctionnement identiques, avec typiquement N=4 pour un cycle printemps + cycle été + cycle automne + cycle hiver. Ces cycles seront reproduits le nombre de fois nécessaires pour couvrir une année typique.

Contraintes limitant la puissance extraite/injectée dans la batterie

Ces contraintes font partie des contraintes de fonctionnement classiques d'une batterie :

[cbat2] :

Pour i Є[0;23] : 0 < Bat_EnergieEntrante[i] < 2000 kW [cbat3] : Pour i Є[0;23] : 0 < Bat_EnergieSortante[i] < 2000 kW
Cette puissance est ici à titre d'exemple limitée à 2000 kW, mais au lieu d'avoir une valeur fixe, on peut définir un paramètre qui peut être introduit dans la liste des paramètres du tableau 1.

Contraintes liant le dimensionnement de la batterie et le cycle de fonctionnement optimal

Contraintes imposant que l'état de charge de la batterie soit inférieur à sa capacité Il s'agit de la contrainte suivante :

[cbat4] : 0 < Bat_EnergieStockee[i]< Bat_Capacite La prise en compte de cette contrainte par l'algorithme d'optimisation va permettre de lier la recherche de taille optimale de la batterie (Bat_Capacite) à la recherche du cycle optimale de gestion de celle-ci (définie par Bat_EnergieStockee

[i], elle-même lié à Bat_EnergieEntrante[i] et Bat_EnergieSortante[i] via [eqbat1]).
Le tableau 1 donne les paramètres définis pour le modèle de batteries (ce sont des variables constantes au cours d'une optimisation, mais qui peuvent varier d'une optimisation à l'autre), tandis que le tableau 2 donne la liste des variables d'optimisation. Cette taille de problème implique d'utiliser des algorithmes détermines de type MILP ou SQP qui sont aptes à gérer un si grand nombre de contraintes et de paramètres en un temps de réponse rapide (le temps étant un paramètre important en phase d'esquisse ou le concepteur se doit de multiplier les essais et les esquisses). Une hybridation avec des algorithmes stochastiques peut toutefois être considérée, en cas de présence de problèmes de convergence dépendant du point du départ (pour des algorithmes de types SQP), ou en cas de problèmes comportant des optima locaux, mais la convergence finale ne peut être assurée que par un algorithme déterministe assurant dans tous les cas une gestion fine des contraintes malgré leur grand nombre. La figure 11 montre la stratégie de gestion optimale de la batterie trouvée pour c_inv_PV=200 €/m² et c_inv_PV=450 €/m². On y observera que :

-La stratégie de gestion de la batterie est toujours déterminée pour plutôt charger de l'électricité du réseau lorsque celle-ci est peu chère, et la décharger localement lorsque le prix de cette électricité du réseau monte -Que la taille optimale de la batterie préconisée est beaucoup plus grande lorsque le coût PV est grand

  1. INTRODUCTION DE LA NOTION D'ESQUISSE ENERGETIQUE POUR LE BATIMENT : DEFINITION CONCEPTUELLE ET PROPOSITION D'UN TYPE D'OUTIL AD'HOC On se propose d'introduire un concept innovant d'esquisse énergétique. Il s'agit d'un concept équivalent à ce qui peut exister en architecture ou en design [HAU-01], mais focalisé sur l'aspect énergétique. On propose, pour instrumenter cette phase, de définir un concept d'outil d'esquisse énergétique utilisant:

Fig. 1 :

 1 Fig. 1 : Importance, place de l'esquisse dans les processus de conception 3.1.1. Illustration de cette importance dans la filière bâtiment L'exemple de la filière bâtiment est particulièrement percutant pour comprendre les enjeux et la complexité relatifs à cette phase d'esquisse, dès lors qu'on rappelle : -qu'un bâtiment est un investissement lourd -réalisé pour du long terme (un bâtiment doit durer plus de 30 ans) -que dès les phases d'esquisses, il conviendrait : o de ne pas se tromper dans la multiplicité des choix de conception possibles sur par exemple les choix de systèmes de production énergétique (PV, co-générateurs, chaudière gaz, chaudière bois, …) o d'avoir une vue globale du coût de type Analyse de cycle de vie incluant l'investissement initial, les coûts d'exploitation et de maintenance, et si possible de démantèlement o en gérant toutes les incertitudes liées à l'évolution des couts des technologies et de l'énergie -ce d'autant plus que la filière bâtiment représente 42% de l'énergie primaire consommée en France et 64% de la consommation d'électricité [PEU]. Il convient donc d'être efficace au plus tôt et dès les phases d'esquisse vu le poids économique, énergétique et environnemental de cette filière

Fig. 3 :

 3 Fig. 3 : Verrou de la phase esquisse : concevoir au plus tôt avec une vue « système et couplée » 5. INTRODUCTION D'UN TYPE «MODELE D'ESQUISSE ENERGETIQUE» POUR LES COMPOSANTS DU SYSTEME

Fig. 4 :

 4 Fig. 4 : Structure des modèles de composants réalisés pour un outil d'esquisse énergétique

  par cette équation que la batterie est énergétiquement reliée au système électrique, à chaque pas horaire du cycle, par les portsBat_EnergieEntrante[i] etBat_EnergieSortante[i] Calcul du bilan et du coût financier de la batterieconnection aux ports économiques du système On va supposer que le coût de la batterie (investissement+maintenance) dépend proportionnellement de la puissance qui en est extrait sur le nombre de cycle que l'on va considérer : c'est un exemple typiquement de modèle ad'hoc en phase d'esquisse, offrant un juste compromis qui ne soit ni trop fin, ni trop grossier. On fait donc la somme heure par heure des puissances sortantes. On multiplie ensuite par le nombre de jour de l'étude (le nombre de cycle) et par le coefficient de proportionnalité permettant d'aboutir au prix final. [eqbat2] :Pour i Є[0;23] : Interne_Qe= NbCycle*Σ Bat_EnergieSortante [i] [eqbat3] : cout_inv_bat=Interne_Qe*c_inv_bat On introduit ainsi la dimensionnement économique et financière. Chaque composant du système doit ainsi renvoyer des informations relatifs à son coût d'investissement et à son coût d'exploitation, constituant les ports économiques que ces composants doivent offrir pour pouvoir calculer un coût d'investissement et d'exploitation global à l'échelle du système. Ici le port économique est couv_inv_bat intégrant le coût d'investissement et d'exploitation pour la batterie.

Figure 5 :

 5 Figure 5 : Composition des modèles (équations + contraintes) à l'échelle du système via les ports (énergétiques et financiers) 8. L'EXEMPLE CONSIDERE: DIMENSIONNEMENT SIMULTANE DE L'ENVELOPPE + LA TAILLE DES SYSTEMES ENERGETIQUES +LA STRATEGIE DE SUPERVISION OPTIMALE L'exemple qui sera présente a été construit pour un bâtiment type gare pouvant être équipé de systèmes énergétiques tels que des panneaux solaires, des cogénérateurs, des batteries, … Il possède un réseau thermique de chauffage et il est connecté au réseau électrique avec une courbe de charge électrique pouvant comporter des capacités de report. Le modèle et l'application ont été développés dans le framework CADES initialement développé au G2ELAB [DEL-07] et industrialisé à présent par Vesta-system (cf. http://www.cades-solutions.com). Ceci a permis de réaliser l'outil d'optimisation en phase d'esquisse qui est visualisé sur la figure 6 selon le processus détaillé sur la figure 7°:

Fig. 6 :

 6 Fig. 6 : Outil pour la phase d'esquisse énergétique réalisant permettant le dimensionnement simultané du bâtiment, de la taille des systèmes énergétiques et de la stratégie de supervsion

  10. DESCRIPTION D'ETUDES ET DE RESULTATS QUE L'ON PEUT PRODUIRE L'OUTIL D'ESQUISSE ENERGETIQUE PROPOSE On montre ici des études et des analyses qui peuvent typiquement réalisées avec le type d'outil d'esquisse que nous proposons. L'objectif ici n'est pas de nous focaliser sur les valeurs quantitatives, que nous ne présenterons pas en détails, ce qui nous emmènerait dans des discussions sur les scénarios et les hypothèses prises en compte, ce qui n'est pas l'objet de ce papier. Notre prétention est ici simplement de montrer le type de résultat et d'étude qu'un concepteur peut mener en phase d'esquisser de sorte à pouvoir prendre des décisions au plus tôt, qui soient le mieux éclairées possible et avec une vision qui soit la plus système et la plus intégrée possible. Sur toutes les optimisations réalisées, la fonction objectif était un coût global du bâtiment intégrant le coût d'investissement, d'exploitation et de maintenance sur une durée de vie de 30 ans. Avec l'approche d'esquisse proposée nous sommes aptes à produire rapidement des résultats comme ceux visualisés sur la figure 7, issus d'une série d'optimisations paramétrées montrant la surface optimale de PV dans laquelle investir en fonction de la technologie PV que l'on fait varier entre 200 et 1100 €/m² avec 20 optimisations réalisées entre ces 2 extrêmes. On notera tout d'abord qu'avec l'approche d'optimisation utilisée (SQP + Calcul de gradient du modèle), l'optimisation est rapide : de l'ordre de 1 minutes pour les 20 optimisations, en assurant une convergence fine (respect précis des contraintes).Ce type d'optimisation paramétré est très intéressant pour gérer l'incertitude, typiquement ici l'incertitude sur le prix de la technologie PV : sachant que le concepteur a cette incertitude à gérer, et qu'il devra prendre des risques en conséquences, une bonne manière d'éclairer sa décision, est d'étudier tous les scénarios de prix possibles, pour ensuite prendre sa décision en fonction des résultats de chacun de ces scénarios.Ainsi ces optimisation paramétrés font apparaitre des effets de seuils : pour un coût PV entre 200€/m² et 450 €/m², l'optimisation nous indique que la surface optimale de PV est de 178 m². Dans ces hypothèses de prix (les autres hypothèses étant constantes par ailleurs), l'optimisation d'esquisse nous indique donc que nous avons intérêt à maximiser la production de production de photovoltaïque en quantifiant l'ordre de grandeur de PV à installer.

Fig. 8 :

 8 Fig. 8: Taille optimale de panneaux solaire (en m²) en fonction du coût de la technologie de panneaux solaire (en k€/m²) -Chaque point correspond à une optimisation réalisée en faisant varier c_inv_PC entre 200 €/m² et 1100 €/m²

Fig. 10 :

 10 Fig. 10: Focus sur une optimisation typique dans la zone 450€/m² à 950 €/m² -La technologie co-générateur devient à présent compétitivé et permet de diminuer l'énergie prélévée sur le réseau (cf. partie hachurée consommée en moins)

Fig. 11 :

 11 Fig. 11: Exemples de stratégie de gestion de batterie (sur 1 journée type de fonctionnement) pour des prix extrêmes de la technologie PV -Les stratégies de gestion de la batterie sont des résultats d'optimisation, le vecteur de « Tarif dynamique utilisé » est un paramètre de l'optimisation La figure 12 donne un exemple d'évolution de la résistance thermique optimale du bâtiment côté est : on constate qu'elle diminue avec des effets de seuils qui semblent à priori peu compréhensibles. En fait il s'agit d'un effet système résultant du mécanisme suivant : dès lors que la technologie PV devient plus chère, la technologie du co-générateur devient plus intéressante, et à chaque que l'on introduit un co-générateur de taille plus grande, on bénéficie en même temps d'une production de chaleur. Cette production de chaleur permet d'assurer ensuite le même confort thermique en permettant un bâtiment moins isolé thermiquement, ce qui s'avère être l'optimum dans le cadre du modèle et des hypothèses formulées.

Fig. 12 :

 12 Fig. 12: Evolution de la résistance thermique des murs Les résultats des scenarios précédents n'ont pas vocation à être généralisés : ils ne peuvent l'être, car ils ont été établis dans le cadre de modèles d'esquisse et d'un nombre important d'hypothèses et de choix de paramètres que nous ne pouvons exposer ici. Leur seule ambition et de montrer le potentiel de l'approche en terme de capacité à construire rapidement des