Giovanni Chierchia

Nelly Pustelnik
email: nelly.pustelnik@ens-lyon.fr

Jean-Christophe Pesquet
email: jean-christophe.pesquet@univ-paris-est.fr

Béatrice

Béatrice Pesquet-Popescu

† N Pustelnik

Epigraphical Projection and Proximal Tools for Solving Constrained Convex Optimization Problems: Part I

published or not. The documents may come

Introduction

As an offspring of the wide interest in frame representations and sparsity promoting techniques for data recovery, there has been an emergence of proximal methods to efficiently solve convex optimization problems arising in inverse problems. Examples of applications of these methods can be found in f-MRI reconstruction [START_REF] Chaari | A wavelet-based regularized reconstruction algorithm for SENSE parallel MRI with applications to neuroimaging[END_REF][START_REF] Guerquin-Kern | A fast wavelet-based reconstruction method for magnetic resonance imaging[END_REF], satellite image restoration [START_REF] Facciolo | Irregular to regular sampling, denoising, and deconvolution[END_REF][START_REF] Hajlaoui | Satellite image restoration in the context of a spatially varying point spread function[END_REF], microscopy image deconvolution [START_REF] Dupé | A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF][START_REF] Jezierska | A primal-dual proximal splitting approach for restoring data corrupted with Poisson-Gaussian noise[END_REF], computed tomography [START_REF] Vandeghinste | Split-Bregman-based sparse view CT reconstruction[END_REF], Positron Emission Tomography [START_REF] Pustelnik | Parallel algorithm and hybrid regularization for dynamic PET reconstruction[END_REF][START_REF] Anthoine | Some proximal methods for CBCT and PET tomography[END_REF], texture-geometry decomposition [START_REF] Aujol | Image decomposition into a bounded variation component and an oscillating component[END_REF][START_REF] Aujol | Structure-texture image decomposition -modeling, algorithms, and parameter selection[END_REF][START_REF] Briceño-Arias | Proximal algorithms for multicomponent image recovery problems[END_REF], machine learning [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF][START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF], stereo vision [START_REF] Chaux | A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation[END_REF], and audio processing [START_REF] Kowalski | Beyond the narrowband approximation: Wideband convex methods for under-determined reverberant audio source separation[END_REF][START_REF] Akyildiz | An analysis prior based decomposition method for audio signals[END_REF]. Proximal algorithms have gained much popularity in solving large-size optimization problems involving non-differentiable (or even non finite) functions. One of the main advantages of these methods is that they are amenable to parallel implementations. For a survey on proximal algorithms and their applications, the reader is referred to [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF]. Note also that some of these methods are closely related to augmented Lagrangian approaches [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF][START_REF] Figueiredo | Restoration of Poissonian images using alternating direction optimization[END_REF][START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF]. Even if proximal algorithms and the associated convergence properties have been deeply investigated [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF][START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF][START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF][START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF], some questions persist in their use for solving inverse problems. A first question is: how can we set the parameters serving to enforce the regularity of the solution in an automatic way? Another question is related to the selection of the most appropriate algorithm within the class of proximal algorithms for a given application. This also raises the question of the computation of the proximity operators associated with the different functions involved in the criterion. Various strategies were proposed in order to address the first question [START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF][START_REF] Hansen | The use of the L-curve in the regularization of discrete ill-posed problems[END_REF][START_REF] Pizurica | Estimating the probability of the presence of a signal of interest in multiresolution single-and multiband image denoising[END_REF][START_REF] Ramani | Monte-Carlo SURE: A black-box optimization of regularization parameters for general denoising algorithms[END_REF][START_REF] Chaari | A hierarchical bayesian model for frame representation[END_REF], but the computational cost of these methods is often high, especially when several regularization parameters have to be set. Alternatively, it has been recognized for a long time that incorporating constraints directly on the solutions [START_REF] Youla | Image restoration by the method of convex projections. Part I -theory[END_REF][START_REF] Trussell | The feasible solution in signal restoration[END_REF][START_REF] Combettes | Inconsistent signal feasibility problems : least-squares solutions in a product space[END_REF][START_REF] Kose | Filtered variation method for denoising and sparse signal processing[END_REF][START_REF] Teuber | Minimization and parameter estimation for seminorm regularization models with I-divergence constraints[END_REF] instead of considering regularized functions may often facilitate the choice of the involved parameters. Indeed, in a constrained formulation, the constraint bounds are usually related to some physical properties of the target solution or some knowledge of the degradation process, e.g. the noise statistical properties. Note also that there exist some conceptual Lagrangian equivalences between regularized solutions to inverse problems and constrained ones, although some caution should be taken when the regularization functions are nonsmooth (see [START_REF] Ciak | Homogeneous penalizers and constraints in convex image restoration[END_REF] where the case of a single regularization parameter is investigated).

In this context, the objective of this paper is to provide an answer to the second question, that is to propose efficient parallel algorithms for solving the following constrained convex optimization problem:

Problem 1.1. minimize x∈H R r=1 g r (T r x) s. t.        H 1 x ∈ C 1 , . . . H S x ∈ C S , (1)
where (i) for every r ∈ {1, . . . , R}, T r is a bounded linear operator from a real Hilbert space H to R Nr ;

(ii) for every r ∈ {1, . . . , R}, g r is a proper lower-semicontinuous convex function from R Nr to ]-∞, +∞];

(iii) for every s ∈ {1, . . . , S}, H s is a bounded linear operator from H to R Ms ;

(iv) for every s ∈ {1, . . . , S}, C s is a nonempty closed convex subset of R Ms .

A large number of inverse problems can be formulated under the form of Problem 1.1. A classical application example is image recovery from blurred and noisy observations. Let z denote the vector of observed data. When the noise is assumed to be zero-mean additive white Gaussian, the restored data can be obtained by solving Problem 1.1 where R = 2, S = 0, g 1 = λ • 1 with λ ∈ ]0, +∞[, T 1 is an analysis frame operator [START_REF] Jacques | A panorama on multiscale geometric representations, intertwining spatial, directional and frequency selectivity[END_REF] allowing us to sparsify signal x, g 2 = • -z 2 , and T 2 denotes the matrix associated with the degradation blur [START_REF] Figueiredo | An EM algorithm for wavelet-based image restoration[END_REF][START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]. On the other hand, a constrained formulation of the same restoration problem leads to Problem 1.1 with R = 1, S = 1, g 1 = • 1 , T 1 is the same frame operator as previously, and

(∀x ∈ H)

H 1 x ∈ C 1 ⇔ H 1 x -z 2 ≤ η 1 .
Here, H 1 denotes the matrix associated with the degradation blur, and η 1 is a positive constant which is typically chosen proportional to the noise variance. This specific constraint formulation has been considered in [START_REF] Afonso | An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems[END_REF].

The present work aims at designing efficient methods in order to deal with Problem 1.1 when the convex constraints are expressed as follows: for every s ∈ {1, . . . , S},

(∀x ∈ H) H s x ∈ C s ⇔ h s (H s x) ≤ η s , (2)
where η s ∈ R and h s is a proper lower-semicontinuous function from R Ms to ]-∞, +∞].

The paper is organized as follows. In Section 2, we motivate the choice of proximal tools and recall some of their theoretical properties. In addition, closed form expressions for specific proximity operators are derived. Then, in order to deal with a constraint expressed under the general form (2), a splitting approach involving an epigraphical projection is proposed. This epigraphical projection technique is described in detail in Section 3. Experiments in the context of image reconstruction are presented in Section 4. Finally, some conclusions are drawn in Section 5.

Notation: Let H be a real Hilbert space endowed with the norm • and the scalar product

• | • . Γ 0 (H) denotes the set of proper lower-semicontinuous convex functions from H to ]-∞, +∞]. Recall that a function ϕ : H → ]-∞, +∞] is proper if its domain dom ϕ = y ∈ H ϕ(y) < +∞ is nonempty. The epigraph of ϕ ∈ Γ 0 (H) is the nonempty closed convex subset of H × R defined as epi ϕ = (y, ζ) ∈ H × R ϕ(y) ≤ ζ and the lower level set of ϕ at height ζ ∈ R is the nonempty closed convex subset of H defined as lev ≤ζ ϕ = y ∈ H ϕ(y) ≤ ζ . A subgradient of ϕ at y ∈ H is an element of its subdifferential defined as ∂ϕ(y) = t ∈ H (∀u ∈ H) ϕ(u) ≥ ϕ(y) + t | u -y . 1 The indicator function ι C ∈ Γ 0 (H) of a nonempty closed convex subset C of H is given by (∀y ∈ H) ι C (y) = 0, if y ∈ C, +∞, otherwise. (3)
The relative interior of a subset C of H is denoted by ri C.

2 Proximal tools

From gradient descent to proximal algorithms

The first methods for finding a solution to an inverse problem were restricted to the use of a differentiable cost function [START_REF] Tikhonov | Tikhonov regularization of incorrectly posed problems[END_REF], i.e. Problem 1.1 where S = 0 and, for every r ∈ {1, . . . , R}, g r denotes a differentiable function.

In this context, gradient-based algorithms, e.g. nonlinear conjugate gradient or quasi-Newton methods, are popular (see [START_REF] Chouzenoux | A Majorize-Minimize strategy for subspace optimization applied to image restoration[END_REF] and references therein for recent developments concerning these approaches). However, in order to model additional properties, sparsity promoting penalizations or hard constraints (S ≥ 1) may be introduced and the diffentiability property is not satisfied anymore. One way to circumvent this difficulty is to resort to smart approximations in order to smooth the involved non-differentiable functions [START_REF] Aubert | Sur la minimisation d'une fonctionnelle non convexe, non différentiable en dimension 1[END_REF][START_REF] Ben-Tal | A smoothing technique for nondifferentiable optimization problems[END_REF][START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms, Part I : Fundamentals[END_REF].

