Alessandro Lazaric
email: alessandro.lazaric@inria.fr

Mohammad Ghavamzadeh
email: mohammad.ghavamzadeh@inria.fr

Bayesian Multi-Task Reinforcement Learning

We consider the problem of multi-task reinforcement learning where the learner is provided with a set of tasks, for which only a small number of samples can be generated for any given policy. As the number of samples may not be enough to learn an accurate evaluation of the policy, it would be necessary to identify classes of tasks with similar structure and to learn them jointly. We consider the case where the tasks share structure in their value functions, and model this by assuming that the value functions are all sampled from a common prior. We adopt the Gaussian process temporal-difference value function model and use a hierarchical Bayesian approach to model the distribution over the value functions. We study two cases, where all the value functions belong to the same class and where they belong to an undefined number of classes. For each case, we present a hierarchical Bayesian model, and derive inference algorithms for (i) joint learning of the value functions, and (ii) efficient transfer of the information gained in (i) to assist learning the value function of a newly observed task.

Introduction

Multi-task learning (MTL) is an important learning paradigm and has recently been an area of active research in machine learning (e.g., [4; 1; 17; 16; 3]). A common setup is that there are multiple related tasks for which we are interested in improving the performance over individual learning by sharing information across the tasks. This transfer of information is particularly important when we are provided with only a limited number of data to learn each task. Exploiting data from related problems provides more training samples for the learner and can improve the performance of the resulting solution.

Most reinforcement learning (RL) algorithms [2; 13] often need a large number of samples to solve a problem and cannot directly take advantage of the information coming from other similar tasks. Nonetheless, recent work has shown that transfer and multi-task learning techniques can be employed in RL to reduce the number of samples needed to achieve nearly-optimal solutions. All approaches to multi-task RL (MTRL) assume that the tasks share similarity in some components of the problem such as dynamics, reward structure, or value function. While some methods explicitly assume that the shared components are drawn from a common generative model [15; 10], this assumption is more implicit in others [14; 9]. In [START_REF] Mehta | Transfer in variable-reward hierarchical reinforcement learning[END_REF], tasks share the same dynamics and reward features, and only differ in the weights of the reward function. The proposed method initializes the value function for a new task using the previously learned value functions as a prior. In [START_REF] Wilson | Multi-task reinforcement learning: A hierarchical Bayesian approach[END_REF], the distribution over the dynamics and the reward functions of the tasks is drawn from a hierarchical Bayesian model (HBM). Due to some similarity to our work, we discuss this method in more details in Section 5. In [START_REF] Lazaric | Transfer of samples in batch reinforcement learning[END_REF], the authors implicitly assume that the tasks are drawn from a common distribution. They propose a method to selectively transfer samples from source tasks to a target task based on the likelihood of the target samples being generated by the models built for the source tasks. Finally, in [START_REF] Taylor | Transfer learning via inter-task mappings for temporal difference learning[END_REF], learning the value function of the target task is expedited using the solution learned in a source task with related, but different, state and action spaces.

In this paper, we study the MTRL scenario in which the learner is provided with a number of tasks with common state and action spaces. For any given policy, only a small number of samples can be generated in each task, which may not be enough to accurately evaluate the policy. In such a MTRL problem, it is necessary to identify classes of tasks with similar structure and to learn them jointly. In our work, we consider a particular class of MTRL problems in which the tasks share structure in their value functions. To allow the value functions to share a common structure, one way would be to assume that they are all sampled from a common prior. We adopt the Gaussian process temporal-difference (GPTD) value function model [START_REF] Engel | Reinforcement learning with Gaussian processes[END_REF] for each task, model the distribution over the value functions using a HBM, and develop solutions to the following problems: (i) joint learning of the value functions, and (ii) efficient transfer of the information acquired in (i) to facilitate learning the value function of a newly observed task. We refer to the above problems as symmetric and asymmetric multi-task learning, respectively. In Section 3, we present a HBM for the case in which all the value functions belong to the same class, and derive an EM algorithm to find MAP estimates of the value functions and the model's hyper-parameters. However, as pointed out in [4; 1], if the functions do not belong to the same class, simply learning them together can be detrimental (negative transfer). It is therefore important to have models that will generally benefit from related tasks and will not hurt performance when the tasks are unrelated. This is particularly important in RL as changing the policy at each step of the policy iteration algorithm (this is true even for the fitted value iteration algorithm) can change the way tasks are clustered together. This means that even if we start with value functions that all belong to the same class, after one iteration the new value functions may be clustered into several classes. In Section 4, we introduce a Dirichlet process (DP) based HBM for the case that the value functions belong to an undefined number of classes, and derive inference algorithms for both the symmetric and asymmetric scenarios. In Section 5, we discuss the similarities and differences with closely related work. In Section 6, we report and analyze experimental results.

Preliminaries

The agent-environment interaction in RL is conventionally modelled as a Markov Decision Process (MDP). A MDP is a tuple M = X , A, R, P where X and A are the state and action spaces, respectively; R is the probability distribution over rewards R; and P is the transition probability distribution. A stationary policy π : X × A → [0, 1] is a mapping from states to action selection probabilities. The MDP controlled by a policy π induces a Markov chain with transition probability distribution P π (x ′ |x) = A P (x ′ |x, a)π(a|x)da. Given a policy π, the (possibly discounted, γ ∈ [0, 1)) return for a state x, D π (x), is a random process defined by

D π (x) = ∞ t=0 γ t R(x t )|x 0 = x, with x t+1 ∼ P π (•|x t ). The value function V π (x) is the expected value of D π (x)
where the expectation is over all possible trajectories and rewards collected along them.

A key problem in RL is to learn the value function of a given policy, which is called policy evaluation [2; 13]. Loosely speaking, in policy evaluation the goal is to find a "close enough" approximation V of the value function V π . Unlike in supervised learning, the target function V π is not known in advance and its values have to be inferred from the observed rewards. Therefore, it is required to define a stochastic generative model connecting the underlying hidden value function with the observed rewards. In this paper, we adopt the GPTD value function model proposed in [START_REF] Engel | Reinforcement learning with Gaussian processes[END_REF], in which the discounted return D is decomposed into its mean V and a random zero-mean residual ∆V , D(x) = V (x) + ∆V (x). Combining this decomposition with the Bellman equation, we obtain

R(x) = V (x) -γV (x ′ ) + ǫ(x, x ′ ), x ′ ∼ P π (•|x), (1)
where ǫ(x, x ′ ) def = ∆V (x) -γ∆V (x ′ ). Suppose we are provided with a set of samples

D = {(x n , x ′ n , r n )} N n=1
, where r n and x ′ n are the reward and the next state observed by following policy π in state x n , respectively. By writing the model of Eq. (1) w.r.t. these samples, we obtain R = HV + E, where H ∈ R N ×2N and Note that if the samples are generated from a single trajectory then x ′ n = x n+1 and H will be of the form defined by Eq. (2.7) in [START_REF] Engel | Reinforcement learning with Gaussian processes[END_REF]. In order to specify a complete probabilistic generative model connecting values and rewards, we need to define a prior distribution for the value function V and the distribution of the noise ǫ. Similar to [START_REF] Engel | Reinforcement learning with Gaussian processes[END_REF], we model the value function as a Gaussian process (GP), and the noise vector as E ∼ N (0, S), where S is the noise covariance matrix modelling the correlation of the noise between different states. In the following we write S = σ 2 P, where σ 2 and P are the variance and the correlation matrix of the noise, respectively. For a more extended discussion about different models of noise we refer readers to [START_REF] Engel | Algorithms and Representations for Reinforcement Learning[END_REF] and Section 8.4. The value function V may be represented either in parametric or non-parametric form. In this paper we use the parametric representation to make the formulation easier to follow, but all the results can be extended to the nonparametric case following similar steps as in Section 5.2 of [START_REF] Yu | Learning Gaussian processes from multiple tasks[END_REF]. In the parametric form, the value function is represented by a finite set of d features φ(•) = φ 1 (•), . . . , φ d (•)

R ⊤ = (r n ) N n=1 ; E ⊤ = ǫ(x n , x ′ n ) N n=1 ; V ⊤ = V (x n ), V (x ′ n ) N n=1 (a) xmn N wm rmn σ 2 (µ, Σ) (α 0 , β 0 ) (µ 0 , k 0 , ν 0 , Σ 0 ) x ′ mn (b) xmn N wm rmn M σ 2 (α0, β0) (µ 0 , k0, ν0, Σ0) (µ, Σ) x ′ mn (c) xmn N cm wm rmn M ∞ vc σ 2 c (µ c , Σc) (α0, β0) τ (µ 0 , k0, ν0, Σ0) x ′ mn
⊤ and a weight vector

w = (w 1 , . . . , w d ) ⊤ as V (•) = φ(•) ⊤ w.
The randomness in V is now due to w being a random vector with Gaussian prior w ∼ N (µ, Σ). The model equation now becomes R = HΦ ⊤ w + E, where Φ = [φ(x 1 ), φ(x ′ 1 ), . . . , φ(x N ), φ(x ′ N )] is a d × 2N matrix. Fig. 1(a) shows the graphical representation of this model. Note that the model shown in Fig. 1(a) is the model that we use for single-task learning (STL) in the experiments of this paper. It is an extension of the original GPTD model by defining a Normal-inverse-Wishart and an inverse-Gamma hyper-priors parametrized by ψ 0 = (µ 0 , k 0 , ν 0 , Σ 0 , α 0 , β 0 ) over the model parameters (µ, Σ, σ 2 ). This allows us to optimize the model parameters given the data.

In the MTRL setting of this paper, the learner is provided with M tasks or MDPs with common state and action spaces M m = X , A, R m , P m , m = 1, . . . , M . Given a fixed policy, N samples are generated in each task, i.e., D m = {(x mn , x ′ mn , r mn )} N n=1 , which may not be enough to have an accurate evaluation of the policy. We consider the case in which the tasks share structure in their value functions. In the parametric value function model discussed above, this can be interpreted as the value functions share the same feature space and their weight vectors are sampled independently from a common prior, i.e., V m (•) = φ(•) ⊤ w m ; w m ∼ N (µ, Σ). In the next two sections, we study two different scenarios: 1) when all the tasks belong to the same class, i.e., they share the same prior, and 2) when they can be clustered into an undefined number of classes.

Single-class Multi-task Learning

In this section, we consider the case where all the tasks belong to the same class, i.e., they share the same distribution over their value functions w m ∼ N (µ, Σ), m = 1, . . . , M ; and the same observation noise σ 2 . The goal in the symmetric form of this problem is to estimate {w m } M m=1 from the data {D m } M m=1 , whereas in the asymmetric case we are interested in estimating the parameters θ = (µ, Σ, σ 2 ) from the data in order to use them as a prior for a newly observed task (e.g., task M M+1 ). We use a parametric HBM for this problem. HBMs allow us to model both the individuality of the tasks and the correlation between them. In HBMs, individual models with task specific parameters are usually located at the bottom, and at the layer above, tasks are connected together via a common prior placed over those parameters. Learning the common prior is a part of the training process in which data from all the tasks contribute to learning, thus making it possible to share information between the tasks usually via sufficient statistics. Then given the learned prior, individual models are learned independently. As a result, learning at each task is affected by both its own data and by data from the other tasks related through the common prior.

The Model

We assume a normal-inverse-Wishart and an inverse-Gamma hyper-priors for (µ, Σ) and σ 2 , respectively.