If one wants to address the original nonsmooth problem without introducing approximation errors, one may apply some specific algorithms e.g. Gauss-Seidel or Uzawa methods, the convergence of which is guaranteed under restrictive assumptions [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF]. Interior point methods [START_REF] Wright | Primal-dual interior-point methods[END_REF] can also be employed for small to medium size optimization problems.

On the other hand, in order to solve convex feasibility problems, i.e. to find a vector belonging to the intersection of convex sets (Problem 1.1 with R = 0), iterative projection methods were developed. The projection onto convex sets algorithm (POCS) is one of the most popular approach to solve data recovery problems [START_REF] Bregman | The method of successive projection for finding a common point of convex sets[END_REF][START_REF] Gurin | Projection methods for finding a common point of convex sets[END_REF][START_REF] Youla | Image restoration by the method of convex projections. Part I -theory[END_REF][START_REF] Combettes | The foundations of set theoretic estimation[END_REF]. A drawback of POCS is that it is not well-suited for parallel implementations. The Parallel Projection Method (PPM) and Method of Parallel Projections (MOPP) are variants of POCS making use of parallel projections. Moreover, these algorithms were designed to efficiently solve inconsistent feasibility problems (when the intersection of the convex set is empty). Thorough comparisons between projection methods have been performed in [START_REF] Combettes | Convex set theoretic image recovery by extrapolated iterations of parallel subgradient projections[END_REF][START_REF] Censor | On the effectiveness of projection methods for convex feasibility problems with linear inequality constraints[END_REF].

Computing the projection P C onto a nonempty closed convex subset C of a real Hilbert space H requires to solve a constrained quadratic minimization problem:

(∀y ∈ H) P C (y) = argmin u∈C u -y . (4)
The distance to C of every point y ∈ H is then given by d C (y) = y -P C (y) . However, it turns out that a closed form expression of the solution to (4) is available in a limited number of instances. One such well-known example is the projection onto a closed half-space:

P C (y) =    y, if y | t ≤ η, y + η -y | t t 2 t, otherwise. (5)
When an expression of the direct projection is not available, the convex set C can be approximated by a half-space, which leads to the concept of subgradient projection. Proposition 2.2. [START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms, Part I : Fundamentals[END_REF] Let η ∈ R and let ϕ ∈ Γ 0 (H). Suppose that C = lev ≤η ϕ = ∅. Let y ∈ H and let t ∈ ∂ϕ(y) be a subgradient of ϕ at y. Then, C is a subset of the half-space

C y η = u ∈ H u -y | t ≤ η -ϕ(y) . (6)
If t = 0, the projection of y onto C y η is a subgradient projection of y ∈ H onto C.

The subgradient projection plays a key role in Polyak algorithm that alternates a subgradient projection and an exact projection [START_REF] Polyak | Minimization of unsmooth functionals[END_REF] (see also [START_REF] Shor | Minimization Methods for Non-differentiable Functions[END_REF] for alternative projected subgradient approaches). An efficient block iterative surrogate splitting method was proposed in [START_REF] Combettes | A block-iterative surrogate constraint splitting method for quadratic signal recovery[END_REF] in order to solve Problem 1.1 when, for every r ∈ {1, . . . , R}, g r = • -z r 2 where z r ∈ R Nr . A main limitation of this method is that the global objective function must be strictly convex. For recent works about subgradient projection methods, the readers may refer to [START_REF] Slavakis | The adaptive projected subgradient method over the fixed point set of strongly attracting nonexpansive mappings[END_REF][START_REF] Bouboulis | Adaptive learning in complex reproducing kernel Hilbert spaces employing Wirtinger's subgradients[END_REF].

A way to overcome this difficulty consists of considering proximal approaches. The key tool in these methods is the proximity operator [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] of a function ϕ ∈ Γ 0 (H), defined as

(∀y ∈ H) prox ϕ (y) = argmin u∈H 1 2 u -y 2 + ϕ(u). (7)
The proximity operator generalizes the notion of projection onto a nonempty closed convex subset C of H in the sense that prox ιC = P C . For every y ∈ H, p = prox ϕ (y) is uniquely defined through the inclusion

y -p ∈ ∂ϕ(p). (8)
Proximity operators enjoy many additional interesting properties. Some of them are recalled next.

Example 2.3. [23] Let τ ∈ ]0, +∞[, β ∈ [1, +∞[, and set ϕ : R → ]-∞, +∞] : ξ → τ |ξ| β . (9)
Then, for every

ξ ∈ R, prox ϕ (ξ) is given by                  sign(ξ) max{|ξ| -τ, 0}, if β = 1, ξ + 4τ 3 . 2 1/3 (ǫ -ξ) 1/3 -(ǫ + ξ) 1/3 , if β = 4 3 , ξ + 9τ 2 sign(ξ) 8 1 -1 + 16|ξ| 9τ 2 , if β = 3 2 , ξ 1+2τ , if β = 2, sign(ξ) √ 1+12τ |ξ|-1 6τ , if β = 3, (10)
where ǫ = ξ 2 + 256τ 3 /729 and sign is the signum function.

It can be noticed that the proximity operator associated with β = 1 reduces to a soft thresholding [START_REF] Donoho | De-noising by soft-thresholding[END_REF].

The class of proximal methods includes primal algorithms [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF][START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF][START_REF] Figueiredo | Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems[END_REF][START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Fornasier | Subspace correction methods for total variation and ℓ 1minimization[END_REF][START_REF] Steidl | Removing multiplicative noise by Douglas-Rachford splitting methods[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF] and primaldual algorithms [START_REF] Chen | A proximal-based decomposition method for convex minimization problems[END_REF][START_REF] Esser | A general framework for a class of first order primal-dual algorithms for convex optimization in imaging science[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Briceño-Arias | A monotone + skew splitting model for composite monotone inclusions in duality[END_REF][START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallel-sum type monotone operators[END_REF][START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF][START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF]. Primal algorithms generally require to compute inverses of some linear operators (typically, R r=1 T * r T r + S s=1 H * s H s ), while primal-dual ones only require to compute (T r ) 1≤r≤R and (H s ) 1≤s≤S , and their adjoints (T * r ) 1≤r≤R and (H * s ) 1≤s≤S . Consequently, primal-dual methods are often easier to implement than primal ones, but their convergence may be slower [START_REF] Pustelnik | Relaxing tight frame condition in parallel proximal methods for signal restoration[END_REF][START_REF] Couprie | Constrained TV-based regularization on graphs[END_REF].

Proximity operators: new closed forms

The projection onto a convex set such as defined in (2) often is non trivial. In Section 3, we will show that this problem can be solved by resorting to a set of epigraphical projections which are easier to compute. The key point is that epigraphical projections are closely related to proximity operators. In the following, we provide some results that will be useful for the calculation of the projection onto the epigraph of a convex function (all the proofs are provided in the appendix).

Proposition 2.4. Let H be a real Hilbert space and let H × R be equipped with the standard product space norm. Let ϕ be a function in Γ 0 (H) such that dom ϕ is open. The projector P epi ϕ onto the epigraph of ϕ is given by:

∀(y, ζ) ∈ H × R P epi ϕ (y, ζ) = (p, θ) (11)
where p = prox1 2 (max{ϕ-ζ,0}) 2 (y), θ = max{ϕ(p), ζ}. (12)
Note that alternative characterizations of the epigraphical projection can be found in [73, Propositions 9.17, 28.28].

Let ϕ be a function in Γ 0 (H) such that dom ϕ is open. From Proposition 2.4, we see that prox 1 2 (max{ϕ-ζ,0}) 2 with ζ ∈ R plays a prominent role in the calculation of the projection onto epi ϕ. We now provide examples of functions ϕ for which this proximity operator admits a simple form.

Proposition 2.5. Let β ∈ [1, +∞[, τ ∈ ]0, +∞[. Assume that (∀y ∈ R) ϕ(y) = τ |y| β . ( 13
)
If ζ ∈] -∞, 0], then for every y ∈ R prox 1 2 (max{ϕ-ζ,0}) 2 (y) =    sign(y) 1 + τ 2 max{|y| + τ ζ, 0}, if β = 1, sign(y)χ 0 , if β > 1, ( 14
)
where χ 0 is the unique solution on [0, +∞[ of the equation

βτ 2 χ 2β-1 -βτ ζχ β-1 + χ = |y|. ( 15
)
If ζ ∈ ]0, +∞[, then, for every y ∈ R, prox 1 2 (max{ϕ-ζ,0}) 2 (y) = y, if |y| ≤ ζ/τ 1/β , sign(y)χ (ζ/τ ) 1/β , otherwise, (16)
where [START_REF] Chaux | A parallel proximal splitting method for disparity estimation from multicomponent images under illumination variation[END_REF].

χ (ζ/τ ) 1/β is the unique solution on [( ζ τ ) 1/β , +∞[ of (
Note that, when β is a rational number, ( 15) is equivalent to a polynomial equation for which either closed form solutions are known or standard numerical solutions exist. Proposition 2.6. Let H be a real Hilbert space and let C be a nonempty convex subset of H. Let β ∈ [1, +∞[, τ ∈ ]0, +∞[, and ζ ∈ R. Then, for every y ∈ H,

prox 1 2 (max{τ d β C -ζ,0}) 2 (y) = y, if y ∈ C, αy + (1 -α)P C (y), otherwise, (17)
where α = prox 1 2 (max{τ |•| β -ζ,0}) 2 dC(y) dC(y)
and the expression of

prox 1 2 (max{τ |•| β -ζ,0}) 2 is provided by Proposition 2.5.
The following result is a consequence of the previous one: Corollary 2.7. Let H be a real Hilbert space. Assume that

(∀y ∈ H) ϕ(y) = τ y -z (18)
where z ∈ H, ζ ∈ R and τ ∈ ]0, +∞[. Then, for every y ∈ R,

prox 1 2 (max{ϕ-ζ,0}) 2 (y) =      z, if y -z < -τ ζ, y, if y -z < ζ τ , z + α(y -z), otherwise, (19)
where α = 1 1+τ 2 1 + τ ζ y-z .
Another result which will be used is the following:

Proposition 2.8. Let M ∈ N * and let ϕ be defined as

(∀y ∈ R) ϕ(y) = 1 2 M m=1 max{τ (m) (ν (m) -y), 0} 2 (20)
where (τ (1) , . . . , τ (M) ) ⊤ ∈ R M , and (ν (1) , . . . , ν

(M) ) ⊤ ∈ R M is such that ν (1) ≤ . . . ≤ ν (M) . Set ν (0) = -∞ and ν (M+1) = +∞.
Then, ϕ ∈ Γ 0 (R), and for every y ∈ R,

prox ϕ (y) = y + m-1 m=1 ν (m) (τ (m) -) 2 + M m=m ν (m) (τ (m) + ) 2 1 + m-1 m=1 (τ (m) -) 2 + M m=m (τ (m) + ) 2 , ( 21
)
where, for every m ∈ {1, . . . , M }, τ

(m) - = min(τ (m) , 0) and τ (m) +
= max(τ (m) , 0), and m is the unique integer in {1, . . . , M + 1} such that

ν (m-1) < y + m-1 m=1 ν (m) (τ (m) -) 2 + M m=m ν (m) (τ (m) + ) 2 1 + m-1 m=1 (τ (m) -) 2 + M m=m (τ (m) + ) 2 ≤ ν (m) (22)
(with the convention

0 m=1 • = M m=M+1 • = 0).

Projection computation

We now turn our attention to convex sets for which the associated projection does not have a closed form, and we show that under some appropriate assumptions, it is possible to circumvent this difficulty. Let C denote such a nonempty closed convex subset of R M and assume that, for every y ∈ R M ,

y ∈ C ⇔ h(y) = L ℓ=1 h (ℓ) (y (ℓ) ) ≤ η (23)
where η ∈ R. Hereabove, the generic vector y has been decomposed into blocks of coordinates as follows y = [(y (1) ) ⊤ size M (1) , . . . , (y

(L) ) ⊤ size M (L) ] ⊤ (24)
and, for every ℓ ∈ {1, . . . , L},

y (ℓ) ∈ R M (ℓ) and h (ℓ) is a function in Γ 0 (R M (ℓ) ) such that ri(dom h (ℓ) ) = ∅.
Actually, any closed convex subset C can be expressed in this way by setting η = 0, L = 1 and h = d C .

By introducing now the auxiliary vector

ζ = ζ (ℓ) 1≤ℓ≤L ∈ R L , Condition (23) can be equivalently rewritten as 2        L ℓ=1 ζ (ℓ) ≤ η, (∀ℓ ∈ {1, . . . , L}) h (ℓ) (y (ℓ) ) ≤ ζ (ℓ) . (25) (26)
Let us now introduce the closed half-space of R L defined as

V = ζ ∈ R L 1 ⊤ L ζ ≤ η , (27)
with 1 L = (1, . . . , 1) ⊤ ∈ R L , and the closed convex set

E = (y, ζ) ∈ R M × R L (∀ℓ ∈ {1, . . . , L}) (y (ℓ) , ζ (ℓ) ) ∈ epi h (ℓ) . (28)
Then, Constraint [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] means that ζ ∈ V , whereas Constraint ( 26) is equivalent to (y, ζ) ∈ E. In other words, the constraint set C can be split into the two constraint sets V and E provided that L additional scalar variables (ζ (ℓ) ) 1≤ℓ≤L are introduced in the original optimization problem. Dealing with additional constraints in the original problem is not a problem for proximal splitting algorithms as far as the projections onto the associated constraint sets can be computed. An example of use of such an algorithm will be described in detail in Section 4 .

In the present case, the projection onto V is given by Proposition 2.1, whereas the projection onto E is given by

(∀(y, ζ) ∈ R M × R L ) P E (y, ζ) = (p, θ) (29)
where θ = (θ (ℓ) ) 1≤ℓ≤L , vector p ∈ R M is blockwise decomposed as p = [(p (1) ) ⊤ , . . . , (p (L) ) ⊤ ] ⊤ like in [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF], and

(∀ℓ ∈ {1, . . . , L}) (p (ℓ) , θ (ℓ) ) = P epi h (ℓ) (y (ℓ) , ζ (ℓ) ). ( 30
)
So, the problem reduces to the lower-dimensional problem of the determination of the projection onto the convex subset epi h (ℓ) of R M (ℓ) × R for each ℓ ∈ {1, . . . , L}. By using the results in Section 2.2, these projections can be shown to have a closed form expression in a number of cases of particular interest. Two examples are now given:

• Euclidean norm functions of the form

∀y (ℓ) ∈ R M (ℓ) h (ℓ) (y (ℓ) ) = τ (ℓ) y (ℓ) (31)
where ℓ ∈ {1, . . . , L} and τ (ℓ) ∈ ]0, +∞[. Proposition 3.1. Assume that h (ℓ) is given by [START_REF] Youla | Image restoration by the method of convex projections. Part I -theory[END_REF].

Then, for every (y

(ℓ) , ζ (ℓ) ) ∈ R M (ℓ) × R, P epi h (ℓ) (y (ℓ) , ζ (ℓ) ) =      (0, 0), if y (ℓ) < -τ (ℓ) ζ (ℓ) , (y (ℓ) , ζ (ℓ) ), if y (ℓ) < ζ (ℓ) τ (ℓ) , α (ℓ) y (ℓ) , τ (ℓ) y (ℓ) , otherwise, (32)
where

α (ℓ) = 1 1 + (τ (ℓ) ) 2 1 + τ (ℓ) ζ (ℓ)
y (ℓ) .

The result of this proposition is a direct application of Proposition 2.4 and Corollary 2.7. As it will be shown in Section 4, this result is useful to deal with multivariate sparsity constraints [START_REF] Wu | Multivariate compressive sensing for image reconstruction in the wavelet domain[END_REF] or total variation bounds [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Aujol | Some first-order algorithms for total variation based image restoration[END_REF], since such constraints typically involve a sum of functions like (31) composed with linear operators corresponding to analysis transforms or gradient operators.

• Infinity norms defined as: for every ℓ ∈ {1, . . . , L} and y

(ℓ) = (y (ℓ,m) ) 1≤m≤M (ℓ) ∈ R M (ℓ) , h (ℓ) (y (ℓ) ) = max |y (ℓ,m) | τ (ℓ,m) | 1 ≤ m ≤ M (ℓ) . (33)
where (τ (ℓ,m) ) 1≤m≤M (ℓ) ∈ ]0, +∞[

M (ℓ)
. Proposition 3.2. Assume that h (ℓ) is given by [START_REF] Combettes | Inconsistent signal feasibility problems : least-squares solutions in a product space[END_REF]. Let (ν (ℓ,m) ) 1≤m≤M (ℓ) be a sequence of reals obtained by sorting (|y (ℓ,m) |/τ (ℓ,m) ) 1≤m≤M (ℓ) in ascending order, and set ν (ℓ,0) = -∞ and ν (ℓ,M (ℓ) +1) = +∞. Then, for every

ζ (ℓ) ∈ R, (p (ℓ) , θ (ℓ) ) = P epi h (ℓ) (y (ℓ) , ζ (ℓ) ) is such that p (ℓ) = (p (ℓ,m) ) 1≤m≤M (ℓ) , with p (ℓ,m) =      y (ℓ,m) , if |y (ℓ,m) | ≤ τ (ℓ,m) θ (ℓ) , τ (ℓ,m) θ (ℓ) , if y (ℓ,m) > τ (ℓ,m) θ (ℓ) , -τ (ℓ,m) θ (ℓ) , if y (ℓ,m) < -τ (ℓ,m) θ (ℓ) , (34)
θ (ℓ) = max ζ (ℓ) + M (ℓ) m=m ν (ℓ,m) (τ (ℓ,m) ) 2 , 0 1 + M (ℓ) m=m (τ (ℓ,m) ) 2 , ( 35
)
and m is the unique integer in {1, . . . , M (ℓ) + 1} such that

ν (ℓ,m-1) < ζ (ℓ) + M (ℓ) m=m ν (ℓ,m) (τ (ℓ,m) ) 2 1 + M (ℓ) m=m (τ (ℓ,m) ) 2 ≤ ν (ℓ,m) . ( 36
)
The proof of this proposition is given in Appendix F.

Note that when τ (ℓ,m) ≡ 1, function h (ℓ) in (33) reduces to the standard infinity norm • ∞ . This proposition can be thus employed to efficient deal with ℓ 1,∞ regularization which has attracted much interest recently [START_REF] Turlach | Simultaneous variable selection[END_REF][START_REF] Quattoni | An efficient projection for ℓ 1,∞ regularization[END_REF][START_REF] Chen | Recursive ℓ 1,∞ lasso[END_REF].

Experimental Results

In this section, we provide numerical examples to illustrate the usefulness of the proposed epigraphical projection method. The presented experiments focus on applications involving projections onto ℓ 1,p -balls where p ∈ {2, +∞}.

As already mentioned in Section 2.1, various algorithms can be used to solve non-smooth convex optimization problems and would potentially benefit from the proposed epigraphical projection technique. In this work, we will employ a primal-dual algorithm, namely the Monotone+Lipschitz Forward Backward Forward (M+LFBF) algorithm, which was recently proposed in [START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallel-sum type monotone operators[END_REF]. This proximal algorithm is able to address a wide class of convex optimization problems without requiring any matrix inversion and it offers a good performance and robustness to numerical errors. Its convergence is guaranteed (under weak conditions) and its structure makes it suitable for implementation on highly parallel architectures.

Degradation model

Set H = R N . Denote by x = (x (n) ) 1≤n≤N ∈ R N the signal of interest, and by z ∈ R N an observation vector such that z = Ax + b. It is assumed that A ∈ R N ×N is a linear operator, and b ∈ R N is a realization of a zero-mean white Gaussian noise. The recovery of x from the degraded observations is performed by following a variational approach which aims at solving the following problem minimize

x∈[µ,µ] N Ax -z 2 s. t. L ℓ=1 Ω ℓ B ℓ F x p ≤ η, (37)
where (µ, µ) ∈ R 2 with µ ≤ µ, η is a real positive constant, and F ∈ R K×N is the linear operator associated with an analysis transform. Furthermore, for every ℓ ∈ {1, . . . , L}, B ℓ ∈ R M (ℓ) ×K is a block-selection linear operator which selects a block of M (ℓ) data from its input vector. 3 For every ℓ ∈ {1, . . . , L}, Ω ℓ denotes an M (ℓ) × M (ℓ) diagonal matrix of real positive weights.

The term Axz 2 is the data fidelity corresponding to the minus log-likelihood of x. The bounds µ and µ allow us to take into account the value range of each component of x.

The second constraint involved in Problem (37) conveys a smoothness condition. It exploits the fact that natural signals usually exhibit a smooth spatial behaviour, except around some locations (e.g. object edges in images) where discontinuities arise. The proposed constraint is based on a weighted ℓ 1,p -norm, where p ∈ {2, +∞}. The associated block-sparsity measure extends many of the smoothness terms used in the literature. It reduces to the weighted ℓ 1 -norm criterion found in [START_REF] Candés | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] when each block reduces to a singleton (i.e. L = K, and, for every ℓ ∈ {1, . . . , L}, M (ℓ) = 1 and B ℓ y = y (ℓ) ). 4 It captures the ℓ 1,2 criteria present in [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Gilboa | Nonlocal Operators with Applications to Image Processing[END_REF][START_REF] Peyré | Group sparsity with overlapping partition functions[END_REF][START_REF] Bayram | Directional total variation[END_REF] when p = 2. It matches the ℓ 1,∞ criterion proposed in [START_REF] Quattoni | An efficient projection for ℓ 1,∞ regularization[END_REF] when p = +∞.

Let us define Λ = [B

⊤ 1 Ω 1 , . . . , B ⊤ L Ω L ] ⊤ and
C 2 = y ∈ R M L ℓ=1 y (ℓ) p ≤ η (38)
where M = M (1) + • • •+ M (L) and the same decomposition as in ( 24) is performed. Then, it can be observed that Problem (37) is a particular case of Problem 1.1 where

S = 2, R = 1, g 1 = • -z 2 , T 1 = A, H 1 = I, C 1 = [µ, µ] N , H 2 = ΛF ,
and C 2 is the above ℓ 1,p -ball.

Algorithmic solution

The main difficulty in solving Problem (37) stems from the second constraint. The point is that most of the applicable algorithms require to compute the projection onto C 2 . Specific numerical methods [START_REF] Van Den | Probing the Pareto frontier for basis pursuit solutions[END_REF][START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF][START_REF] Quattoni | An efficient projection for ℓ 1,∞ regularization[END_REF][START_REF] Fadili | Total variation projection with first order schemes[END_REF] have been developed for this purpose. The aim of this section is to propose an alternative method based on the splitting principle presented in Section 3. So doing, the resulting problem can be efficiently addressed by proximal algorithms. The two possible approaches are now detailed.

• Epigraphical method -The principle of this method is to decompose C 2 into the closed half-space defined by ( 27) and the closed convex set defined by ( 28) with h (ℓ) = • p . The advantage of this decomposition is that the projections onto V and E have closed forms. Indeed, P V is given by Proposition 2.1, while P E is given by Proposition 3.1 for p = 2 and Proposition 3.2 for p = +∞. We are then able to solve Problem ( 37) by means of the M+LFBF algorithm. The associated iterations are given in Algorithm 1.

Algorithm 1 Epigraphical method: M+LFBF for solving Problem ( 37)

Initialization       (v [0] , ν [0] ) ∈ R M × R L (x [0] , ζ [0] ) ∈ R N × R L θ = 2 A 2 + max{ ΛF , 1} ǫ ∈]0, 1 θ+1 [ For i = 0, 1, . . .                   γ i ∈ [ǫ, 1-ǫ θ ] x [i] , ζ [i] = (x [i] , ζ [i] ) -γ i 2A ⊤ (Ax [i] -z) + F ⊤ Λ ⊤ v [i] , ν [i] (p [i] , ρ [i] ) = P [µ,µ] N ( x [i] ), P V ( ζ [i] ) ( v [i] , ν [i] ) = (v [i] , ν [i] ) + γ i (ΛF x [i] , ζ [i] ) (a [i] , α [i] ) = ( v [i] , ν [i] ) -γ i P E γ -1 i v [i] , γ -1 i ν [i] ( v [i] , ν [i] ) = (a [i] , α [i] ) + γ i (ΛF p [i] , ρ [i] ) (v [i+1] , ν [i+1] ) = (v [i] -v [i] + v [i] , ν [i] -ν [i] + ν [i] ) ( x [i] , ζ [i] ) = (p [i] , ρ [i] ) -γ i 2A ⊤ (Ap [i] -z) + F ⊤ Λ ⊤ a [i] , α [t] (x [i+1] , ζ [i+1] ) = (x [i] -x [i] + x [i] , ζ [i] -ζ [i] + ζ [i] )
• Direct method -For completeness, we also consider the projection onto C 2 with the algorithm in [START_REF] Van Den | Probing the Pareto frontier for basis pursuit solutions[END_REF] when p = 2, 5 or the iterative algorithm in [START_REF] Quattoni | An efficient projection for ℓ 1,∞ regularization[END_REF] when p = +∞. 6 Then, M+LFBF can still be used to solve Problem [START_REF] Jacques | A panorama on multiscale geometric representations, intertwining spatial, directional and frequency selectivity[END_REF].

Note that, according to the general result in [68, Theorem 4.2], the sequence x [i] i∈N generated by Algorithm 1 is guaranteed to converge to a (global) minimizer of Problem (37). 5 Code available at www.cs.ubc.ca/∼mpf/spgl1 6 Code available at www.lsi.upc.edu/∼aquattoni 10

Smoothness constraint

We consider two smoothness constraints based on recent Non-Local TV measures [START_REF] Gilboa | Nonlocal Operators with Applications to Image Processing[END_REF][START_REF] Zhang | Bregmanized nonlocal regularization for deconvolution and sparse reconstruction[END_REF]. They constitute particular cases of the one considered in [START_REF] Jacques | A panorama on multiscale geometric representations, intertwining spatial, directional and frequency selectivity[END_REF] when L = N . They are described for 2D data in the following.

• ℓ 2 -Non-Local TV -This constraint has the form

N ℓ=1 n∈N ℓ ⊂W ℓ ω ℓ,n (x (ℓ) -x (n) ) 2 1/2 ≤ η, ( 39
)
where N ℓ is the neighbourhood support of ℓ, and W ℓ is the set of positions n ∈ {1, . . . , N } \ {ℓ} located into a Q × Q window centered at ℓ, where Q ∈ N is odd. This constraint is a particular case of the one considered in [START_REF] Jacques | A panorama on multiscale geometric representations, intertwining spatial, directional and frequency selectivity[END_REF] where K = (Q 2 -1)N and F is a concatenation of discrete difference operators F q1,q2 with (q 1 , q 2 ) ∈ {-(Q -1)/2, . . . , (Q -1)/2} 2 \ {(0, 0)}. More precisely, for every (q 1 , q 2 ), F q1,q2 is a 2D filter with impulse response: for every (n

1 , n 2 ) ∈ Z 2 , f (n1,n2) q1,q2 =      1, if n 1 = n 2 = 0, -1, if n 1 = q 1 and n 2 = q 2 , 0, otherwise. (40)
In addition, for every ℓ ∈ {1, . . . , N }, M (ℓ) ≤ Q 2 -1, B ℓ selects the components of F x corresponding to differences (x (ℓ) -x (n) ) n∈N ℓ , and the positive weights (ω ℓ,n ) n∈N ℓ are gathered in the diagonal matrix Ω ℓ .

• ℓ ∞ -Non-Local TV -We consider the following constraint

N ℓ=1 max n∈N ℓ ω ℓ,n |x (ℓ) -x (n) | ≤ η. ( 41
)
We proceed similarly to the previous constraint, except that the ℓ ∞ -norm is now substituted for the ℓ 2 -norm.

Note that the classical isotropic TV constraint (designated by ℓ 2 -TV in the following) constitutes a particular case of the ℓ 2 -NLTV one, where each neighbourhood N ℓ only contains the horizontal/vertical neighbouring pixels (M (ℓ) = 2) and the weights are ω ℓ,n ≡ 1. Similarly, the ℓ ∞ -TV constraint is a special case of the ℓ ∞ -NLTV one.

Weight estimation and neighbourhood choice

To set the weights, we got inspired from the Non-Local Means approach originally described in [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]. Here, for every ℓ ∈ {1, . . . , N } and n ∈ N ℓ , the weight ω ℓ,n depends on the similarity between patches built around the pixels ℓ and n of the image. Since our degradation process involves some missing data, a two-step approach has been adopted. In the first step, the ℓ 2 -TV approach is used in order to obtain an estimate x of the target image. This estimate is subsequently used in the second step to compute the weights through a self-similarity measure, yielding

ω ℓ,n = ω ℓ exp -δ -2 B ℓ F ℓ x -B n F n x 2 , (42)
where

δ ∈ R \ {0}, ω ℓ ∈ ]0, +∞[, B ℓ (resp. B n ) selects a Q × Q patch centered at position ℓ (resp. n) and F ℓ (resp. F n )
is a linear processing of the image depending on the position ℓ (resp. n). The constant ω ℓ is set so as to normalize the weights (i.e. n∈N ℓ ω ℓ,n = 1).

The measure in [START_REF] Aubert | Sur la minimisation d'une fonctionnelle non convexe, non différentiable en dimension 1[END_REF] generalizes the one proposed in [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], which corresponds to the case when F ℓ (resp. F n ) reduces to a Gaussian function with mean ℓ (resp. n). In the present work, we consider the foveated selfsimilarity measure recently introduced in [START_REF] Foi | Foveated self-similarity in nonlocal image filtering[END_REF], due to its better performance in denoising. This approach can be derived from ( 42) by setting F ℓ (resp. F n ) to a set of low-pass Gaussian filters whose variances increase as the spatial distance from the patch center ℓ (resp. n) grows.

For every ℓ ∈ {1, . . . , N }, the neighbourhood N ℓ is built according to the procedure described in [START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF]. In practice, we limit the size of the neighbourhood, so that M (ℓ) ≤ M .

Simulation results

The purpose of this section is twofold. First, the quality of images reconstructed with our variational approach is evaluated for different choices of regularization constraints and comparisons are made with a state-of-the-art method. Secondly, the convergence speed of the proposed epigraphical technique is evaluated w.r.t. the direct methods.

In the following experiments, if not specified otherwise, the degradation matrix A is a decimated convolution which consists of a 3 × 3 uniform blur followed by a decimation which randomly deletes 60% of the pixels (N = 0.4 × N). The standard deviation of the additive white Gaussian noise is equal to σ = 10. Since we deal with natural images, the data range bounds are µ = 0 and µ = 255. For the smoothness constraint, we set Q = 11, Q = 5, δ = 35 and M = 14.

Extensive tests have been carried out on several standard images of different sizes. The SNR and SSIM [START_REF] Wang | Mean squared error: love it or leave it?[END_REF] results obtained by using the various previously introduced TV-like constraints are collected in Table 1. In addition, a comparison is performed between our method and the Gradient Projection for Sparse Reconstruction (GPSR) method [START_REF] Figueiredo | Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems[END_REF], which also relies on a variational approach. The constraint bound for both methods was hand-tuned in order to achieve the best SNR values. The best results are highlighted in bold. A visual comparison is made in Figure 1, where two representative images are displayed. These results demonstrate the interest of considering non-local smoothness measures. Non-Local TV with ℓ 1,2norm indeed proves to be the most effective constraint with gains in SNR and SSIM (up to 1.82 dB and 0.042) with respect to ℓ 2 -TV, which in turn outperforms GPSR. The better performance of NLTV seems to be related to its ability to better preserve edges and thin structures present in images. In terms of computational time, GPSR is about twice faster than ℓ 2 -NLTV. Our codes were developed in MATLAB 7 , the operators F and F ⊤ being implemented in C using mex files.

In order to complete the analysis, we report in Figure 2 SNR/SSIM comparisons between ℓ 2 -NLTV and ℓ 2 -TV for different blur and noise configurations. These plots show that ℓ 2 -NLTV provides better results regardless of the degradation conditions.

In the second part of the section, we focus on the convergence speed of epigraph and direct methods. All the results refer to the culicoidae image cropped at 256 × 256 (N = 256 2 ), since a similar behaviour was observed for other images. The stopping criterion is set to: x [i+1]x [i] ≤ 10 -4 x [i] . For the ℓ 1,p -ball projectors needed by the direct method, we used the software publicly available on-line [START_REF] Van Den | Probing the Pareto frontier for basis pursuit solutions[END_REF][START_REF] Quattoni | An efficient projection for ℓ 1,∞ regularization[END_REF]. • Total Variation -Tables 2 and3 report a comparison between the direct and epigraphical methods for different values of η and for ℓ 2 -TV and ℓ ∞ -TV, respectively. For more readability, these are expressed as a multiplicative factor of the ℓ p -TV-semi-norm of the original image. It can be noticed that the parameter η influences both the quality of the results and the convergence speed. For the ℓ 1,2 -norm, the epigraphical method is 3.5 faster than the direct one. For the ℓ 1,∞ -norm, it is 65 times faster.

Figs. 3-a and 3-b show the relative error ∞] as a function of the computational time, where x [∞] denotes the solution computed after a large number of iterations (typically, 5000 iterations). The dashed line presents the results for the direct method while the solid line refers to the epigraphical one. These plots show that the epigraphical approach is faster despite it requires more iterations in order to converge. This can be explained by the computational cost of the subiterations required by the direct projections onto the ℓ 1,p -ball. • ℓ 2 -Non-Local Total Variation -Table 4 reports the results for ℓ 2 -NLTV. Different combinations of neighbourhood size Q and bound value η are considered. To set the weights, the first TV estimate is computed with η = 0.8. The best SNR improvement (1.82 dB over ℓ 2 -TV) is observed when a relatively small neighbourhood is used (Q = 11). This behaviour was already observed in [START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF]. It may be related to the bias-variance trade-off commonly encountered in Non-Local filters [START_REF] Lebrun | Secrets of image denoising cuisine[END_REF][START_REF] Kervrann | Local adaptivity to variable smoothness for exemplar-based image regularization and representation[END_REF].

x [i] -x [∞] / x [
In Figure 3-c, a plot similar to those in Figs. 3-a and3-b show the convergence profile. The epigraphical method requires about the same number of iterations as the direct one in order to converge. This results in a time reduction (1.8 times), as a single iteration of the epigraphical method is faster than one iteration of the direct method.

• ℓ ∞ -Non-Local Total Variation -Table 5 and Figure 3-d show the results obtained with the ℓ ∞ -NLTV constraint. Note that ℓ ∞ -NLTV requires more iterations than ℓ 2 -NLTV in order to converge. In what concerns the convergence speed, the epigraphical method is 19 times faster.

Conclusions

We have proposed a new epigraphical technique to solve constrained convex optimization problems with the help of proximal algorithms. In this paper (Part I), our attention has been turned to block-sparsity constraints based on weighted ℓ 1,p -norms with p ∈ {2, +∞}. The obtained results demonstrate the better performance of non-local measures in terms of image quality. Our results also show that the ℓ 1,2 -norm has to be preferred over the ℓ 1,∞ -norm for image recovery problems. However, it would be interesting to consider alternative applications of ℓ 1,∞ -norms such as regression problems [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. Furthermore, the experimental part indicates that the epigraphical method converges faster than the approach based on the direct computation of the projections via standard iterative solutions. Full implementations in C of the proposed algorithms and parallelization of our codes should even allow us to accelerate them [START_REF] Gaetano | Parallel implementations of a disparity estimation algorithm based on a proximal splitting method[END_REF]. Note that, although the considered application involves two constraint sets, the proposed approach can handle an arbitrary number of convex constraints. In a companion paper (Part II), we will show that the epigraphical approach can also be used to develop approximation methods for dealing with more general convex constraints. Comparison between epigraphical method (solid line) and direct method (dashed line):

x [i] -x [∞]
x [∞] in dB vs time A Proof of the Proposition 2.4

where the Lagrange multiplier α is such that

α(ϕ(p) -θ) = 0. ( 46
)
Since the value α = 0 is not allowable (since it would lead to p = y and θ = ζ), it can be deduced from the above equality that ϕ(p) = θ. In addition, differentiating the Lagrange functional in (45) w.r.t. ξ yields

ϕ(p) = θ = ζ + α ≥ ζ. (47)

  (a) Culicoidae. (b) Degraded. (c) Zoom. (d) GPSR, SNR: 17.03 dB. (e) Lena. (f) Degraded. (g) Zoom. (h) GPSR, SNR: 20.26 dB. (i) ℓ 2 -TV, SNR: 20.80 dB. (j) ℓ∞-TV, SNR: 20.25 dB. (k) ℓ 2 -NLTV, SNR: 22.62 dB. (l) ℓ∞-NLTV, SNR: 22.38 dB. (m) ℓ 2 -TV, SNR: 23.18 dB. (n) ℓ∞-TV, SNR: 22.77 dB. (o) ℓ 2 -NLTV, SNR: 24.18 dB.(p) ℓ∞-NLTV, SNR: 24.14 dB.

Figure 1 :

 1 Figure 1: Image restoration examples (noise parameters: blur = 3 × 3, σ = 10, decimation = 60%)

Figure 2 :

 2 Figure 2: SNR and SSIM values for ℓ 2 -NLTV (vertical axes) and ℓ 2 -TV (horizontal axes), for the culicoidae image. The plots show the results obtained for σ ∈ {5, 10, . . . , 50}, where lower SNR or SSIM values correspond to higher σ values. No decimation is applied in this experiment.

  ℓ∞-NLTV.

Figure 3 :

 3 Figure3: Comparison between epigraphical method (solid line) and direct method (dashed line):x [i] -x [∞]

For

  every (y, ζ) ∈ H × R, let (p, θ) = P epi ϕ (y, ζ). If ϕ(y) ≤ ζ, then p = y and θ = ζ = max{ϕ(p), ζ}. shows that (12) holds. Let us now consider the case when ϕ(y) > ζ. From the definition of the projection, we get (p, θ) = argmin (u,ξ)∈epi ϕ uy 2 + (ξζ) 2 . (44) From the Karush-Kuhn-Tucker theorem [97, Theorem 5.2], 8 there exists α ∈ [0, +∞[ such that (p, θ) = argmin ζ) 2 + α(ϕ(u)ξ)

Table 1 :

 1 SNR dB and SSIM results of our method and GPSR (noise parameters: blur = 3 × 3, σ = 10, 25.47 -0.823 24.70 -0.808 26.31 -0.836 25.87 -0.823 22.14 -0.734 Man 256 2 19.24 -0.725 18.96 -0.714 19.66 -0.741 19.51 -0.736 17.11 -0.629 Peppers 512 2 23.69 -0.801 23.25 -0.786 24.80 -0.829 24.45 -0.813 21.94 -0.709 Barbara 512 2 16.74 -0.653 16.64 -0.642 17.02 -0.673 16.99 -0.652 15.97 -0.562 Hill 512 2 22.18 -0.723 21.89 -0.715 22.55 -0.735 22.43 -0.733 20.21 -0.637 Culicoidae 1024 2 20.84 -0.855 20.49 -0.812 23.25 -0.885 22.57 -0.810 17.19 -0.725

	decimation = 60%)
	SNR (dB) -SSIM	N	ℓ 2 -TV	ℓ ∞ -TV	ℓ 2 -NLTV	ℓ ∞ -NLTV	gpsr
	Culicoidae	256 2	20.80 -0.855 20.25 -0.853 22.62 -0.897 22.38 -0.897	17.03 -0.738
	Lena	256 2	23.18 -0.783 22.77 -0.769 24.18 -0.812 24.14 -0.812	20.26 -0.678
	Cameraman	256 2	20.06 -0.774 19.68 -0.755 20.71 -0.801 20.17 -0.743	17.92 -0.673
	Boat	256 2	20.25 -0.739 19.74 -0.718 21.13 -0.770 20.77 -0.741	18.06 -0.649
	House	256 2

Table 2 :

 2 Results for the ℓ 2 -TV constraint and different values of η

	η	SNR (dB) -SSIM	direct	epigraph.	speed up
			# iter. sec.	# iter. sec.
	0.6	19.70 -0.839	116 5.74	146 2.23	2.58
	0.7	20.62 -0.862	132 7.14	151 2.29	3.11
	0.8	20.80 -0.855	160 9.17	171 2.59	3.54
	0.9	20.45 -0.826	195 11.95	196 2.98	4.01

Table 3 :

 3 Results for the ℓ ∞ -TV constraint and different values of η

	η	SNR (dB) -SSIM	direct	epigraph.	speed up
			# iter. sec.	# iter. sec.
	0.6	19.58 -0.839	195 373.84	230 6.21	60.22
	0.7	20.25 -0.853	206 417.71	231 6.39	65.34
	0.8	20.24 -0.835	233 486.17	245 6.86	70.82
	0.9	19.80 -0.800	263 551.03	267 7.44	74.06

Table 4 :

 4 Results for the ℓ 2 -NLTV constraint and some values of η and Q

	η	SNR (dB) -SSIM	direct	epigraph.	speed up
			# iter. sec. # iter. sec.
		Neighbourhood size: Q = 3
	0.7	22.26 -0.895	70 5.82	77 3.07	1.90
	0.8	22.41 -0.893	72 6.39	75 3.00	2.13
	0.9	22.06 -0.875	88 7.97	89 3.58	2.23
		Neighbourhood size: Q = 5
	0.7	22.44 -0.900	70 7.13	74 4.37	1.63
	0.8	22.58 -0.898	72 7.85	75 4.44	1.77
	0.9	22.25 -0.880	87 9.71	88 5.24	1.85
		Neighbourhood size: Q = 11
	0.7	22.50 -0.901	72 7.52	76 4.51	1.67
	0.8	22.61 -0.897	75 8.21	78 4.64	1.77
	0.9	22.27 -0.879	89 9.83	91 5.41	1.82

Table 5 :

 5 Results for the ℓ ∞ -NLTV constraint and some values of η and Q

	η	SNR (dB) -SSIM	direct	epigraph.	speed up
			# iter. sec.	# iter. sec.
		Neighbourhood size: Q = 5
	0.6	21.78 -0.889	229 766.21	251 25.61	29.92
	0.7	22.28 -0.894	213 653.31	226 22.87	28.56
	0.8	22.10 -0.876	207 596.90	216 21.89	27.26
		Neighbourhood size: Q = 5
	0.6	21.92 -0.894	231 914.03	252 43.32	21.10
	0.7	22.38 -0.897	219 749.32	232 39.73	18.86
	0.8	22.21 -0.880	212 673.60	224 38.05	17.70
		Neighbourhood size: Q = 11
	0.6	21.95 -0.894	236 940.38	256 43.75	21.49
	0.7	22.38 -0.897	222 760.49	234 39.84	19.09
	0.8	22.20 -0.879	216 689.07	227 38.80	17.76

When ϕ is Gâteaux-differentiable at y, ∂ϕ(y) = {∇ϕ(y)} where ∇ϕ(y) is the gradient of ϕ at y.

Note that the inequality in[START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] can be also replaced by an equality although it makes little difference in our approach.

This means that there exist distinct indices m 1 , . . . , m M (ℓ) in {1, . . . , K} such that, for every y = (y (k) ) 1≤k≤K ∈ R K , B ℓ y = (y (m j ) ) 1≤j≤M (ℓ) .

In this case, the constraint in[START_REF] Jacques | A panorama on multiscale geometric representations, intertwining spatial, directional and frequency selectivity[END_REF] reduces to a sum of absolute values.

R2011b version on an Intel Xeon CPU at 2.80 GHz and

GB of RAM.

By considering u 0 ∈ dom ϕ and ξ 0 > ϕ(u 0 ), the required qualification condition is obviously satisfied.

Hence, (p, θ) given by [START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms, Part I : Fundamentals[END_REF] is such that p = argmin u∈H ϕ(u)≥ζ uy 2 + (ϕ(u)ζ) 2 (48)

Furthermore, as ϕ(y) > ζ, we have inf u∈H ϕ(u)≤ζ uy 2 = P lev ≤ζ ϕ (y) -

where we have used the fact that P lev ≤ζ ϕ (y) belongs to the boundary of lev ≤ζ ϕ which is equal to u ∈ H ϕ(u) = ζ since ϕ is lower-semicontinuous and dom ϕ is open [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Corollary 8.38]. We have then

Altogether, ( 48) and ( 51) lead to

which is equivalent to (12) since

B Proof of the Proposition 2.5

Since (max{ϕζ, 0}) 2 is an even function, prox 1 2 (max{ϕ-ζ,0}) 2 is an odd function [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]Remark 4.1(ii)]. In the following, we thus focus on the case when y ∈ ]0, +∞[.

where, according to [98, Corollary 2.5], p ≥ This allows us to deduce that = χ 0 .

Let us now focus on the case when

yields p = y. On the other hand if y > (ζ/τ ) 1/β , as the proximity operator of a function from R to R is continuous and increasing [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF]Proposition 2.4]

) allows us to deduce that p is the unique value in [(ζ/τ ) 1/β , +∞[ satisfying [START_REF] Combettes | A block-iterative surrogate constraint splitting method for quadratic signal recovery[END_REF]. It can be concluded that, when ζ ∈]0, +∞[, [START_REF] Kowalski | Beyond the narrowband approximation: Wideband convex methods for under-determined reverberant audio source separation[END_REF] holds.

C Proof of the Proposition 2.6

Let us notice that 1 2 (max{τ

According to [25, Proposition 2.7], for every y ∈ H,

where α = prox ψ dC (y) dC (y)

. In addition, we have

E Proof of the Proposition 2.8

The function ϕ belongs to Γ 0 (R) since for every m ∈ {1, . . . , M }, v → max{τ (m) (ν (m)v), 0} is finite convex and (•) 2 is finite convex and increasing on [0, +∞[. In addition, ϕ is differentiable and it is such that, for every v ∈ R and every k ∈ {1, . . . , M + 1},

For every y ∈ R, as p = prox ϕ (y) is characterized by [START_REF] Pustelnik | Parallel algorithm and hybrid regularization for dynamic PET reconstruction[END_REF], there exists m ∈ {1, . . . , M + 1} such that ν (m-1) < p ≤ ν (m) and

This yields [START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF], and we have : [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] ⇔ ν (m-1) < p ≤ ν (m) . The uniqueness of m ∈ {1, . . . , M + 1} satisfying this inequality follows from the uniqueness of prox ϕ (y).

F Proof of the Proposition 3.2

For every (y (ℓ) , ζ (ℓ) ) ∈ R M (ℓ) × R, in order to determine P epi h (ℓ) (y (ℓ) , ζ (ℓ) ) we have to find min

. . .

For every θ (ℓ) ∈ [0, +∞[, the inner minimization is achieved when, for every j ∈ {1, . . . , M (ℓ) }, p (ℓ,m) is the projection of y (ℓ,m) onto [-τ (ℓ,m) θ (ℓ) , τ (ℓ,m) θ (ℓ) ], which is given by [START_REF] Kose | Filtered variation method for denoising and sparse signal processing[END_REF]. Then, the problem reduces to minimize

which is also equivalent to calculate prox ϕ+ι [0,+∞[ (ζ (ℓ) ), where ϕ is such that

By using now the fact that prox ϕ+ι [0,+∞[ = P [0,+∞[ •prox ϕ ([24, Proposition 12]) and by invoking Proposition 2.8, the expression of the optimal solution in (35) follows.


Berangere Suiphon
email: berangere.suiphon@g-scop.grenoble-inp.fr

Zineb Simeu-Abazi
email: zineb.simeu-abazi@g-scop.grenoble-inp.fr

Eric Gascard
email: eric.gascard@imag.fr

Implementation of a fault diagnosis method for timed discrete-event systems

Keywords: fault diagnosis, detection, isolation, characteristic times, timed event-discrete system, timed automaton

To date, fault diagnosis is more and more studied due to its essential role in equipment availability and human safety. This paper proposes a model-based method for event-discrete systems, using the characteristic times in its normal behavior. Indeed, a failure can be detected is these times are not respected. Then the failure is isolate using fault signature analysis. A failure can be distinguished from each other by a set of variables, called diagnostic variables. The choice of these particular variables is essential to diagnose efficiently a system. The diagnosis performance can also be quantified through two parameters : the detection and isolation times.

Introduction

Fault diagnosis plays an important role in equipment availability. In fact, it enables to improve performances, but also to minimize harmful consequences which can be catastrophic for equipment and human safety. This explains why industry devotes a significant part of its work to fault diagnosis. Industrial systems are more and more intricate and consequently more and more difficult to diagnose.

To date, many diagnostic methods have been implemented. They can be divided into two main types : model-free and model-based methods. Model-free method consists in analyzing monitoring variables to diagnose the system. In model-based method, a model of the system is established. This model includes the knowledge of normal and faulty system behavior. In this approach, the diagnosis can be decomposed into two steps : detection and isolation.

In the case of discrete-event systems (DES), the most used approach is the model-based method. The diagnosis is often based on logical models, like Petri Nets. In [START_REF] Hashtrudi Zad | Fault Diagnosis in Finite-State Automata and Timed Discrete-Event Systems[END_REF], inputs and outputs are monitored in order to detect and isolate the source of the fault. But quickly, it appeared that DES models are sometimes insufficient to diagnose all faults. Indeed, the order of the system steps is sometimes not affected by a fault, but the step times are. In this case, it is necessary to Steps for modeling a system as a Timed-Automaton add a new variable : time. That is why the recent diagnostic methods are based on timed Petri nets, stochastic automata [START_REF] Supavatanakul | Identification of timed discrete-event models for diagnosis[END_REF], timed automata [START_REF] Hristov | A time discrete-event approach to diagnosis of the damadics actuator benchmark[END_REF], template languages [START_REF] Pandalai | Template Languages for Fault Monitoring of Timed Discrete Event Processes[END_REF] or Semi-Markov processes.

In this paper, a model-based method is presented for timed discrete-event systems. This approach is based on using characteristic times of the system. After analysis and modeling of the system, a diagnostic automaton (or diagnoser ) is build. The model and the diagnoser are Timed Automata (TA). In this first part, the subject and the context of the paper have been introduced. The second part deals with the tools used to analyze and model a system. In a third part, the necessary steps for building the diagnoser are detailed. The results and validation of our method are exposed in the fourth part. A conclusion with different perspectives of this work is given in the end of the paper.

Analysis and modeling tools

The approach described in this paper is a model-based method. It means that the fault diagnosis goes by the differences between the real system and its model. In this part, we expose the successive steps to model the system (Figure 1). The first step is to analyze the system, first by a temporal analysis (1) in order to estimate the characteristic times of the system, and second by a functional analysis (2) using FAST diagram. It allows to model the system in its normal behavior [START_REF] Hashtrudi Zad | Fault Diagnosis in Finite-State Automata and Timed Discrete-Event Systems[END_REF]. Then, to consider the fault effects, the fault analysis (4) is made, using the FMECA. Finally, the global model (5) of the system is build.

Temporal analysis : search of characteristic times

The diagnostic method detailed in this paper is based on characteristic times of the system. They correspond often to the stage starts (Figure 2). If these times aren't respected, a fault is detected. This subsection explains how to collect these particular moments.

The first step consists in establishing the system dynamics. The objective is to know how the system behaves in the course of time. Many methods can be used : mathematical tools, like differential equations, or symbolic approaches, for example timed Petri nets or timed Markov processes. Once the method chosen, the dynamics is simulated. The simulation allows to collect the characteristic times. To make it easier to understand, you can refer to the following example.This instance will be used in the whole article.

Example : An autonomous electric vehicle

We consider an electric vehicle (Figure 3) which moves without operator, but guiding by variables (white lines on the ground, internal programming of a route, ...). An on-board computer gives instructions to the actuators (engines) according to the sensor values (accelerometer, gyroscope, laser, ...) and to the program it contains. In this example, the system route (on the right in Figure 3) is made up of five stages and two breaks (50 s). Five sensors (C1 to C5), placed on the route, enable to know if the vehicle is present in the end of a stage (1 if present, 0 else). The vehicle moves at a constant speed (5m/s). The interesting variables are the speed and the direction angle of the vehicle. In Cartesian coordinates (Figure 4), dynamics correspond to the following equations :

   x(t) = t 0 cos(α).v.dτ y(t) = t 0 sin(α).v.dτ
By simulating the dynamic equations, it is possible to have the characteristic times (Figure 5).

Functional analysis : using a FAST diagram

The functional analysis is used to determine the main functions, as well as the secondary ones. It allows to know exactly what the system have to do (and how). Among the mass of functional analysis tools, we have chosen to explain the FAST diagram (Figure 6). It consists in breaking a system down into its main functions. Then each function can be broken down into basic functions. For each one, the building solution is specified. This method is often used to described complex systems as exhaustive as possible.

Example : An autonomous electric vehicle

In the case of the autonomous electric vehicle, the analysis will be made for a quarter of the system. In fact, the vehicle can be decomposed into four identical parts, each made up of a axle-xheel block (wheel, suspension, engine), a steering block (steering jack, steering bar, bracket) and a power block (batteries). The FAST diagram of the subsystem is shown on the right in the figure 6. For each component, it is now easy to know its role and the associated function.

Modeling of the normal behavior

Thanks to the temporal and the functional analysis, the system is well-known enough to build the model of the normal behavior. We have chosen to model the system as a TA. A timed-automaton is a finite-state machine extended with a set of clocks. In this way, the automaton states change depending on variable value and on time, that is to say that the transition from a state to another can be a comparison between a sensor and a desired value and/or a comparison between the clock and an integer. Besides, the clocks can be reset during the run of the automaton. A simple example is given in figure 7 with four states and three transitions and where α, β are control variables and a, b sensors values.

In order to model the system in its normal behavior, we need to represent the operative and control parts (Figure 8). The operative part, or process, is composed of the dynamics (studied previously in 2.1) and the state of the sensors, depending on the system variables. In the "Sensors" block, the sensors values are changed according to system values. The control part gives the order to the operative part in function. The order depends on the state of the system. The transition between states depends on the sensor values and the time it occurs. Therefore, this "Control" block has to be a TA. We propose in the following example a modeling using Matlab. Matlab is a development environment, often used for numerical calculation, and which has two useful libraries for automation : Simulink and Stateflow. Thanks to them, it is possible to model finite-state machine and TA.

Example : An autonomous electric vehicle

In the case of the autonomous electric vehicle, the control variables are the speed and the angle. From this control values, the dynamics calculates the coordinates of the vehicle (system variables). Then the sensors on the route (C1 to C5) are updated. If the vehicle is present in the same coordinates that a sensor, its value is "1", else it is "0". The process is modeled using Matlab-Simulink (on the left in Figure 9). The "Process" block contains two Matlab programs for dynamics and sensor changes. For more informations, see [START_REF] Suiphon | Implémentation d'une méthode de diagnostic des systèmes à événements discrets[END_REF]. As for the control part, a Stateflow chart is build (on the right in Figure 9). It has an initialization state and then the seven states of the route (five movements and two breaks). The transition depends on the value sensors (to finish the first state and get in the second, the sensor C2 has to be "up", ...) and on time (a break is finished if the break lasts 50s). In each state, the control values for speed and angle are defined. If a simulation of the model is started up, the displayed vehicle behavior is the one in Figure 5.

Fault analysis : the FMECA approach

Now take an interest in the failures and their modelings. In order to model the behavior of a system in a fault state, the faults have to be listed and their consequences and causes clarified. The FMECA (Failure Mode, Effects and Criticality Analysis) is a powerful tool to describe and analyze the faulty behavior. It appears as a table (Table 1). This is a Component-FMECA : it means that the faults are listed for each component. For each fault, the FMECA gives the faulty component, the altered functions (using the FAST diagram), the effects on the system (locally and on the global behavior) and its probable causes. Then the Criticality level (C) is determined as the product of the occurrence (Occ), the severity (Sev) and the detectability (Det) :

C = Occ. × Sev. × Det.
A grading scale is provided for occurrence, severity and detectability. These ranking scales are ranging from 1 to 5 : for instance, if the occurrence is 1/5, that means that there is a few risk that it occurs ; if severity is 5/5, the consequences on the system are very serious ; if detectability is 1/5, it is easy to detect the failure. Therefore the higher the criticality, the more critical the fault is for the system. It is possible to add columns to the FMECA in order to specify the maintenance operations, the prevention tasks to avoid having the fault or the tests to identify the failure.

Example : An autonomous electric vehicle

In the case of the autonomous vehicle system, there are 18 failures related to the components of the subsystem and 10 related to the add of the five sensors. The table 2 shows some rows of the FMECA. Fig. 10. Modeling of the global system behavior

Modeling of the global behavior

With the knowledge of the fault behavior (thanks to the FMECA) and the base of the modeling (in the normal behavior), the final stage is to model the global system, that is to say in normal but also in faulty behavior. We are going to consider a single-fault scenario, which means that only one fault can occur on the system. In order to model the global system, the first step consists in choosing a failure. Then the second step is to consider the consequences on the basic model.

Firstly, the failure has to be chosen. Two methods can be used : a "random-time" method and a "randomfailure-and-time" method. The "random-time" method consists in choosing the failure to inject but the time it occurs is random. In the "random-failure-and-time" method, both failure and time are randomly chosen. In this case, the occurrence of a failure depends on its probability (present in component data sheets). For more information, refer to figure 10 ("Fault Choice" block) in the example and [START_REF] Chalagiraud | Diagnostic des défaillances par automates temporisés[END_REF].

Secondly, the model behavior has to change because of the injection of the fault. When a failure occurs, it can alter sensor values (if it is a fault on a sensor) or control variables. Therefore, the injection of fault in the system alters only the process in the model (Figure 10). If the fault analysis is precise enough about local failure effects, this step of modeling is quite easy.

Example : An autonomous electric vehicle

As said previously, two steps are necessary : failure choice and injection (Figure 11). First at all, for the failure choice, all we have to do is to tell how many failures can happen and their probabilities. Then, it is possible to choose the method for the simulation thanks to the switch. About the injection of failure effects in the model, the problem is broken down into two parts : effects on sensor values and on control variables. In the Matlab "Fault on sensors" block, each sensor are handled separately (Figure 12). Here, sensors are binary, consequently if a fault occurs on a sensor, it stays either in "down" state or in "up" state ; that is why the programming is made of logical operations. For failures affecting control variables, the "Fault on process variables" block altered separately each variable : here the speed and the angle (Figure 13). We consider that speed have several alteration degrees : low increase, no change, low decrease, high decrease and no speed. In the same manner, the angle can be unchanged (from the normal route), low deflected, high deflected and "straight ahead" (unchanged while it should). In the program, depending on the fault number, the variables have a definite alteration degree ; accordingly they take randomly values on a definite set. Finally, the global model of the system has been built successfully : the system is modeled in its normal behavior but also when a failure occurs. Now we will use this model in order to diagnose failures in the system. With only the system behavior and its characteristic times, the objective is to detect and isolate the fault.

Diagnosis by using of signature analysis

In this part, detection and isolation are first explained. Then the notion of fault signature is introduced, which is essential in the diagnostic method in this paper. Finally, the diagnoser building is explained for the example of the autonomous electric vehicle.

Detection principle

When a failure occurs in the system, the first step is to detect it. It is possible by comparing the system with its model in normal behavior for a same control. If a difference between the two behaviors is noticed, a failure is detected (Figure 14).

Isolation principle

Once the failure is detected, it need to be isolated (Figure 15). Firstly, the faulty component has to be found : it is the localization. Then a component can be associated with several failures. So all these failures are studied to find the one, which occurs in the system : it is the identification.

Fault signature analysis

In order to isolate the failure, the method described in this paper uses the notion of fault signature. It considers that a failure can be distinguished from each other by the values of some system variables. The definition of fault signature in this paper is based on [START_REF] Aseervatham | Diagnostic à base de modèles intégrant les modèles de fautes[END_REF]. More precisely, a fault affects some variable values in a definite manner. It can be represented as a diagnostic matrix.

On figure 16, you can see an example of diagnostic matrices. M 1 is the diagnostic matrix for a variable x 1 which can take three values V 1 to V 3 , in a system with five failures d 1 to d 5 . If x 1 can take the value V j in presence of the failure d i , the matrix element m ij = 1 (i th line, j th column), else m ij = 0. Then, the signature of the failure d i corresponds to the i th row of the matrix.

Moreover, it is interesting to use several system variables in order to refined the diagnosis. Indeed, in the example on figure 16, if x 1 = V 2 , the possible failures are d 1 and d 3 . But to distinguish the two failures, it is necessary to add a new variable which does not take the same value for both failures. By adding the variable x 2 and its In conclusion, the system variables used for the diagnosis have to be well-chosen in order to isolate correctly all failures.

Example of diagnoser building

In this subsection, the method to build a diagnoser using fault signatures is introduced. The approach is directly applied to the example of the autonomous electric vehicle, programed with Matlab.

The diagnoser considers time in order to detect failures, consequently the detection has to be a TA (Figure 17). After detection, the failure is isolated thanks to diagnostic matrices and the variable values.

Let's take an interest in the detection (Figure 18). In this figure, the program is broken down into six parts. The first block is the control model of the system. In the second part, the cycle of the system is represented again but if the characteristic times are not respected, the TA leaves the cycle to enter in a detection state. At this moment, the detection occurs.

Once the fault is detected, the system variable values (angle and speed) are calculated (3 rd block in figure 18). Indeed, we consider that the system variables used for isolation are unobservable, but can be calculated using observable variables. Once the calculus made, the values are sent out of the TA. Moreover, the TA contains three blocks to save the injection, detection and isolation times. They will serve to validate the diagnostic method (Section 3). Now, let's take an interest in the isolation (Figure 17). The isolation block has as entries the diagnostic matrices and the variable values from the TA. For each variable, the column of its diagnostic matrix which corresponds in the calculated value is recovered. By applying a logical "AND" operation between these column vectors, the result of the diagnosis is returned. In the best case, the result column vector contains only one "1" and the failure is isolated. Else, the result is a set of possible failures including the one which occurred.

Results and validation

The diagnostic method explained in this paper has been validated. Indeed, for each failure injection, the diagnoser is able to detect and isolate it, as you can see in the table 3. For some failures, the diagnoser is unable to distinguish failures and return a set of failure instead of only one. That means that the variables used for the diagnosis are not enough precise and another one would be necessary to distinguish the failures. However, the set of failures in result contains always the fault which occurred and it targets a precise block of the system : in the case of the 15 th failure (i.e. an out-of-order encoder for the steering), all the failures proposed in the diagnoser result belongs to the steering block. The diagnoser allows nevertheless to locate the failure and helps the operate for the maintenance. About the speed of the diagnostic, in this method, it depends on the injection time. If the injection occurs at the beginning of a new state, the detection can be a bit long : the failure is detected when the state should change in the normal behavior or in the worst case (for the sensor faults) after an entire cycle. However, the isolation is quite fast. It depends essentially on the calculation of variables values.

Conclusion

The diagnostic method presented in this paper is interesting, because it can be applied on all timed eventdiscrete systems easily. Some improvements can nevertheless be brought about : the choice of system variables used for the diagnosis could be studied in order to avoid the uncertainty of the diagnoser. By adding some sensors in precise place on the system, it could be possible to add diagnostic variables and refine the results. In every case, using the characteristic times of a system is a good idea to diagnose failures. It is an interesting way of work.

Fig. 1 .

 1 Fig. 1. Steps for modeling a system as a Timed-Automaton

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. Example of characteristic times (T1 to T4)

Fig. 5 .

 5 Fig. 5. Characteristic times of the vehicle on the route

Fig. 8 .

 8 Fig. 8. Modeling of the normal behavior using a TA

Fig. 11 .

 11 Fig. 11. Fault choice and injection for the autonomous electric vehicle

Fig. 14 .

 14 Fig. 14. Failure detection

Fig. 16 .

 16 Fig. 16. Example of diagnostic matrices for variables x1 and x2

Fig. 18 .

 18 Fig. 18. The detection as a TA 1: control model -2: detection -3: variable calculation -4,5,6: injection, detection and isolation times

Table 1

 1 Structure of a FMECA

	Failure	Item	Failure Mode	Function		Causes	Effects	Occ.	Sev.	Det.	Criticality
	Mode Id.
	1	Component block	Failure of the component	Altered tion(s)	func-	Causes problem, ...) (extern aggression, intern	Effects (in the block, the global system) in the sub-system, in	*/5	*/5	*/5	(product)

Table 2

 2 Some rows of the FMECA of the autonomous electric vehicle

	Failure	Item	Failure Mode	Function		Causes	Effects	Occ.	Sev.	Det.	Criticality
	Mode Id.
							-Transmitted mechan-
							ical energy does not
	1	Axle-Wheel block	Damaged wheel	Transform tric energy mechanical ergy	elec-in en-	nails, pebble, tire wear, ...	correspond to request -Reduced speed -Altered wheel motion tion -Altered vehicle direc-	5	3	1	15
							-Important energetic
							consumption
						Defective suspen-
	6	Axle-Wheel block	Out-of-order suspension	Store energy		sion jack, dfail-lant, broken sus-pension jack ele-	No store of the load energy	3	5	2	30
	11 17	Steering block Power	Crooked steering bar Empty bat-tery	Convey cycle energy turning System power	ment Too high pressure on the bar Forgetting load, manufacturing fault, wear	Unsufficient cycle energy No power in the sys-turning tem, no motion	2 5	4 3	4 1	32 15
	18	Power	No power	System power	Short circuit	No power in the sys-tem, no motion	4	3	3	36

Table 3

 3 Results of the diagnosis method : detection and isolation times and diagnoser results for several fault injections

	Injected failure number	Injection time	Detection time	Isolation time	Diagnostic time	Diagnoser result (failure number)
	1	197.5	250.5	254	56.5	1
	15	139.5	250.5	254	114.5	9, 11, 13, 14 or 15