p(θ|ψ 0 ) = p(µ, Σ) × p(σ 2 ) = N (µ; µ 0 , Σ/k 0 ) IW(Σ; ν 0 , Σ 0 ) × IG(σ 2 ; α 0 , β 0 ). (2)
These distributions are the conjugate priors for multivariate Gaussian distributions p(w m |µ, Σ) and p(R m |w m , σ 2 ) = N (HΦ ⊤ m w m , σ 2 P), respectively. This leads to the following generative model for the data, {D m }. Fig. 1(b) shows the graphical representation of this model. The details of the model can be found in Section 8.

Single-Class Model: Given the hyper-parameters

ψ 0 = (µ 0 , k 0 , ν 0 , Σ 0 , α 0 , β 0 ),
1. The parameters θ = (µ, Σ, σ 2 ) are sampled once from the hyper-prior as in Eq. (2), 2. For each task M m (value function V m ), the weight vector is sampled as w m ∼ N (µ, Σ),

3. Given {(x mn , x ′ mn )} N n=1 , we have R m = HΦ ⊤ m w m + E, where E ∼ N (0, σ 2 P), m = 1, . . . , M .

Inference

This model can be learned by optimizing the penalized likelihood p {R m }|{(x mn , x ′ mn )}, θ p(θ) w.r.t. the parameters θ = (µ, Σ, σ 2 ) using an EM algorithm. In the rest of the paper, we refer to this algorithm as SCMTL, for single-class multi-task learning.

E-step: Since the posterior distribution of the latent variables p({w

m }|{D m }, θ) is a prod- uct of M Gaussian posterior distributions p(w m |D m , θ) = N (µ ′ 0m , Σ ′ 0m
), for each task m, we compute the mean and covariance as

µ ′ 0m = Σ ′ 0m 1 σ 2 Φ m H ⊤ P -1 R m + Σ -1 µ , Σ ′ 0m = 1 σ 2 Φ m H ⊤ P -1 HΦ ⊤ m + Σ -1 -1
.

M-step:

We optimize θ to maximize the penalized expected log-likelihood of complete data log p {D m }, {w m }|θ over the posterior distribution estimated in the E-step and obtain the new parameters

µ new = 1 M + k0 k0µ 0 + M m=1 µ ′ 0m , Σnew = 1 M + ν0 + d + 2 k0(µ -µ 0 )(µ -µ 0 ) ⊤ + Σ0 + M m=1 (µ ′ 0m -µ)(µ ′ 0m -µ) ⊤ + Σ ′ 0m , σ 2 new = 1 M N + 2(1 + α0) 2β0 + M m=1 tr P -1 HΦ ⊤ m Σ ′ 0m ΦmH ⊤ + Rm -HΦ ⊤ m µ ′ 0m P -1 Rm -HΦ ⊤ m µ ′ 0m ⊤ .

Multi-class Multi-task Learning

In this section, we consider the case where the tasks belong to an undefined number of classes. Tasks in the same class {M m | c m = c} share the same distribution over their value functions w m ∼ N (µ c , Σ c ), and the same observation noise σ 2 c . We use a nonparametric HBM for this problem. In the HBM proposed in this paper, the common prior is drawn from a Dirichlet process (DP). As addressed in the statistics literature (see, e.g., [START_REF] Mukhopadhyay | Dirichlet process mixed generalized linear models[END_REF]), DP is powerful enough to model the parameters of individual classes, to fit them well without any assumption about the functional form of the prior, and to automatically learn the number of underlying classes.

The Model

We place a DP(τ, G 0 ) prior over the class assignment and the class parameters. The concentration parameter τ and the base distribution G 0 can be considered as priors over the number of classes and the class parameters θ c = (µ c , Σ c , σ 2 c ), respectively. G 0 is specified as the product of a d-dimensional normal-inverse-Wishart and a 1-dimensional inverse-Gamma distributions, with parameters ψ 0 = (µ 0 , k 0 , ν 0 , Σ 0 , α 0 , β 0 ), (see Eq. 2). We employ the stick-breaking representation of the DP prior [START_REF] Ishwaran | Gibbs sampling methods for stick-breaking priors[END_REF], and define a task-to-class assignment variable (c m1 , . . . , c m∞ ) for each task m, whose elements are all zero except that the cth element is equal to one if task m belongs to class c. Given the above, the data {D m } can be seen as drawn from the following generative model, whose graphical representation is shown in Fig. 1(c).

Multi-Class Model:

Given the hyper-parameters (τ, ψ 0 ),

1. Stick-breaking view: Draw v c from the Beta distribution Be(1, τ ), compute π c = v c c-1 i=1 (1 -v i ), and independently draw θ c ∼ G 0 , c = 1, . . . , ∞, 2. Task-to-class assignment: Draw the indicator (c m1 , . . . , c m∞ ) from a multinomial distri- bution M ∞ (1; π 1 , . . . , π ∞ ), m = 1, . . . , M ,

The weight vector is sampled as

w m ∼ N (µ cm , Σ cm ), m = 1, . . . , M , 4. Given {(x mn , x ′ mn )} N n=1 , we have R m = HΦ ⊤ m w m + E
, where E ∼ N (0, σ 2 cm P), m = 1, . . . , M .

Inference

We are interested in the posterior distribution of the latent variables Z = {wm}, {cm}, {θc} given the observed data and the hyper-parameters τ and ψ 0 , i.e., p Z|{Dm}, τ, ψ0 ∝ p {Dm}|Z, τ, ψ0 p(Z|τ, ψ0). In the following we outline the main steps of the algorithm used to solve this inference problem, which we refer to as MCMTL, for multi-class multi-task learning (see Fig. 2). MCMTL combines the SCMTL algorithm of Sec. 3.2 for class parameters estimation, with a Gibbs sampling algorithm for learning the class assignments [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF]. The main advantage of such combination is that at each iteration, given the current estimate of the weights, we take advantage of the conjugate priors to derive an efficient Gibbs sampling procedure.

More formally, given an arbitrary initial class assignment {c m }, a distinct EM algorithm is run on each class c = 1, . . . , C (with C the current estimate of the number of classes) and returns M distributions N (µ ′ 0m , Σ ′ 0m ). Given the weights estimated at the previous step, w m = µ ′ 0m , the Gibbs sampling solves the DP inference by drawing samples from the posterior distribution p {c m }|{R m }, { w m }, τ, ψ 0 . In particular, the state of the Markov chain simulated in the Gibbs sampling is the class assignment {c m }, i.e., the vector of the classes each task belongs to. At each iteration, each component c m of the state is updated by sampling from the following distribution where M -m,c is the number of tasks in class c except task m, and b is a normalizing constant. While the first term in Eq. ( 3) is the probability of task m to belong to an existing class c, the second term returns the probability of assigning task m to a new class. Thanks to the conjugate base distribution G 0 , the integrals in Eq. ( 3) can be solved analytically. In fact

If c = c m ′ , m ′ = m: p cm = c|{c m ′ }, Rm, wm, τ, ψ0 = b M-m,c M -1 + τ p Rm, wm|θc p θc|{c m ′ }, ψ0 dθc, else: p cm = c m ′ , m ′ = m|{c m ′ }, Rm, wm, τ, ψ0 = b τ M -1 + τ p(Rm, wm|θ)p(θ|ψ0)dθ, (3)
p(R m , w m |θ)p(θ|ψ 0 ) = p(R m | w m , σ 2 )p( w m |µ, Σ) × p(µ, Σ|µ 0 , k 0 , ν 0 , Σ 0 )p(σ 2 |α 0 , β 0 ) ∝ N (µ; µ ′ 0 , Σ/k ′ 0 )IW(Σ; ν ′ 0 , Σ ′ 0 ) × IG(σ 2 ; α ′ 0 , β ′ 0 ), (4)
where

ψ ′ 0 = (µ ′ 0 , k ′ 0 , ν ′ 0 , Σ ′ 0 , α ′ 0 , β ′ 0 )
are the posterior parameters of G 0 given the weight w m and the rewards R m (see Section 8 for their definition). Using the posterior hyper-parameters, the second integral in Eq. ( 3) can be written as

p(R m , w m |θ)p(θ|ψ 0 )dθ = k 0 πk ′ 0 d 2 |Σ 0 | ν0/2 |Σ ′ 0 | ν ′ 0 /2 Γ ν ′ 0 2 Γ ν ′ 0 -d 2 ×(2π|P|) -N 2 β α0 0 Γ(α 0 ) Γ(α ′ 0 ) β ′α ′ 0 0 . (5)
In the first integral of Eq. ( 3), the density function p(θ c |{c m ′ }, ψ 0 ) is the posterior probability over the class parameters θ c given the data from all the tasks belonging to c m according to the current class assignment {c m ′ }. Similar to Eq. ( 4), we compute the posterior hyper-parameters ψ 0c of the normal-inverse-Wishart and inverse-Gamma distributions given { w m ′ } and {R m ′ }, with m ′ = m and c m ′ = c m . Finally, the integral can be analytically calculated as in Eq. ( 5), where the hyperparameters ψ 0 and the posterior hyper-parameters ψ ′ 0 are replaced by ψ 0c and ψ ′ 0c , respectively.

Symmetric vs. Asymmetric Learning

The MCMTL algorithm returns both the distribution over the weights for each task and the learned hierarchical model (task-class assignments). While the former can be used to evaluate the learning performance in the symmetric case, the latter provides a prior for learning a new task in the asymmetric scenario.

Symmetric Learning.

According to the generative model in Section 4.1, the task weights are distributed according to the normal distribution N (µ ′ 0m , Σ ′ 0m ), where µ ′ 0m and Σ ′ 0m are the posterior mean and covariance of the weight vector w m returned by the MCMTL algorithm. Since V m (x) = φ(x)

⊤ w m , the value of V m at a test state x * is distributed as

p V m (x * )|x * , µ ′ 0m , Σ ′ 0m = N φ(x * ) ⊤ µ ′ 0m , φ(x * ) ⊤ Σ ′ 0 -1 m φ(x * ) .
If MCMTL successfully clusters the task, we expect the value function prediction to be more accurate than learning each task independently.

Asymmetric Learning. In the asymmetric setting the class of the new task is not known in advance.

The inference problem is formalized as p w M+1 |R M+1 , ψ 0 , {c m } M m=1 , where w M+1 and R M+1 are the weight vector and rewards of the new task M M+1 , respectively. Similar to Section 4.2, this inference problem cannot be solved in closed form, thus, we must apply the MCMTL algorithm to the new task. The main difference with the symmetric learning is that the class assignments {c m } and weights { w m } for all the previous tasks are kept fixed, and are used as a prior over the new task learned by the MCMTL algorithm. As a result, the Gibbs sampling reduces to a one-step sampling process assigning the new task either to one of the existing classes or to a new class. If M M+1 belongs to a new class, then the inference problem becomes p (w M+1 |R M+1 , ψ 0 ), that is exactly the same as in STL. On the other hand, if M M+1 belongs to class c, the rewards and weights {R m ′ }, {w m ′ } of the tasks in class c can be used to compute the posterior hyper-parameters ψ ′ 0c as in Eq. ( 4), and to solve the inference problem p(w M+1 |R M+1 , ψ ′ 0c ).

Related Work

In RL, the approach of this paper is mainly related to [START_REF] Wilson | Multi-task reinforcement learning: A hierarchical Bayesian approach[END_REF]. Although we both use a DP-based HBM to model the distribution over the common structure of the tasks, in [START_REF] Wilson | Multi-task reinforcement learning: A hierarchical Bayesian approach[END_REF] the tasks share structure in their dynamics and reward function, while we consider the case that the similarity is in the value function. There are scenarios in which significantly different MDPs and policies may lead to very similar value functions. In such scenarios, the method proposed in this paper would still be able to leverage on the commonality of the value functions, thus performing better than single-task learning.

Moreover in [START_REF] Wilson | Multi-task reinforcement learning: A hierarchical Bayesian approach[END_REF], the setting is incremental, i.e., the tasks are observed as a sequence, and there is no restriction on the number of samples generated by each task. The focus is not on joint learning with finite number of samples, it is on using the information gained from the previous tasks to facilitate learning in a new one. This setting is similar to the asymmetric learning considered in our work.

In supervised learning, our work is related to [START_REF] Yu | Learning Gaussian processes from multiple tasks[END_REF] and [START_REF] Xue | Multi-task learning for classication with dirichlet process priors[END_REF]. In [START_REF] Yu | Learning Gaussian processes from multiple tasks[END_REF], the authors present a single-class HBM for learning multiple related functions using GPs. Our single-class model of Section 3 is an adaptation of this work for RL using GPTD. Besides, our multi-class model of Section 4 extends this method to the case with an undefined number of classes. In [START_REF] Xue | Multi-task learning for classication with dirichlet process priors[END_REF], a DP-based HBM is used to learn the extent of similarity between classification problems. The problem considered in our paper is regression, the multi-class model of Section 4 is more complex than the one used in [START_REF] Xue | Multi-task learning for classication with dirichlet process priors[END_REF], and the inference algorithms of Section 4 are based on Gibbs sampling, where a variational method is used for inference in [START_REF] Xue | Multi-task learning for classication with dirichlet process priors[END_REF].

Experiments

In this section, we report empirical results applying the Bayesian multi-task learning (BMTL) algorithms presented in this paper to a regression problem and a benchmark RL problem, inverted pendulum. We compare the performance of single-task learning (STL) with single-class multi-task learning (SCMTL), i.e., all tasks are assumed to belong to the same class, and multi-class multi-task learning (MCMTL), i.e., tasks belong to a number of classes not known in advance. By STL, we refer to running the EM algorithm of Section 3.2 for each task separately. The reason to use the regression problem in our experiments is that it allows us to evaluate our BMTL algorithms when the tasks are generated exactly according to the generative models of Sections 3 and 4.

A Regression Problem

In this problem, tasks are functions in the linear space spanned by a feature space φ(x) = 1, x, x 2 , x 3 , x 4 , x 5 ) ⊤ on the domain X = [-1, 1]. The weights for the tasks are drawn from four different classes, i.e., four 6-dim multivariate Gaussian distributions, with the parameters shown in Fig. 3(a). The noise covariance matrix S = diag(σ 2 ) for all the algorithms. We evaluate the performance of each BMTL algorithm by computing its relative mean squared error (MSE) improvement over STL : (M SESTL -M SEBMTL )/M SESTL. The MSEs are computed over N ′ = 1000 test samples. All the reported results are averaged over 200 independent runs. Min Max c1 37.35% ± 1.2 64.78% ± 0.9 c2 29.86% ± 0.9 67.64% ± 0.8 c3 66.69% ± 0.8 91.81% ± 0.6 c4 56.47% ± 1.2 75.78% ± 0.3 In the first experiment, we draw all the tasks from class c 2 . Fig. 3(c) shows the performance of SCMTL for different number of tasks (M ) and samples per task (N ). SCMTL achieves an improvement over STL that varies from 29.86% ± 0.9 for N = 100 and 20 tasks to 67.64% ± 0.8 for 100 tasks with only 20 samples each. The results indicate that SCMTL successfully takes advantage of the samples coming from all the tasks to build a more accurate prior than the one obtained by considering each task separately as in STL. However, the advantage of SCMTL over STL declines as N is increased. In fact, as STL converges, SCMTL cannot make further improvement. We repeated the same experiment for the other classes. The minimum and maximum performance of SCMTL for all the classes (all obtained for N = 100, M = 20 and N = 20, M = 100, respectively) are summarized in Fig. 3(b).

Number of Samples

Number of Tasks

In the second experiment, we draw the tasks randomly from the four classes. We first apply SCMTL to this problem. Fig. 4(a) shows the SCMTL's performance. As it can be seen, the results are worse than those in the first experiment (Fig. 3(c)), varying from 30.15% ± 4.8 to 54.05% ± 1.2. By clustering all the tasks together, SCMTL takes advantage of all the available samples, thus, performs better than STL. However, when the tasks are drawn from significantly different distributions, it learns a very general prior which does not allow a significant improvement over STL. We then apply MCMTL to this problem. MCMTL's performance (Fig. 4(b)) is significantly better than SCMTL's (Fig. 4(a)), and it varies from 45.64% ± 5.6 to 77.65% ± 0.8. In order to evaluate how well MCMTL classifies the tasks, we also compare its performance to a version of MCMTL in which each task is assigned to the right class in advance. The difference between the two algorithms is statistically significant only for N = 20 (with the maximum of 5.08% for M = 20), in which the noise on the samples makes it more difficult to discriminate between the distributions generating the tasks, and thus, to classify them correctly.

Finally, we compare the performance of SCMTL and MCMTL in the asymmetric setting. At the end of each run of the symmetric problem, we draw 100 new test tasks at random from the same four classes used to generate the training tasks. We run the asymmetric algorithm described in Section 4.3 on each of the test tasks separately. Fig. 4(c) shows the performance of SCMTL and MCMTL for different number of training tasks and N fixed to 20. The results indicate that MCMTL performs relatively better than SCMTL as the number of training tasks increases.

Inverted Pendulum

The experiments of Section 6.1 indicate that when the tasks are generated exactly according to the generative models of Sections 3 and 4, the BMTL methods can significantly improve the performance of a regression problem w.r.t. STL. As discussed in Section 2, the policy evaluation step of policy iteration can be casted as a regression problem, thus, similar improvement can be expected.

In this section, we compare our BMTL algorithms with STL in the problem of learning a control policy for balancing an inverted pendulum. Dynamics, reward function, and basis functions are the same as in [START_REF] Lagoudakis | Least-squares policy iteration[END_REF]. Each task is generated by drawing the parameters of the dynamics (pole mass, pole length, cart mass, and noise on the actions) from Gaussian distributions with means and variances summarized in Fig. 5(a). The distribution over the two classes is uniform. It is worth noting that, unlike the regression experiments, here we have no guarantee that the weights of the value functions will follow the generative models assumed by the BMTL methods. We use policy iteration with 10 iterations and the noise correlation matrix P -1 = Φm ⊤ (ΦmΦm ⊤ ) -1 Φm for all the algorithms. 1 In STL, each policy evaluation step is solved using the EM algorithm of Section 3.2 for each task separately, where in BMTL, it is solved by running SCMTL or MCMTL over all the tasks. All the results are averaged over 150 independent runs. Fig. 5(b) shows the performance of the policy learned by STL, SCMTL, and MCMTL for M = 10 tasks and different (up to 500) number of samples per task. Note that STL converges at about 1200 samples per task with an average performance of 2473 ± 61.9 balanced steps. As it can be seen, both BMTL methods outperform STL, and MCMTL achieves a better performance than SCMTL as the number of samples is increased. Since SCMTL forces all the tasks to have weights generated from a common distribution, it learns a very general prior, and thus, it cannot approximate the value functions as accurate as MCMTL, which is able to correctly discriminate between the tasks in classes c 1 and c 2 . In order to show how the performance changes with different number of tasks, we compute the area ratio [START_REF] Taylor | Transfer learning via inter-task mappings for temporal difference learning[END_REF] on the first 500 samples as ρ BMTL = ABMTL-ASTL ASTL , where A STL (A BMTL ) is the area under the learning curve of STL (BMTL) from 100 to 500 samples. Fig. 5(c) shows that MCMTL has significantly better area ratio than SCMTL for all values of M except very small ones.

Conclusions

We presented hierarchical Bayesian models (HBMs) and inference algorithms for multi-task reinforcement learning (RL) where the tasks share structure in their value functions. To the best of our knowledge, this is the first work that models value function similarity using HBMs. In particular, we considered two cases, where all the value functions belong to the same class, and where they belong to an undefined number of classes. In these cases, we modelled the distribution over the value functions using a parametric HBM and a Dirichlet process (DP) based non-parametric HBM, respectively. For each case, we derived inference algorithms for learning the value functions jointly and to transfer the knowledge acquired in the joint learning to improve the performance of learning the value function of a new task.

We first applied our proposed Bayesian multi-task learning (BMTL) algorithms to a regression problem, in which the tasks are drawn from the generative models used by the BMTL methods. The results indicate that BMTL algorithms achieve significant improvement over single-task learning (STL) in both symmetric and asymmetric settings. We then applied our BMTL algorithms to a benchmark RL problem, inverted pendulum. Although the tasks are no longer generated according to the models used by the BMTL algorithms, they still outperform STL. In our DP-based model we used Gibbs sampling, the most common simulation tool for Bayesian inference. We plan to look into variational techniques for Bayesian inference as an alternative approach.

Appendix

Details of the Single-class Multi-task Model

• R m p(R m |w m , θ) = p(R m |w m , σ 2 ) = N (HΦ ⊤ mn w m , σ 2 P) • w m p(w m |θ) = p(w m |µ, Σ) = N (µ, Σ) = (2π) -d/2 |Σ| -1/2 exp - 1 2 (w m -µ) ⊤ Σ -1 (w m -µ)
• θ p(θ|ψ 0 ) = p(µ, Σ; µ 0 , k 0 , ν 0 , Σ 0 ) × p(σ 2 ; α 0 , β 0 ) = N (µ; µ 0 , Σ/k 0 ) IW(Σ; ν 0 , Σ 0 ) × IG(σ 2 ; α 0 , β 0 )

= (2π/k 0 ) -d/2 |Σ| -1/2 exp - k 0 2 (µ -µ 0 ) ⊤ Σ -1 (µ -µ 0 ) (normal) × B|Σ 0 | ν0/2 |Σ| -(ν0+d+1)/2 exp - 1 2 tr Σ 0 Σ -1 (inverse-Wishart) × β α0 0 Γ(α 0 ) 1 σ 2 α0+1 exp - β 0 σ 2 (inverse-Gamma)
where

B -1 = 2 ν0d/2 π d(d-1)/4 d j=1 Γ ν 0 + 1 -j 2 .
p(µ, Σ, σ 2 |w m , R m , ψ 0 ) = N (µ; µ ′ 0 , Σ/k ′ 0 ) IW(Σ; ν ′ 0 , Σ ′ 0 ) × IG(σ 2 ; α ′ 0 , β ′ 0 ) (6)
where the posterior hyper-parameters

ψ ′ 0 = (µ ′ 0 , k ′ 0 , ν ′ 0 , Σ ′ 0 , α ′ 0 , β ′ 0 ) are µ ′ 0 = M k 0 + M w + k 0 µ 0 k 0 + M , (7) k
′ 0 = k 0 + M , (8) ν
′ 0 = ν 0 + M , (9)
Σ ′ 0 = Σ 0 + Q 0 + k 0 M k 0 + M -m,c ( w -µ 0 )( w -µ 0 ) ⊤ , (10)
α ′ 0 = α 0 + N M 2 , (11)
β ′ 0 = β 0 + 1 2 M m=1 (HΦ ⊤ m w m -R m ) ⊤ P -1 (HΦ ⊤ m w m -R m ) , (12)
where

w = 1 M M m=1 w m , Q 0 = M m=1 (w m -w)(w m -w) ⊤ .

Gibbs Sampling

In this section, we report the equations used in the Gibbs sampling described in Section 4.2 of the paper. At each iteration of the MCMTL algorithm of Figure 2 in the paper, the Gibbs sampling is fed with observations ({ w m }, {R m }), where w m = µ ′ 0m . In particular, we use the the Gibbs sampling with conjugate-priors (Algorithm 3 in [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF]).

We begin with the probability of task m belonging to a new class. Given observation ( w m , R m ) and hyper-prior ψ 0 , the non-normalized probability can be written as

p(R m , w m |θ)p(θ|ψ 0 )dθ = k 0 πk ′ 0 d/2 |Σ 0 | ν0/2 |Σ ′ 0 | ν ′ 0 /2 Γ ν ′ 0 2 Γ ν ′ 0 -d 2 × (2π|P |) -N/2 β α0 0 Γ(α 0 ) Γ(α ′ 0 ) β ′α ′ 0 0
, where hyper-parameters ψ ′ 0 = (µ ′ 0 , k ′ 0 , ν ′ 0 , Σ ′ 0 , α ′ 0 , β ′ 0 ) are computed as in Section 8.2 for observation (R m , w m ).

Similarly, let m ′ : c m ′ = c, m ′ = m, then the posterior over parameters for class c is

p(µ c , Σ c , σ 2 c |{ w m ′ }, {R m ′ }, ψ 0 ) = N (µ c ; µ 0c , Σ c /k 0 )IW(Σ c ; ν 0c , Σ 0c ) × IG(σ 2 c ; α 0c , β 0c ) (13)
with the following posterior hyper-parameters (which play the role of prior hyper-parameters for the class)

µ 0c = M -m,c k 0 + M -m,c w + k 0 µ 0 k 0 + M -m,c , (14)
k 0c = k 0 + M -m,c , (15)
ν 0c = ν 0 + M -m,c , (16)
Σ 0c = Σ 0 + Q 0 + k 0 M -m,c k 0 + M -m,c ( w -µ 0 )( w -µ 0 ) ⊤ , (17)
α 0c = α 0 + N M -m,c 2 , (18)
β 0c = β 0 + 1 2 m ′ (HΦ ⊤ m w m -R m ) ⊤ P -1 (HΦ ⊤ m w m -R m ) , (19)
where M -m,c is the number of tasks belonging to class c except task m, w = 1 M-m,c m ′ w m ′ , and

Q 0 = m ′ ( w m ′ -w)( w m ′ -w)
⊤ . As a result, the integral for the probability of task m belong to class c becomes

p R m , w m |θ p θ|{c m ′ }, ψ 0 dθ = p R m , w m |θ p θ|ψ 0c dθ = k 0c πk ′ 0c d/2 |Σ 0c | ν0c/2 |Σ ′ 0c | ν ′ 0c /2 Γ ν ′ 0c 2 Γ ν ′ 0c -d 2 × (2π|P |) -N/2 β α0c 0c Γ(α 0c ) Γ(α ′ 0c ) β ′α ′ 0c 0 ,
where ψ ′ 0c is computed as in Equations ( 2)-( 7) but using ψ 0c as prior instead of ψ 0 .

Noise Correlation Models

We analyze the equations for w m for two different noise correlation models. In particular we show that depending on the covariance matrix, both GPTD and BMTL can be seen as extensions of either Bellman residual minimization (BRM) or LSTD.

We first analyze the general formulation with S = σ 2 P be the covariance matrix of the noise.

Let θ = (µ, Σ, σ 2 ) be the model parameters. The expected value of the weights in BMTL is written as

w m = 1 σ 2 Φ m H ⊤ P -1 HΦ ⊤ m + Σ -1 -1 1 σ 2 Φ m H ⊤ P -1 R m + Σ -1 µ , (20)
where S = σ 2 P is the noise covariance matrix modeling the correlation of the noise at different states. We call σ 2 and P the noise variance and the noise correlation matrix, respectively.

In case µ = 0 and Σ = I, we obtain a general form for the posterior mean of the weights in the parametric form of GPTD (Equation (4.4.41) in [START_REF] Engel | Algorithms and Representations for Reinforcement Learning[END_REF])

w m = Φ m H ⊤ P -1 HΦ ⊤ m + σ 2 I -1 Φ m H ⊤ P -1 R m . (21)
The (parametric) GPTD of [5; 6] is obtained by setting P = HH ⊤ in Equation (21). As it was shown in Section 4.5 of [START_REF] Engel | Algorithms and Representations for Reinforcement Learning[END_REF], by setting σ 2 → 0 and P -1 = Φ m ⊤ GΦ m in Equation ( 21), where G is an arbitrary d × d symmetric positive-definite matrix, we can derive a new set of GPTD algorithms that are based on the LSTD(0) algorithm. As it was discussed in Section 4.5.3 of [START_REF] Engel | Algorithms and Representations for Reinforcement Learning[END_REF]

, a reasonable choice for G is G = Φ m Φ m ⊤ -1
.

Experiment Setups

In the following we list all the details about the setup used in the experiments of the paper. We report the hyper-prior parameters (τ, ψ 0 ) and the parameters used in the inference algorithm. In particular, ǫ EM is the threshold used in the stopping condition of the EM algorithm, n MCMTL is the maximum number of iterations of the outer loop of MCMTL, and n Gibbs is the number of steps of the Gibbs sampling. In none of the experiments the parameters have been systematically optimized.

The Regression Problem

In Figure 6 we report examples of the functions (tasks) and the generated samples for each of the four different classes used in the experiments of Section 6.1 of the paper. The parameters used in the experiments are reported in Tables 1. As it can be noticed the prior is not very informative and it has not been optimized for this specific problem. Since the four classes are quite well separated, the value of the concentration τ is not critical for the success of the MCMTL algorithm. Finally, the length of the Gibbs sampling changes with the number of tasks, so that when many tasks are involved, a longer MCMC simulation is performed and a more accurate estimation of p({c m }|{ w m }, {R m }, τ, ψ 0 ) is computed. Table 1: Hyper-prior parameters and parameters of the inference algorithm in the regression problem.

Inverted Pendulum

In the inverted pendulum the state space X ∈ R 2 is a 2-dimensional space with variables (θ, θ), the vertical angle and the angular velocity, respectively. The action space is A = {-50, 0, 50}, where a = 50 means that a right force of 50 Newtons is applied to the cart. Each action is perturbed by a uniform noise in [-η, η] (as reported in the paper for each task η is drawn from a Gaussian distribution N (12, 0.1)). The discount factor is γ = 0.9. At each step of policy iteration, Q(s, a) is approximated in a linear space spanned by 10 RBF features for each action as in [START_REF] Lagoudakis | Least-squares policy iteration[END_REF]. The critical parameters controlling the nonlinear dynamics of the system are m the mass of the pendulum, M the mass of the cart, and l the length of the pendulum. To illustrate the impact of these parameters on the optimal value function and the corresponding optimal policy, we show the value functions and the policies for two sample tasks drawn from classes c 1 and c 2 in Figure 7. As it can be noticed different values for m, M , and l, induce significantly different value functions and policies. The parameters used in the experiments are reported in Tables 2. m have the same x n , a n , x ′ n components, while the fourth component is changed according to π k-1 (x ′ n ). At the end of policy iteration, the performance of the learned policy is evaluated by taking the average over 30 episodes with a maximum of 3,000 steps each.

Figure 1 :

 1 Figure 1: Graphical representations for (a) the single-task model -an extension of GPTD by defining a hyper-prior over the parameters, (b) the single-class multi-task model of Section 3, and (c) the multi-class multi-task model of Section 4.

MCMTL

  {Rm}, τ, ψ0 Initialize {cm} repeat for c = 1, . . . , C do Initialize θc repeat for m : cm = c do p(wm|Rm, θc) = N (µ ′ 0m , Σ ′ 0m ) (E-step) end for Optimize θc (M-step) until convergence end for Set wm = µ ′ 0m , m = 1, . . . , M p({cm}|{ wm}, {Rm}, τ, ψ0) until convergence return { wm} and {cm}

Figure 2 :

 2 Figure 2: The inference algorithm for the multi-class multi-task learning (MCMTL) scenario.

Figure 3 :

 3 Figure 3: (a) class parameters, (b) minimum and maximum improvement of SCMTL over STL in each class, (c) relative MSE improvement of SCMTL over STL when all the tasks are drawn from class c 2 .

Figure 4 :

 4 Figure 4: Results for the case that the tasks are drawn randomly from the four classes: (a) relative MSE improvement of SCMTL over STL, (b) relative MSE improvement of MCMTL over STL, (c) asymmetric performance of MCMTL and SCMTL for N = 20.

Figure 5 :

 5 Figure 5: Results for the inverted pendulum problem: (a) distributions of the parameters of the dynamics, (b) comparing the performance of STL, SCMTL, and MCMTL in terms of the number of balanced steps for M = 10, (c) comparing the performance of SCMTL and MCMTL in terms of the area ratio on the first 500 samples.

8. 2

 2 Posterior Distribution of the Parameters with the Normal-Inverse-Wishart × Inverse-Gamma Prior Taking advantage of the conjugate prior the posterior distribution over θ = (µ, Σ, σ 2 ) given observations {(w m , R m )} M m=1 and hyperprior ψ 0 is

Figure 6 :

 6 Figure 6: Examples of functions and samples drawn from each of the 4 classes (N = 100, M = 4).

Figure 7 :Table 2 :

 72 Figure7: Examples of value functions and corresponding policies for two tasks with parameters m 1 = 1.03, M 1 = 5.64, l 1 = 0.5, η 1 = 12.09 and m 2 = 3.0, M 2 = 10.04, l 2 = 2.01, η 2 = 11.91, drawn from classes c 1 and c 2 , respectively. The tasks are solved using STL with N = 2, 500. In the policy plots blue, green, and red colors correspond to actions -50, 0, and 50, respectively.

This is the noise correlation matrix of LSTD(0) in the parametric GPTD form (see Section 8.4).


J.-M Billiot

F Corset

E Fontenas

Continuum Percolation in the Relative Neighborhood Graph

In the present study, we establish the existence of nontrivial site percolation threshold in the Relative Neighborhood Graph (RNG) for Poisson stationary point process with unit intensity in the plane.

Introduction

Percolation theory is very useful to describe various physical phenomena. In particular, there are important connections with phase transition problems [START_REF] Grimmet | Percolation, Second Edition[END_REF][START_REF] Georgii | The random geometry of equilibrium phases[END_REF].

The interest for percolation problems has grown rapidly during the last decades: see Lyons and Peres [START_REF] Lyons | Probability on Trees and Networks[END_REF] for percolation on trees and networks and Meester and Roy [START_REF] Meester | Continuum Percolation[END_REF] for continuum percolation and the references therein. In 1996, Häggström and Meester [START_REF] Häggström | Nearest Neighbor and Hard Sphere Models in Continuum Percolation[END_REF] proposed results for continuum percolation problems for the k-nearest neighbor graph under Poisson process. In a recent paper, Balister and Bollobás [START_REF] Balister | Percolation in the k-nearest neighbor graph[END_REF] give bounds on k for the k-nearest neighbor graph for percolation with several possible definitions.

Site and bond percolation as well as phase transition for several statistical mechanics models have also been studied by Häggström [START_REF] Häggström | Markov random fields and percolation on general graphs[END_REF] for general graphs and in particular for the Delaunay triangulation for the Poisson point process. It also shows for graphs with a bounded degree if the bond percolation threshold is not trivial then it implies a phase transition for the Ising model on this graph. For general results on Delaunay graphs and Voronoi tesselations see Moller [START_REF] Møller | Lectures on Random Voronoi Tessellations[END_REF]. For recent results on percolation on these graphs see Balister et al. [START_REF] Balister | Percolation in a voronoi process[END_REF][START_REF] Balister | Bond percolation with attenuation in high dimensional voronoi tilings[END_REF].

Benjamini and Schramm [START_REF] Benjamini | Percolation beyong Z d , many questions and a few answers[END_REF] proposed a comprehensive study on general graphs, with special focus on Cayley graphs, quasi-transitive graphs and planar graphs. Recently, Procacci and Scoppola [START_REF] Procacci | Infinite graphs with a nontrivial bond percolation threshold: some sufficient conditions[END_REF] proposed sufficient conditions on infinite graphs to deduce a non trivial bond percolation threshold. Among these assumptions, they assume that the dual graph is bounded degree. It is interesting to relax this condition in order to deal with proximity random graphs (which in general have a dual not bounded degree) like the skeletons on point Poisson processes of the plane. These random graphs includes the Gabriel graph and the Relative Neighborhood Graph (RNG) which are important for many applications.

Remark that continuous models defined on nearest neighbors graphs are interesting for small temperature as alternative of standard models on regular networks, because it allows vibrations and deformations of the network and may be find an application in physics of the solid state. In particular one example is the study of the order-disorder transition of binary alloys or ionic cristals. It is well-known that Delaunay graph or Voronoi regions (rather called Wigner-Seitz grid and Brillouin zone in physics framework) take a fondamental place for the understanding of the electrical current, waves propagation and phase transitions observed by Bragg diffraction of X rays.

Another domain of application should be found in cancerology for the study of the growth of tumour when the cancer cells suddently begin to invade healthy tissue. The Delaunay graph is well adapted for such study as explained in [START_REF] Bertin | Spatial Delaunay Gibbs Point Processes[END_REF] but the RNG should give some more information. More precisely, in histology we have some slides with marked cells (cancer cells and normal cells in first approximation) and probably, a connection between percolation and the grade of a given cancer will be helpful to give an aid for the diagnostic of a pathologist.

The existence and unicity at small activity of nearest neighbors stationary Gibbs states can be found in Bertin et al. [START_REF] Bertin | Existence of Delaunay Pairwise Gibbs Point Processes with Superstable Component[END_REF][START_REF] Bertin | Existence of "Nearest-Neighbour" Gibbs Point Models[END_REF][START_REF] Bertin | k-Nearest-Neighbour Gibbs Point Processes[END_REF]. In Bertin et al. [START_REF] Bertin | Phase Transition in Nearest-Neighbour Continuum Potts Models[END_REF], the phase transition in the Delaunay continuum Potts model is established. It is a generalization of the Lebowitz-Lieb model as described in Georgii and Häggström [START_REF] Georgii | Phase transition in continuum Potts models[END_REF] where the soft repulsion between several species of particules acts on the Delaunay graph. What is the good Delaunay subgraph on which the repulsion is strong enough to maintain a phase transition? In terms of percolation, it means: is bond percolation maintained in this subgraph? Bertin et al. [START_REF] Bertin | Continuum Percolation in the Gabriel Graph[END_REF] gave an answer for the Gabriel graph. Another well-known subgraph of Delaunay graph is the Minimum Spanning Tree (MST). The structure of spanning forest generalized minimum spanning tree for infinite graphs. Constructed with the greedy algorithm, this graph exactly is described in the book of Meester and Roy [START_REF] Meester | Continuum Percolation[END_REF]. The link between branching number and percolation on trees is proved by Lyons [START_REF] Lyons | Random walks and percolation on trees[END_REF]. Connection between minimum spanning forest and occupied and vacant percolation is strong (for more details on simultaneous uniqueness see Alexander [START_REF] Alexander | Percolation and minimal spanning forests in infinite graphs[END_REF]). Thus, the minimum spanning forest is a tree with one infinite path a.s. in two dimensions when the points are distributed under a stationary Poisson process. It comes that a.s. the site or bond critical threshold are equal to 1.

The present study gives an answer for the RNG, well known in computational geometry, when the points are distributed under a stationary Poisson point process with unit intensity in the plane. We adapted a powerful method of the rolling ball proposed by Balister and Bollobás [START_REF] Balister | Percolation in the k-nearest neighbor graph[END_REF] relying on 1-independent bond percolation on Z 2 . Then if we control the probability of having less than a fixed number of points in a given region and considering the event that all the sites are open in this region, we can proceed similarly as in Häggström [START_REF] Häggström | Markov random fields and percolation on general graphs[END_REF].

The paper is organized as follows. The first section is devoted to some definitions and notations. Next, the main result on site and bond percolation on the RNG is given. Then, we give the proof of the main result, by using a result of bond percolation in the 1-independent case in Z 2 and the rolling ball method. We conclude on possible extensions of this work.

Notations and definitions

Let |A| denotes the Lebesgue measure if the set A is a bounded Borel set of R 2 , and the counting measure if A is a discrete set.

Given a finite box Λ ⊂ R 2 , we denote by Π Λ the Poisson point process on the locally finite set of points in Λ denoted by Ω Λ with intensity 1 i.e.

ΩΛ f dΠ Λ = exp (-|Λ|) ∞ n=0 1 n! Λ n f ({x 1 , . . . x n })dx 1 , . . . dx n
for any bounded measurable function f on Ω Λ .

Let Ω the set of locally finite subsets of R 2 . We have to consider only the configurations ϕ ∈ Ω which are in general position (four points on the same circle do not occur and no three points are colinear) in order to ensure the existence and unicity of the Delaunay graph. One can notice that, for any stationary point processes, the probability of the set of tesselations in general position is equal to one [START_REF] Møller | Lectures on Random Voronoi Tessellations[END_REF]. Let us recall the definitions of Delaunay, Gabriel and relative neighborhood graphs.

Definition 3

The Relative Neighborhood graph Rn(ϕ) of a configuration ϕ in Ω is defined as the set of edges {u, v} ⊂ ϕ such that the intersection of the disks with center u and v with radius uv does not contain any point of the configuration ϕ.

These graphs are planar in R 2 . Furthermore, the RNG is a subgraph of the Gabriel graph which is a subgraph of the Delaunay Graph (see figure 1).

Main result

We first deal with the site percolation on the RNG for the Poisson point process Φ. We introduce the Bernoulli process Θ(Φ) providing the type picking mechanism (1 for open and 0 for closed) of the points (or sites) in Φ. Obviously, Φ = (Φ, Θ(Φ)) can be seen as a marked Poisson process. The probability measure of Φ is given by : P(dϕ, dθ) = Π(dϕ)µ p ϕ (dθ) We now introduce p site c (Rn, Π) defined as the lowest p for which the probability of the event that there exists an infinite open cluster in the RNG relative to the marked point Poisson process Φ is equal to 1. By ergodicity of Π, this previous event, invariant by translation, is a trivial event. The marked point Poisson process Φ may exhibit some percolation phenomenon with critical value p site c (Rn, Π). We want to prove in the following theorem that this is a non trivial critical value, i.e., The following section is devoted to the proof of this theorem.

Proof

First, we point out that the method proposed by Häggström [START_REF] Häggström | Markov random fields and percolation on general graphs[END_REF] and used for k nearest neighbor graph and Delaunay graph, as well as the adaptation by Bertin et al. [START_REF] Bertin | Continuum Percolation in the Gabriel Graph[END_REF] to the Gabriel graph does not applied to the RNG. However, it is not sufficient to control the probability of having at least one point and less than a fixed number of points in each small box K (see figure 2).

Indeed, we may choose some configuration of points such that the length and number of points of a Rnpath between two points are arbitrary large, see figure 2. wt is an edge of the RNG but uv is not because for example the point q belongs to the vacuity region of this edge. Moreover, such methods are based on comparison with independent bond percolation on the grid Z 2 . As, in Balister and Bollobás [START_REF] Balister | Percolation in the k-nearest neighbor graph[END_REF] in the case of k nearest neighbor graph, we adapt the method of rolling ball (see figure 3) to the RNG. We procceed in two steps:

1. Controling the probability of some suitable configurations of points under Poisson point process.

2. Consider the Bernoulli site percolation in such configurations.

To prove that continuous percolation occurs, we shall compare the process to various bond percolation models on Z 2 . In these models, the states of the edges will not be independent. However they will satisfy the following definition: Definition 4 A bond percolation model is 1-independent if whenever E 1 and E 2 are sets of edges at graph distance at least 1 from each another (i.e., if no edge of E 1 is incident to an edge of E 2 ) then the state of the edges in E 1 is independent of the state of the edges in E 2 .

We shall use the following result in Balister et al. [START_REF] Balister | Continuum percolation in the square and the disk[END_REF]. Let us first consider the case of percolation in the RNG. Write u ∼ v if uv is an edge of the underlying graph Rn(ϕ). For percolation we need to find an infinite path, i.e., a sequence u 1 , u 2 . . . with u i ∼ u i+1 for all i. Consider the rectangular region consisting of two adjacent squares S 1 , S 2 shown in figure 3. Both S 1 and S 2 have side length 2r + 2s, where r and s are to be chosen later. We define the basic good event E S1,S2 to be the event that every vertex u 1 in the central disk C 1 of S 1 is joined to at least one vertex v in the central disk C 2 of S 2 by a Rnpath, regardless of the state of the Poisson process outside of S 1 and S 2 . Now consider the following percolation model on Z 2 . Each vertex (i, j) ∈ Z 2 corresponds to a square [Ri, R(i + 1)] × [Rj, R(j + 1)] ∈ R 2 , where R = 2r + 2s, and an edge is open between adjacent vertices (corresponding to squares S 1 and S 2 ) if both the corresponding basic good events E S1,S2 and E S2,S1 hold. Note that this is indeed a 1-independent model on Z 2 since the event E S1,S2 depends only on the Poisson process within the region S 1 and S 2 , and thus sets of edges at distance at least one apart in Z 2 depend on the Poisson process in disjoint regions of R 2 . Any open path in Z 2 corresponds to a sequence of basic good events E S1,S2 , E S2,S3 . . . that occur, where S i is the square associated with a site in Z 2 . Every vertex u 1 of the original Poisson process that lies in the central disk C 1 of S 1 now has an infinite path leading away from it, since one can find points u i in the central disk of S i and paths from u i-1 to u i inductively for every i > 1. In particular, each such u 1 lies in an infinite component. Moreover, such vertices exist in C 1 , so there is an infinite component. One can choose r and s so that the probability that the intersection of a basic good event is large and then we will apply the theorem 2.

In order to bound the probability that this intersection of a basic good event fails, we shall use the following rolling ball method. Let C 1 , C 2 , and L be as in Figure 3. (L is the region between the two disks C 1 and C 2 .) We need to define E S1,S2 the event that for every point v ∈ C 1 ∪ L, there is a u such that: a) v ∼ u; b) d(u, v) ≤ s; and c) u ∈ D v , where D v is the disk of radius r inside C 1 ∪ L ∪ C 2 with v on its C 1side boundary (the dotted disk in Figure 3). Note in particular that (b) implies that the condition u ∼ v in (a) is independent of the Poisson process outside of S 1 ∪ S 2 . This is because both u and v are at distance at least s from the exterior of S 1 ∪ S 2 , so the event that (u, v) is an edge of the RNG only depends on the points within S 1 ∪ S 2 . We denote ĒS1,S2 the complementary of E S1,S2 . The probability of ĒS1,S2 is bounded by the expected number of points u for which above conditions (a)-(c) fail. Thus, we have

Π( ĒS1,S2 ) ≤ 2r(2r + 2s)p Rn,r,s (1)
where p Rn,r,s is the probability that (a)-(c) fail for some fixed v. Notice that this probability is independent of the location of v in C 1 ∪ L.

Lemma 1

We can choose r and s such that Π( ĒS1,S2 ) is arbitrary small.

Proof:

Let D(x, α) be the disk of radius α and of center x. To calculate this last integration, we have :

             arccos( α
         |L(α, r, θ)| = α 2 θ + (α 2 -r 2 ) arcsin α sin θ √ r 2 + α 2 -2αr cos θ + αr sin θ |L(α, r, arccos α 2r -π/3)| = α 2 2 ( π 3 - √ 3
2 ) -r 2 arcsin α 2r + rα For instance, the bound involved in inequality 1 gives around 10 -40 with r = s = 8000.

For our purpose, we also need to control the probability for the Poisson point process of having at least one point in C 1 and less than m points in C 2 ∪ C 1 ∪ L.

We denote

F C1 = {Φ(C 1 ) ≥ 1} and A m = {Φ(C 2 ∪ C 1 ∪ L) ≤ m}. Notice that as E S1,S2 ∩ F C1 ∩ A m ⊂ E S1,S2 then E S1,S2 ∩ E S2,S1 ∩ F C1 ∩ F C2 ∩ A m ⊂ E S1,S2 ∩ E S2,S1 .
We obtain that This bound becomes negligible whenever r ≫ 1 and m > e(8 + π)r 2 using Stirling formula.

Π(E S1,S2 ∩ E S2,S1 ) ≥ 1 -Π( ĒS1,S2 ) + Π( ĒS2,S1 ) + Π( FC1 ) + Π( FC2 ) + Π( Ām ) where Π( FC1 ) = Π( FC2 ) = e -
Set ǫ = 0.1361. Similarly as in Häggström [START_REF] Häggström | Markov random fields and percolation on general graphs[END_REF] but adapted in the 1-independent case, let B r,s be the event that all the sites are open in C 2 ∪ C 1 ∪ L with probability p = 1 -ǫ 2m and also define C r,s = E S1,S2 ∩ E S2,S1 ∩ B r,s . Then, we have:

P (C r,s ) =
Π(E S1,S2 ∩ E S2,S1 ∩ F C1 ∩ F C2 ∩ A m ) ≥ 1 -ǫ/2
We conclude with theorem 2 that p site c (Rn, Π) ≤ 1 -ǫ 2m < 1.

Concluding remarks

This kind of proof also apply in the case of bond or site percolation of the k nearest neighbor graph: it is direct consequence of [START_REF] Balister | Percolation in the k-nearest neighbor graph[END_REF]. It is sufficient to use what they called p U , p B , p I or p O for several possible definitions of percolation. As suggested in [START_REF] Balister | Percolation in the k-nearest neighbor graph[END_REF] with a high confidence, k = 3 is the critical out-degrees for percolation on the k nearest neighbor graph. We notice that in the RNG the number of neighbors is bounded by 6. So a point have neighbors in several directions. It may be interesting to study a family of Delaunay subgraphs defined on sequence of vacuity regions such that the number of neighbors of each point is the lowest as possible but keep good connectivity properties for percolation purposes. Otherwise a challenge would be to extend the method of the rolling ball when points are distributed under a Gibbs point process for instance a hard-core point process.

Definition 1 Definition 2

 12 The Delaunay graph Del 2 (ϕ) of a configuration ϕ in Ω is the set of edges of the unique triangulation Del 3 (ϕ) in which the interior of the circle circumscribed by every triangle of Del 3 (ϕ) does not contain any point of ϕ. The Gabriel graph Gab(ϕ) of a configuration ϕ in Ω is defined as the set of edges {u, v} ⊂ ϕ such that the open circle with {u, v} as diameter does not contain any point of the configuration ϕ.

Figure 1 :

 1 Figure 1: Vacuity regions for Gabriel and Relative neighborhood graphs for an edge uv

Theorem 1 p site c (

 c Rn, Π) < 1 and p bond c (Rn, Π) < 1.

Figure 2 :

 2 Figure 2: Example of construction of configuration with arbitrary large Rn-path.

Figure 3 : 2

 32 Figure 3: The Rolling Ball Method

p 2 4r 2 + (2r 2 -s 2 ) arcsin s 2r + s 2 π/ 2 .

 2222 Rn, r, s ≤ e -|Dv∩D(v,s)| + Dv ∩D(v,s) p r, s (u) du where p r, s (u) = e -|Dv∩D(v, d(u,v))| (1 -e -|D(v, d(u,v))∩D(u, d(u,v))\Dv | ) is the probability that u is the closest point to v inside D v , but that (u, v) is not an edge of the RNG. To calculate this upper bound, note that |D v ∩ D(v, s)| = -rs 1 -s By choosing polar coordinates (α, θ) of u, it comes: Dv∩D(v,s) p r, s (u) du = 2 s 0 αe -|Dv ∩D(v,α)| arccos(α/(2r)) 0 J(α, θ) dθ dα where J(α, θ) = 1 -e -|D(v, α)∩D((α,θ), α)\Dv | .

1 -e -α 2 Figure 4 :

 124 Figure 4: The Lune L(α, r, θ) = D(u, α)\D v is the dotted area

2 1 -α 2 4r 2 .

 22 Thus,Dv ∩D(v,s) p r, s (u) du = 2s 0 αe -|Dv ∩D(v,α)| × arccos α 2r -(α, r, θ)| dθ -2e -|L(α, r, arccos( α 2r )-π/3)| α 21 -e -α 2 π/6 dα.TakeDv ∩D(v,s) p r, s (u) du = 2 (v,α)|-|L(α,r,arccos( α 2r )-π/3)| 1 -e -α 2 π/6 dα we conclude that p Rn, r, s ≤ e -|Dv∩D(v,s)| + Dv ∩D(v,s) p r, s (u) du (v,α)|-|L(α,r,arccos( α 2r )-π/3)| 1 -e -α 2 π/6 dα which can be bounded by p Rn, r, s ≤ 1 -2 s 0 α e -|Dv ∩D(v,α)|-|L(α,r,arccos( α 2r )-π/3)| × arccos α 2r -π/3 + 2(1-e -α 2 π/6 )

  πr 2 and Π( Ām ) = k>m (2r(2r + 2s) + πr 2 ) k k! e -2r(2r+2s)-πr 2 . Choosing r = s, we have the following bound Π( Ām ) ≤ ((8 + π)r 2 ) m+1 (m + 1)! .

ES 1 ,S 2 ∩ES 2 ,S 1 ΠES 1 ,S 2

 122112 (dϕ) µ p ϕ (B r,s ) ≥ ∩ES 2 ,S 1 ∩FC 1 ∩FC 2 ∩Am Π (dϕ) µ p ϕ (B r,s ) ≥ (1 -ǫ/2)p m > 1 -ǫ = 0.8639because we can choose r, s, m (as preceding) such that


Rainer Nitsche
email: rainer.nitsche@de.bosch.com

Matthias Bitzer
email: matthias.bitzer2@de.bosch.com

Mahmoud El Khaldi

Gérard Bloch
email: gerard.bloch@esstin.uhp-nancy.fr

Fault Tolerant Oxygen Control of a Diesel Engine Air System

Keywords: Fault tolerant control, Diesel engine, Air system, Model-based trajectory planning

This paper is devoted to the fault tolerant control problem of a Diesel engine air system having a jammed Exhaust Gas Recirculation (EGR) valve. The fault tolerant control is based on replaning the trajectory in order to track a new controlled variable which is the oxygen concentration in the intake manifold instead of the fresh air mass flow. The trajectory planning is based on an inverse model approach, utilizing the fundamental thermodynamic relations of the air system.

INTRODUCTION

Diesel engines have made significant progress in terms of performance over the last years and have become very popular and increasingly attractive for passenger cars thanks to their low fuel consumption and high torque at low speed. However, they are significant contributors to overall emissions, especially of nitrogen oxides (NO x ) and particulate matter (PM, also called soot).

The air system of a Diesel engine, shown in Figure 1, is responsible for introducing air into the engine cylinders in order to ensure the fuel combustion. The functioning of the air system is the following. Fresh air entering the air system is compressed by the compressor unit of the turbocharger. The air leaves the compressor with higher pressure, and its temperature increases as a result of the compression. A charge air cooler (CAC), not represented in the Figure, cools down the compressed air before it reaches the throttle valve (TVA; or throttle blade, TB). The fresh air is then mixed with recirculated exhaust gas in the intake manifold and the mixture goes into the engine1 . The waste of the combustion is then expelled from the engine, and a part of the exhaust gas is sent to the intake manifold over the exhaust gas recirculation (EGR) valve to be mixed with fresh air. The rest of the exhaust gas passes through the variable nozzle turbine (VNT) unit of the turbocharger, where its energy is used to drive the compressor wheel via the shaft, which connects the turbine and the compressor wheels.

The main purpose of the exhaust gas recirculation is to reduce the nitrogen oxide (NO x ) emissions. Indeed, mixing the exhaust gas with the fresh air reduces the oxygen concentration in the intake manifold. Thus, since there is less oxygen to react with the fuel in the engine, fewer NO x are produced during the combustion. On the other hand, a lower oxygen rate in the engine leads to an incomplete combustion, and induces the production of more carbon hybrid and particulate matter (PM).

The EGR valve is the main component in the exhaust gas recirculation system. The valve is located between the intake and exhaust manifolds, then a part of the exhaust gas can flow back into the engine. The rate of exhaust gas mixed with the fresh air in the intake manifold is controlled by acting on the cross section area of the valve. However, as any physical system, EGR valves are subject to faults. The main fault that can affect an EGR valve is the blocking at a certain position, which increases the emissions of pollutants.

The paper will deal thus with the fault tolerant control (FTC) of the Diesel engine air system having a jammed EGR valve, with the aim of reducing pollutant emissions. Section 2 introduces the overall faultless control and fault tolerant control structures of the air system. The different parts of the fault tolerant control are described in Sections 3, 4 and 5. Some simulations results are presented in Section 6, before concluding remarks.

AIR SYSTEM CONTROL STRUCTURE

System models and faultless control structure

The air system is controlled by using a particular modelbased approach. Such an approach performs three tasks: planning of trajectories with bounding conditions, calculation of nominal control variables through the inverted system and stabilization of the system around the desired trajectory by closed-loop control. The advantages include the consideration of physical coupling and effects, the simultaneous control of desired values and the fulfillment of high dynamic requirements.

A schematic representation of the air system is presented in Figure 2 with the main describing variables. This system is modelled with the so-called Robert Bosch Standard Simulation Library (RBSL). The resulting air system RBSL model consists of four volumes connected by algebraic relations and has eleven states: -the boost pressure p 21 , -the air mass in the boost pressure volume m 21 , -the intake manifold pressure p 22 , -the air mass in the intake manifold m 22 , -the mass of oxygen in the intake manifold m O2,22 , -the rate of fresh air in the intake manifold r F A,22 , -the exhaust manifold pressure p 3 , -the air mass in the exhaust manifold m 3 , -the downstream pressure of the turbine (VNT) p 4 , -the downstream gas mass of the VNT m 4 , -the angular speed of the turbocharger ω.

The manipulated variables, i.e. the control inputs, are:

the effective area of the throttle valve (TVA) A ef f T V A , -the effective area of the EGR valve A ef f EGR , -the effective area of the VNT A ef f V N T , and the controlled variables: -the fresh air (FA) mass flow ṁF A , -the rate of EGR gas in the intake manifold r EGR,22 , -the boost pressure p 21 .

The model-based control architecture of the air-path is shown in Figure 3. It consists of three modules:

Fault tolerant control structure

The fault-tolerant control strategy aims at keeping this control structure in the faulty case, and not designing new controllers. As the considered fault is an actuator fault (EGR valve jam), there is no possibility to adapt the controller by simply changing its parameters for the faulty case. Instead, alternative actuators which have similar interactions with the system have to be used [START_REF] Blanke | Diagnosis and Fault-Tolerant Control[END_REF]. The alternative actuator used here is the throttle valve, which is controlled by the MCC.

For the faulty case, it is furthermore desired to switch to a new controlled variable which is the oxygen concentration in the intake manifold, because this variable is directly linked with the pollutants emissions.

The resulting closed loop structure is depicted in Figure 4 with references to the main corresponding equations. The first block, detailed in Section 3, calculates the desired oxygen concentration in the intake manifold r d O2,EGR . In order to track this setpoint, the second task, detailed in Section 4, calculates, by inversion of the intake manifold model, the new desired fresh air mass flow ṁd * F A . A reachability analysis (Section 5) allows the MCC to receive ṁd,final

Reachability

Look Up

Tables and ( 20)

Inverse Intake Manifold (16) Oxygen Concentra- tion (3)

Diagnosis

(not treated here)

Estimation Jam Position

(not treated here) The known time variant parameters required for (3) and ( 16) are collected in the vectors ζ(t) and ξ(t). Note that the diagnosis task and the estimation of the jam positition [START_REF] Scheu | Fehlertolerante Frischluftmassenregelung eines PKW-Diesel-Motors[END_REF] are not treated in this paper. Therefore, this jam position is assumed here to be perfectly known.

✲ ✲ ✲ ✲ ✲ ✲ ✲ r ✻ ❄ ❄ ❄ ✛ ✛ ❄ ✻ m d F A p d 2 r d O2,22 ṁd * F A ṁd,final F A EGR, TV VNT A * EGR ζ(t) ξ(t) λ min

NEW CONTROLLED VARIABLE

The new controlled variable, i.e., the oxygen concentration r O2,22 in the intake manifold, is the ratio of the oxygen mass flow ṁO2,eng over the total gas mass flow ṁeng , going in the engine:

r O2,22 = ṁO2,eng ṁeng , (1)
or, as the gas entering the engine is exclusively composed of fresh air and EGR gas:

r O2,22 = ṁF A r O2,F A + ṁEGR r O2,EGR ṁF A + ṁEGR , (2)
where r O2,F A and r O2,EGR are the oxygen ratios in the fresh air and EGR gas, respectively.

This gives the desired oxygen concentration in the intake manifold:

r d O2,22 = ṁd F A r O2,F A + ṁd EGR r d O2,EGR ṁd F A + ṁd EGR , (3)
where ṁd F A is the desired fresh air mass flow (given by the Electronic Control Unit, ECU), and ṁd EGR is the desired EGR mass flow and r d O2,EGR the desired oxygen ratio in the EGR gas, which can be calculated in the following way.

The enthalpy balance in the intake manifold can be written as

c p,M ix ṁ22 T 22 = c p,EGR ṁd EGR T EGR + c p,F A ṁd F A T 21 , (4
) where c p,M ix , c p,EGR and c p,F A are the heating capacities of the mixture (EGR gas + fresh air), EGR gas and fresh air, respectively, T 22 is the intake manifold temperature, T EGR the EGR gas temperature, and T 21 the temperature of the fresh air over the throttle valve.

Assuming that the upstream pressure p 21 of the throttle valve is the same as the pressure p 22 in the intake manifold [START_REF] Scheu | Fehlertolerante Frischluftmassenregelung eines PKW-Diesel-Motors[END_REF], the ideal gas law gives:

ṁ22 T 22 = p 22 Veng R M ix = p d 21 Veng R M ix , ( 5
)
where Veng is the variation of the engine volume, and R M ix is the gas constant of the mixture (EGR gas + fresh air).

Therefore, the desired EGR mass flow ṁd EGR needed to compute (3) can be estimated from ( 4) and ( 5) by

md EGR = c p,M ix p d 21 Veng R M ix -c p,F A ṁd F A T 21 c p,EGR T EGR . (6)
Moreover, with the overall mass balance and the oxygen mass balance in the engine, the desired oxygen ratio r d O2,EGR in the EGR gas, also needed to compute (3), can be obtained by:

r d O2,EGR = r O2,F A ṁd F A -r stoic ṁfuel ṁd F A + ṁfuel , (7)
where r stoic is the stoichiometric air to fuel coefficient and ṁfuel is the amount of injected fuel into the engine.

MODEL-BASED TRAJECTORY PLANNING OF THE DESIRED FRESH AIR MASS FLOW

According to Figure 4, the desired oxygen concentration in the intake manifold can now be used to calculate the new desired fresh air mass flow ṁd * F A in the faulty case. In other words, a new trajectory for the desired fresh air mass flow has to be planned, depending on the known jam-position of the EGR valve and the desired oxygen concentration in the intake manifold. For this, a model-based approach is used, where the inversion of the intake manifold model is based on a flatness analysis in the sense of [START_REF] Fliess | On differentially flat nonlinear systems[END_REF] and [START_REF] Rothfuß | Flachheit: Ein neuer Zugang zur Steuerung und Regelung nichtlinearer Systeme[END_REF].

Inversion of the intake manifold model

For a model-based trajectory planning, the dynamics of the intake manifold can be described with the following state, control input, and output definition:

x = [m O2,22 , m O2,22 ] T , u = ṁF A , y = r O2,22 , (8)
where m O2,22 and m O2,22 are respectively the masses of oxygen and all other gases in the intake manifold, with the external time-varying parameters q = [ ṁEGR , ṁeng , r O2,F A , r O2,EGR , n eng ] T , (9) with n eng the engine speed.

State Equations. The first state equation is obtained by considering the mass flow balance in the intake manifold:

ṁO2,22 = ṁO2,EGR + ṁO2,F A -ṁO2,eng , which leads to ṁO2,22 = ṁEGR r O2,EGR + ṁF A r O2,F A -ṁeng r O2,22 , and finally to ẋ1 = q 1 q 4 + u q 3 -y q 2 . (10) The second state equation is obtained as follows:

ṁO2,22 = ṁ22 -ṁO2,22 = ṁEGR + ṁF A -ṁeng -ẋ1 , leading to ẋ2 = q 1 + uq 2 -ẋ1 , (11) or ẋ2 = q 1 + uq 2 -(q 1 q 4 + u q 3 -y q 2 ).

Output Equation .The system output y is equal to

y = r O2,22 = m O2,22 m O2,22 + m O2,22
,

which gives y = x 1 x 1 + x 2 . ( 12
)
System Inversion. The input u can be expressed as a function of the output y and its successive derivatives. Indeed, (12) gives

x 1 = y(x 1 + x 2 ) = y m 22 .
As neng) ṁeng , where the known function neng) , with V 22 the intake manifold volume, considers the amount of air consumption for the engine, the first state becomes

m 22 = V22 R M ix f (
V22 R M ix f (
x 1 = V 22 R M ix f (n eng ) q 2 y. (13)
Moreover, (12) gives

x 2 = x 1 y -x 1 and ẋ2 = ẋ1 y -ẏ x 1 y 2 -ẋ1 . (14)
Identifying ( 14) with (11) leads to

q 1 + u -q 2 = ẋ1 y -ẏ x 1 y 2 .
Then, using ( 10) and ( 13), simple calculations give finally the output u as

u = q 1 (q 4 -y) - V22 R M ix f (neng) q 2 ẏ y -q 3 . ( 15
)
New desired fresh air mass flow. This leads directly to the expression of the new desired fresh air mass flow:

ṁd * F A = (r O2,EGR -r d O2,22 ) ṁEGR r d O2,22 -r O2,F A - V22 R M ix f(neng) ṁeng ṙd O2,22 r d O2,22 -r O2,F A . ( 16
)
To compute ( 16), one needs to calculate the EGR mass flow ṁEGR , the oxygen concentration r O2,EGR in the EGR gas and the mass flow ṁeng going into the engine.

EGR mass flow

To calculate the mass flow over the jammed EGR valve, the orifice equation can be used:

ṁEGR = p 3 A * EGR 2Π -1 (1 -Π -1 ) T EGR cooler,out R , ( 17
)
where

Π = p 3 p 22 , T EGR cooler,out = T 3 -η HE (T 3 -T coolingf luid ),
and A * EGR is the effective area of the EGR valve estimated by a Kalman filter of a diagnosis algorithm developed by [START_REF] Scheu | Fehlertolerante Frischluftmassenregelung eines PKW-Diesel-Motors[END_REF].

Calculation of the engine mass flow

To estimate the engine mass flow ṁeng , the following equation can be used:

ṁeng = p meas 22 T ASMod 22 R V eng 2 ψ (n eng ) , (18)
where p meas 22 is the measured intake manifold pressure, T ASMod 22 is the ASMod estimation of the intake manifold temperature, p meas 22 is the measured intake manifold pressure, T ASMod 22 is the ASMod estimation of the intake manifold temperature, V eng is the engine volume, ψ (n eng ) is a correction term representing the engine volumetric efficiency, R is the perfect gases constant, and c p the heating capacity of the intake gas at constant pressure is the measured intake manifold pressure.

Calculation of the oxygen rate in the exhaust gas

To estimate the oxygen rate r O2,EGR in the exhaust gas, ( 7) can be used by replacing the desired values by the real values:

r O2,EGR = r O2,F A ṁmeas F A -r stoic ṁfuel ṁmeas F A + ṁfuel , (19)
with r stoic = 14.5.

Test of the model-based trajectory planning in the faultless case

Equation ( 16) has been implemented in the simulation environment and tested using a special test cycle in the faultless case. Since the faultless case is considered, the old (given by the ECU) and the new desired fresh air mass flow should be equal. Figure 5

REACHABILITY ANALYSIS OF THE DESIRED FRESH AIR MASS FLOW AND CLOSED-LOOP STRUCTURE

The next step (see Figure 4) is to check if the values of the calculated new desired fresh air mass flow ṁd * F A (16) are physically reachable (higher limit), and not below a certain limit (lower limit). This reachability analysis is decribed below.

Higher limit. The fresh air mass flow depends on several variables: engine speed, injected fuel quantity, EGR valve position, turbocharger VNT position, TVA position. However, it is obvious that the throttle valve (TVA) has to be fully opened to reach a maximal value. Therefore, a look-up table is created with only these 4 first variables as inputs, and the maximal reachable fresh air mass flow as output. The output values of this table are the fresh air mass flow values obtained in the faultless case from appropriate simulation runs, where

• the engine speed varies between 1500rpm and 4000rpm with a step of 250; • the injected fuel quantity varies between 10mg/hub and 55mg/hub with a step of 5; • the EGR valve position varies between 0% (fully opened) and 100% (fully closed) with a step of 10; • the VNT position varies between 0% and 100% with a step of 2.

Lower limit. The lower limit of the fresh air mass flow is governed by the air to fuel ratio λ which has to be greater than a certain limit λ min :

λ = ṁF A r stoic ṁfuel ≥ λ min = 1.3, (20)
where r stoic = 14.5 is the stoichiometric air/fuel ratio.

Figure 6 shows an example of the desired fresh air mass flow in the faultless case, with the maximal reachable and minimal allowed values. new desired h ig h er limit lower limit Figure 6. The new desired fresh air mass flow ṁd * F A , with the maximal reachable and minimal allowed values.

APPLICATION TO DIFFERENT JAM POSITIONS

Several simulation results are given here to illustrate the behavior of the proposed fault tolerant control.

Faulty case 25%. When the EGR valve is blocked at 25% position (widely opened), the first plot in Figure 7 shows that the system would require much more fresh air (blue line) than the physical limit (green line) in order to reach the desired oxygen concentration in the intake manifold.

A closer look at the second plot, which compares the measured oxygen concentrations in the intake manifold, shows that, most of the time, the control with model-based trajectory planning of the desired fresh air mass flow gives similar results than without it, except on the low operating points where the FTC is able to reach the desired oxygen concentration in the intake manifold.

Faulty case 50%. As in the previous case, the system needs sometimes more fresh air than the maximal physically reachable value in order to reach the desired oxygen concentration in the intake manifold (Figure 8, first plot). However, in spite of the physical limitation of the system, the model-based trajectory planning of the desired fresh air mass flow makes the system reach the desired oxygen concentration in the intake manifold in the most operating regions (Figure 8, second plot), except in the full-load regions where the desired oxygen concentration requires to have only fresh air in the intake manifold, and that is not possible when the EGR valve is blocked at a middle position.

Faulty case 75%. In this case, the EGR valve is blocked at 75% position, which means almost closed. Thus, less EGR gas comes into the intake manifold and the system needs less fresh air to reach the desired oxygen concentration in the intake manifold. However, the system is blocked by the λ min -limit (20) and therefore it is constrained to have more fresh air than desired in the intake manifold (Figure 9, first plot). As a consequence, the measured oxygen concentration in the intake manifold with FTC is greater than the desired one; however, it is lower than without FTC (Figure 9, second plot).

CONCLUSION

The model-based trajectory planning of the desired fresh air mass flow has shown many advantages. It allows to keep the same control structure, which already works pretty well in the faultless case. Furthermore, regarding the simulation results for different jam positions of the EGR valve, it has come out that this fault-tolerant control approach gives satisfactory results. Indeed, it makes the system to get closer to the desired oxygen concentration in the intake manifold than without any replanning of the fresh air mass flow trajectory.

Another additional actuator could be used in parallel to the throttle valve to improve the FTC performances. Indeed, the Variable Nozzle Turbine influences the boost pressure and therefore the oxygen concentration in the intake manifold. This would require a trajectory planning of the desired boost pressure [START_REF] El Khaldi | Fault-Tolerant Control of the Diesel Engine Air system[END_REF], in an analogous way of the replanning of the desired fresh air mass flow.

Figure 1 .

 1 Figure 1. Air system of the Diesel engine.

Figure 2 .

 2 Figure 2. Schematic representation of the air path.

Figure 3 .

 3 Figure 3. Current air system control unit.

Figure 4 .

 4 Figure 4. Overall closed-loop structure of the proposed fault tolerant control.

Figure 5

 5 Figure 5. Old ṁd F A and new ṁd * F A desired fresh air mass flows, in the faultless case.

  flow [g/s]

Figure 7 .

 7 Figure 7. Faulty case 25% at t = 50s. First plot: new desired fresh air mass flow with the maximal and minimal reachable values. Second plot: oxygen concentrations in the intake manifold, desired and measured (with and without FTC trajectory planning).

Figure 8 .Figure 9 .

 89 Figure 8. Faulty case 50% at t = 50s. First plot: new desired fresh air mass flow with the maximal and minimal reachable values. Second plot: oxygen concentrations in the intake manifold, desired and measured (with and without FTC trajectory planning).

  . The two-degree of freedom design is based on an inverse model of the turbocharger for the feedforward part of the controller and a PIstructure for the feedback part. The MBC uses a family of PID controllers, each of them providing a satisfactory control for a particular operating point (gain scheduling). It can work either in open or closed loop, depending on the operating point. In closed loop mode, the MBC regulates the boost pressure to avoid getting too high pressure that may damage the turbocharger.• The Model-based Charge Control (MCC) controls the fresh air mass flow ṁF A and/or the EGR rate r

	• The Air System Model (ASMod) estimates the air
	system signals that cannot be measured with the
	available sensors. From the modelling of the system, it
	also dynamically corrects the measured signals, when
	the measured value is in the range of low sensitivity
	of the sensor. The outputs of the ASMod are used by
	the model-based controllers, briefly described next.
	• The Model-based Boost Control (MBC) controls
	the boost pressure p 21 by manipulating the VNT input A ef f V N T

EGR , depending on the operating point, by acting on the throttle and EGR valves. Again, a two-degree of freedom design is used. The MCC contains Kalman filters which estimate the current areas of the throttle and EGR valves, an inverse model of the intake manifold to calculate the desired areas, and PI-controllers.

In this paper, the term engine stands for the combustion chamber (cylinders) and not for the overall engine.


C Thomas

P Blanc

T Ranchin

L Wald

GEOMETRIC QUALITY OF IMAGES: ESTIMATION OF THE MTF

INTRODUCTION

The knowledge of the resolution of a sensor is an important characteristic which allows the analysis of their performance. The spatial resolution remains an obscure domain and was tried to be many times defined. The resolution is what characterize the capacity of an imaging system to describe how much details are visible in the image, depending on their size. Image quality is often characterized by well-known measures and criteria, usually linked with radiometry, like mean, contrast, brightness, noise variance, radiometric resolution or granularity. Only a few tools are available to characterize the spatial performance of an imaging device. The most often used are Point Spread Function-like functions [START_REF] Li | Spatial quality evaluation of fusion of different resolution images[END_REF], gradient or high-pass filters in general [START_REF] Zhou | A wavelet transform method to merge Landsat TM and SPOT panchromatic data[END_REF], Ground Sampled Distance (GSD) [START_REF] Vrabel | Multispectral imagery band sharpening study[END_REF][START_REF] Wald | Data Fusion: Definitions and Architectures.Fusion of Images of Different Spatial Resolutions[END_REF], Full Width at Half Maximum (FWHM), or visual analyses [START_REF] Irars | Civil NIIRS Reference Guide, Imagery Resolution Assessment and Reporting Standards Committee[END_REF]. The problem encountered with selecting a certain band width by original image numerical filtering is that only the extracted spatial structures are highlighted. The spatial features of other spatial frequency are neglected. The idea would be to use a function or a criteria which characterize the overall geometrical quality of the images delivered by an instrument. The problem of specifying resolution and perceived sharpness was solved with the introduction of the Modulation Transfer Function (MTF), a precise measurement made in frequency domain. By definition, the MTF is the Fourier transform of the instrument answer to a pulse. It subdues the contrast of the image with the frequency, caused by finite detector size, sensor motion, diffraction, aberrations, atmospheric scattering, turbulence, and electronic effects [START_REF] Holst | Electro-optical imaging system performance[END_REF]. Only the mean of the image is perfectly conserved, and the MTF tends to gather all numerical accounts around this mean. If the spatial frequency of an object is higher than the cut-off frequency, its image will be characterize by a uniform enlightenment. In the field of Earth observation, MTF is an information that, in most cases, satellite providers don't deliver easily. Representative of the intrinsic quality of their device products, the shape of the MTF and more particularly the cut-off frequency is not clearly given or even defined. Distortions and limit of sensors are well guarded. For instance, the documentation of SPOT Image (Anonymous 1986) remains very vague on the subject of MTF. They give specifications for the cut-off frequency of SPOT P images without informing about the way they calculated it and no information was found about multispectral images. Moreover, some confusions between terms like, MTF, resolving power, cut-off frequency, minimum allowable MTF at the Nyquist frequency, create inconvenient ambiguity. Some are looking for the highest resolution pattern where detail was visible, whereas others define a percentage of attenuation of 10, 5 or even 2% of the original contrast, but no standard threshold is determined. The objective of this paper is to use the MTF in order to assess the overall image geometrical quality or resolution of an imaging system. As this function is not always clearly described by manufacturers, the MTF should be directly estimated from raw satellite images. By using the MTF, two systems can readily be compared. In this paper, a method for this estimation is presented. An example is provided, dealing with IKONOS imagery. Ryan et al. (2003) estimated MTF from real images where constructed targets were settled in order to check the matching of the original IKONOS imagery MTF specifications. We compare our estimation of MTF to their results. Some perspectives on the use of MTF as a quality criteria are proposed.

BACKGROUND

Several methods exist to determine the MTF. The MTF is the amplitude spectrum of the system answer to an perfect impulsion, or point spread function PSF. For an optic system, the stain of diffusion which stems from a perfect luminous point (called Dirac impulsion) corresponds to the PSF of the system. The image of this object will be at least one pixel wide. Instead, it will normally consist of a spot of several pixels, brightest in the centre and progressively darker away from the centre, and is deeply related to sampling. It describes how the sharp edge has been spread out by the imaging process. Now, the resolution can be also defined as the width within which the PSF drops to half the maximal value, called full width half maximum FWHM, but give only one single value of the entire image spatial quality. The PSF need not be symmetrical, so there may be different spatial resolutions in different directions.

Drawing a point or even a cloud of points is conceivable for a scanner of a numerical camera. This becomes far more unconceivable for satellite images, even with constructed targets. A solution to since former problem was to simulate the pulse from the available information. The traditional method of calculating the MTF involves taking an edge scan of an image perpendicular to the edge (edge spread function ESF). The edge is not an instantaneous transition and appears to be more gradual, and give some information about how the imaging system treats edges. To create the pulse, the edge is derivated leading to the Line Spread Function or PSF if a profile perpendicular to the edge is plotted. The Fourier transform (FT) of the PSF provides an estimate of the MTF in one orientation. As an edge will only contain a few points and to locate the inflexion point (and consequently the location of the pulse response) we would need to oversample the edge. This can be done by the use of a long enough tilted slop instead of an horizontal or vertical one. In our natural environment, however, finding a perfect tilted sloped such as long bridges, big roads, large buildings and other large size linear features, is not very common. So Ryan et al. (2003) chose to use artificial targets located in large fields, where they controlled the angle of orientation of the tarp to evaluate MTF of IKONOS sensors. Nevertheless, despite the efficiency of the method, logistical difficulties with artificial targets and lack of control of natural ones, make target-free approaches to image-quality evaluation desirable.

EXPERIMENT

Following the idea of Luxen et Förstner (2002), a systematic method to read MTF is defined from any image. The needed PSF is derived from the sharpness of the edges located in the image. In order to alleviate the lack of perfect contours in nature, a statistical approach of the problem was opted. A number of significant amplitude edges is selected and oversampled in each direction with a bicubic interpolation, and in both direction (vertical and horizontal). This selection is performed using the two corresponding images in the directional multiresolution decomposition of [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF]. Then, a numerical derivation is applied to the edges, with the following numerical filter : [-1, 0, 1], and the corresponding transposed for the other direction, in order to get the pulse function. Finally, Fourier transform is applied to the pulse and the MTF for one profile is obtained. For a better estimation, all the maximum edges are processed the same way. The MTF at the Nyquist frequency is expected to be lower in the vertical direction, caused by the scanning system motion.

EVALUATION AND PERSPECTIVES

A comparison between our estimation of the MTF from an IKONOS image over the area of Frederiktown (Canada) and those obtained by Ryan et al. (2003) is presented. A study of differences existing in the MTF between the two directions in the image is proposed. Based on MTF estimation, some quality criteria are envisaged and discussed. These criteria will add on the process of describing the intrinsic quality of a sensor.