Jean-Franois Aujol

Guy Gilboa
email: guy.gilboa@ee.technion.ac.il

Nicolas Papadakis
email: nicolas.papadakis@math.u-bordeaux.fr

Theoretical Analysis of Flows Estimating Eigenfunctions of One-homogeneous Functionals for Segmentation and Clustering

Nonlinear eigenfunctions, induced by subgradients of one-homogeneous functionals (such as the 1-Laplacian), have shown to be instrumental in segmentation, clustering and image decomposition. We present a class of flows for finding such eigenfunctions, generalizing a method recently suggested by Nossek-Gilboa. We analyze the flows on grids and graphs in the time-continuous and time-discrete settings. For a specific type of flow within this class, we prove convergence of the numerical iterations procedure and prove existence and uniqueness of the time-continuous case. Several examples are provided showing how such flows can be used on images and graphs.

Introduction

Eigenvalue analysis of linear operators is by now very well understood theoretically and has shown to be an essential framework for the analysis and understanding of many scientific and engineering problems. Consequently, a vast research was devoted to numerically solve eigenvalue problems [START_REF] Saad | Iterative methods for sparse linear systems[END_REF][START_REF] Demmel | Applied numerical linear algebra[END_REF]. In recent years, there is a growing interest in nonlinear eigenvalue problems, which are based on nonlinear operators. Such problems appear in image processing [START_REF] Benning | Ground states and singular vectors of convex variational regularization methods[END_REF][START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF][START_REF] Moeller | Learning nonlinear spectral filters for color image reconstruction[END_REF], computer vision [START_REF] Zeune | Multiscale segmentation via bregman distances and nonlinear spectral analysis[END_REF], classification and learning [START_REF] Bresson | Multiclass total variation clustering[END_REF][START_REF] Merkurjev | Global binary optimization on graphs for classification of high-dimensional data[END_REF][START_REF] Hein | An inverse power method for nonlinear eigenproblems with applications in 1-spectral clustering and sparse pca[END_REF]. In these problems the nonlinear operators are derived from norms, semi-norms or in general one-homogeneous functionals, where the operator is essentially a subgradient element. In this paper we present a class of flows that converge to nonlinear eigenfunctions of one-homogeneous functionals.

We are interested in solving the following nonlinear eigenvalue problem

λu ∈ ∂J(u), (1)
where J is a convex one-homogeneous functional, ∂J(u) is the subdifferential and λ denotes the eigenvalue. We refer to u admitting (1) as an eigenfunction of J. More details and precise definitions are given in the following section.

A thorough investigation of such eigenfunctions was conducted for the case of the total-variation (TV) functional in the continuous setting. Meyer already observed in [START_REF] Meyer | Oscillating patterns in image processing and in some nonlinear evolution equations[END_REF] that for the ROF problem [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] (T V -L 2 square) for the case of a disk, the solution is the same disk with reduced contrast. In a series of studies [START_REF] Andreu | Minimizing total variation flow[END_REF][START_REF] Andreu | Some qualitative properties for the total variation flow[END_REF][START_REF] Bellettini | The total variation flow in R N[END_REF] shapes which preserve their shape under the TV gradient flow were characterized (termed calibrable sets). It was shown that convex characteristic sets in R 2 with a certain bound on their curvature, are all eigenfunctions of TV. It was realized in a more general manner (see e.g. [START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF]) that eigenfunctions of one-homogeneous functional preserve their shape under three convex regularization methods -gradient flow, minimization with L 2 square and inverse-scale-space [START_REF] Burger | Nonlinear inverse scale space methods[END_REF] (the time continuous form of Bregman iterations [START_REF] Osher | An iterative regularization method for total variation based image restoration[END_REF]). Thus one can view eigenfunctions essentially as atoms of the regularizer, having spatial features which are well preserved in the regularization procedure (up to some contrast change).

The above insights lead to attempts to decompose signals and images into distinct components based on eigenvalue analysis [START_REF] Benning | Ground states and singular vectors of convex variational regularization methods[END_REF][START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF][START_REF] Schmidt | Inverse scale space decomposition[END_REF]. For the gradient flow with respect to one-homogeneous functionals, eigenfunctions decay linearly with respect to the time (flow) parameter and disappear at a finite time point. Thus taking the second time derivative of the solution of the flow yields a single response in time. This characteristic behavior was used to formulate a decomposition technique based on TV, referred to as spectral-TV decomposition [START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF], where certain nonlinear TV filters were defined in an analog manner to Fourier analysis. It was shown how one can extract desired features (and in particular eigenfunctions) in a range of scales (corresponding to eigenvalues) with high accuracy and with full contrast preservation. The method was later generalized to one-homogeneous functionals in [START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF] where certain properties, like orthogonality of the decomposed components, were shown in specific settings. Applications related to denoising [START_REF] Moeller | Learning nonlinear spectral filters for color image reconstruction[END_REF], texture manipulation [START_REF] Horesh | Separation surfaces in the spectral tv domain for texture decomposition[END_REF][START_REF] Benning | Nonlinear spectral image fusion[END_REF] and segmentation of medical data [START_REF] Zeune | Multiscale segmentation via bregman distances and nonlinear spectral analysis[END_REF] were suggested.

In [START_REF] Benning | Ground states and singular vectors of convex variational regularization methods[END_REF] nonlinear eigenfunctions for inverse problems (termed ground states) were investigated, the respective generalized Rayleigh quotient were analyzed and analytic examples of anisotropic TV were shown. Eigenfunctions related to the total-generalized-variation (TGV) functional [START_REF] Bredies | Total generalized variation[END_REF] and to infimal convolution TV [START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF] were investigated in [START_REF] Müller | Advanced image reconstruction and denoising: Bregmanized (higher order) total variation and application in pet[END_REF][START_REF] Benning | Higher-order tv methods:enhancement via bregman iteration[END_REF] and properties of particular eigenfunctions of TGV were shown theoretically and numerically. Examples of certain eigenfunctions for different extensions of TV to color images were presented in [START_REF] Duran | Collaborative total variation: a general framework for vectorial tv models[END_REF].

In the field of machine learning it was shown [START_REF] Bresson | Total variation, cheeger cuts[END_REF][START_REF] Bühler | Spectral clustering based on the graph plaplacian[END_REF] that the Cheeger cut problem can be solved by solutions of the 1-Laplacian eigenvalue problem and consequently by minimizations of the total-variation functional on graphs. This was later developed in several studies for classification, clustering and segmentation in the binary-and multiple-class case [START_REF] Bresson | Convergence and energy landscape for cheeger cut clustering[END_REF][START_REF] Bresson | Multiclass total variation clustering[END_REF][START_REF] Merkurjev | Global binary optimization on graphs for classification of high-dimensional data[END_REF][START_REF] Hein | An inverse power method for nonlinear eigenproblems with applications in 1-spectral clustering and sparse pca[END_REF]. Solutions of the Cheeger problem by using projections was shown in [START_REF] Carlier | Approximation of maximal cheeger sets by projection[END_REF]. Uniqueness and regularity of Cheeger sets in R N were analyzed in [START_REF] Caselles | Some remarks on uniqueness and regularity of cheeger sets[END_REF]. A flow, based on the MBO scheme [START_REF] Merriman | Diffusion generated motion by mean curvature[END_REF], to refine graph-Laplacian eigenvectors for classification based on a diffuse interface model was proposed in [START_REF] Merkurjev | An mbo scheme on graphs for classification and image processing[END_REF]. An algorithm to construct particular TV eigenfunctions on graphs with certain regularity, referred to as nonlocal disks, was shown in [START_REF] Aujol | Fundamentals of non-local total variation spectral theory[END_REF].

In this work we present a family of new nonlinear flows, which considerably generalize the initial work of [START_REF] Nossek | Flows generating nonlinear eigenfunctions[END_REF]. Moreover for a specific type of flow a comprehensive theoretical analysis is provided. Our proposed flows are very general, and can be evolved on both graphs and grids to solve various eigenvalue problems.

Main contributions

The main contributions of this paper are as follows:

1. We first analyze the flow of [START_REF] Nossek | Flows generating nonlinear eigenfunctions[END_REF]. Then a generalized α-flow is proposed for finding eigenfunctions. It is based on different normalizations between the function and its subgradient. A thorough analysis is presented along with a time discrete formulation of iterative convex optimizations to realize the flow.

2. For the specific case of α = 1 we are able to present a complete theory of the flow, including proof of convergence of the discrete case and existence and uniqueness of the time-continuous case.

The plan of the paper is the following. We first introduce some basic material for one homgeneous functionals in Section 2. We then analyse the flow of [START_REF] Nossek | Flows generating nonlinear eigenfunctions[END_REF] in Section 3. We introduce a generalized α-flow for finding eigenfunctions in Section 3.2. Section 4 is devoted to the particular case when α = 1 in the previous flow. For this specific choice of α, we are able to prove existence and uniqueness of a solution, as well as the convergence to the solution of a numerical scheme. In Section 5, we illustrate our theoretical analysis with some numerical examples.

One homogeneous functionals

In this section, we outline some basic properties for one homogeneous functionals.

Introduction

We consider an absolutely one homogeneous functional J that takes as input a function u : x ∈ Ω → R defined on a domain Ω ⊂ R 2 . Ω can either be a discrete domain of size |Ω| = N or an open convex bounded set with Lipschitz boundary. u are elements of some Hilbert space X (e.g. X can be L 2 (Ω)) embeded with some inner product . , . J : X → R {+∞} is assumed to be proper, convex and lower semi-continuous (lsc). Absolutely one-homogeneous functionals satisfy

J(cu) = |c|J(u), ∀c ∈ R, ∀u ∈ X. (2)
The functional J in finite dimensions can be, for instance, of the general form:

J(u) = N i=1   N j=1 w ij |u i -u j | q   1/q , (3)
for q ≥ 1, with w ij ≥ 0 (usually symmetric weights are assumed w ij = w ji ). This formulation can be understood as a typical one-homogeneous functional on weighted graphs. In this case u i is the value of the function u at node i on the graph and w ij is the weight between node i and node j. As grids of any dimension can be realized by specific graph structures, this formulation applies to standard grids as well. Thus (3), with appropriate weights, can be the spatial discrete version of anisotropic TV (q = 1), isotropic TV (q = 2) and anisotropic or isotropic nonlocal TV. We recall the subgradient definition for general convex functionals

p ∈ ∂J(u) ⇔ J(v) -J(u) ≥ p, v -u , ∀v.
We also note the relation to the convex conjugate J * J(u) = sup p u, p -J * (p).

Below we state some properties of one-homogeneous functionals.

Property 1. A function J defined in (3) admits:

(a) If p ∈ ∂J(u), then J(u) = p, u , (b) If p ∈ ∂J(u), then J(v) ≥ p, v , ∀v.
Notice in particular that from (b) we get that ∂J(u) ⊂ ∂J(0) ∀u ∈ X.

Property 2. The convex conjugate J * of a one-homogeneous functional is the characteristic function of the convex set {∂J(0)}. Moreover, when Ω is included in a finite dimensional space, we have [START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF]:

∃C > 0 s.t. ||p|| 2 ≤ C, ∀p ∈ ∂J(0). (4)
From the equivalence of norms, we have that if u is of zero mean, there exists a constant κ > 0 for which ||u|| 2 ≤ κJ(u), ∀u such that u, 1 = 0.

(

) 5
The nullspace of the functional is defined by

N (J) = {u ∈ X | J(u) = 0} . (6)
The properties below are shown in [START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF].

Property 3. An absolutely one-homogeneous functional J is a seminorm and its nullspace is a linear subspace.

Property 4. If a unit constant function u = 1 is in N (J) then any subgra- dient p admits p, 1 = 0.
We use 2 and 1 norms of u defined as ||u|| 2 = u, u and ||u|| 1 = u, sign(u) .

Eigenfunctions of J

In this work, we are interested in the eigenfunctions of functionals J that are defined as follows.

Definition 1 (Eigenfunction of J). An eigenfunction of J is a function that satisfies the eigenvalue problem (1), so that J(u) = λu, u = λ||u|| 2 2 and λ = J(u)

||u|| 2 2 ≥ 0.
An interesting insight on the eigenvalue λ can be gained by the following proposition. We define K = {∂J(0)} to be the set of possible subgradients for any u. Indeed if p ∈ ∂J(u) then p ∈ ∂J(0). We first note that an eigenfunction that admits λu ∈ ∂J(u) has zero mean from Property 4 above. Next, as illustrated in Figure 1 we have the following result.

Proposition 1. For any non constant eigenfunction u, we have ∀µ ≥ λ,

λu = Proj K (µu), where Proj K is the orthogonal projection onto K = {∂J(0)} Proof. If u is a non constant eigenfunction, λu is on the boundary of K. As K is bounded (||p|| 2 ≤ C, ∀p ∈ K) then for all µ > λ, µu / ∈ K.
Let us denote as v the orthogonal projection of µu onto K. For all w ∈ K and w = v, v satisfies:

1 2 ||v -µu|| 2 2 < 1 2 ||w -µu|| 2 2 1 2 ||v|| 2 2 -2µ v, u < 1 2 ||w|| 2 2 -µ w, u
In particular, if we assume by contradiction that v = λu then:

1 2 ||v|| 2 2 -µ v, u < 1 2 ||λu|| 2 2 -µ λu, u 1 2 ||v|| 2 2 -µ v, u < 1 2 ||λu|| 2 2 -µJ(u) 1 2 ||v|| 2 2 < 1 2 ||λu|| 2 2 since J(u) ≥ v, u . We thus have ||v|| 2 < ||λu|| 2 which yields ||v|| 2 ||u|| 2 < λ. We denote ṽ = ||v|| 2 ||u|| 2 u and observe that ||µu -ṽ|| 2 2 = µ 2 ||u|| 2 2 + ||v|| 2 2 -2µ||u|| 2 ||v|| 2 ≤ µ 2 ||u|| 2 2 + ||v|| 2 2 -2µ u, v = ||µu -v|| 2 2 , (7)
and

||µu -ṽ|| 2 = µ - ||v|| 2 ||u|| 2 ||u|| 2 > (µ -λ) ||u|| 2 = ||µu -λu|| 2 (8)
From ( 7) and ( 8) we get ||µu-λu|| 2 < ||µu-v|| 2 so v can not be the orthogonal projection of µu onto K.

A flow for finding eigenfunctions of J

In this section, following the method introduced in [START_REF] Nossek | Flows generating nonlinear eigenfunctions[END_REF], we study flows for estimating eigenfunctions of one-homogeneous functions J satifying Property 4.

Introduction

In order to find eigenfunctions of J, Nossek and Gilboa have introduced the flow [START_REF] Nossek | Flows generating nonlinear eigenfunctions[END_REF]:

u(0) = u 0 , u t = u ||u|| 2 -p ||p|| 2 , p ∈ ∂J(u). (9)
Proposition 2. Assume that there exists a solution u of the flow [START_REF] Benning | Ground states and singular vectors of convex variational regularization methods[END_REF]. Then the following property holds:

d dt 1 2 ||u(t)|| 2 2 ≥ 0 (10)
Moreover, we have:

||u(t)|| 2 ≤ ||u 0 || 2 + t (11)
Proof. Recalling that p, u ≤ ||p|| 2 ||u|| 2 , this flow ensures that:

d dt 1 2 ||u(t)|| 2 2 = u, u t = u, u ||u|| 2 - p ||p|| 2 = ||u|| 2 - u, p ||p|| 2 ≥ 0
We can also remark that d dt

1 2 ||u(t)|| 2 2 ≤ ||u(t)|| 2 so that ||u(t)|| 2 ≤ ||u 0 || 2 + t.
Finally, if u 0 is of zero mean, Property 4 ensures that u(t) is of zero mean, for all t > 0.

Proposition 3. Assume that there exists a solution u of the flow [START_REF] Benning | Ground states and singular vectors of convex variational regularization methods[END_REF]. Then the following property holds:

d dt J(u(t)) ≤ 0 for almost every t. (12)
Moreover, t → J(u(t)) is non increasing for all t ≥ 0.

Proof. We make use of Lemma 3.3 page 73 in [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF] which states that t → J(u(t)) is an absolutely continuous function (see also Lemma 4.1 in [START_REF] Apidopoulos | On a second order differential inclusion modeling the FISTA algorithm[END_REF]). Moreover, recalling that p, u ≤ ||p|| 2 ||u|| 2 , this flow ensures that we have for almost every t (using again Lemma 3.3 of [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]):

d dt J(u(t)) = p, u t = p, u ||u|| 2 - p ||p|| 2 = u, p ||u|| 2 -||p|| 2 ≤ 0.
This inequality holds for almost every t, and since t → J(u(t)) is an absolutely continuous function, we deduce that it is a non increasing function.

The PDE (9) converges iff u t = 0 so that

p = ||p|| 2 ||u|| 2 u ∈ ∂J(u) ⇒ p = J(u) ||u|| 2 2 u
and u is an eigenfunction of J with eigenvalue λ = J(u)

||u|| 2 2 .

Generalized flow

Let us now see the previous flow (9) as a specific instance of a more general framework. We define a flow for α ∈ [0; 1] as:

u(0) = u 0 , u t = J(u) ||u|| 2 2 α u -J(u) ||p|| 2 2 1-α p, p ∈ ∂J(u). (13)
Notice that for α = 1/2, we retrieve the flow of Nossek and Gilboa [START_REF] Benning | Ground states and singular vectors of convex variational regularization methods[END_REF], up to a normalization with J 1/2 (u).

Proposition 4. For u 0 of zero mean and ∀α ∈ [0; 1], the trajectory u(t) of the PDE (13) satisfies the following properties:

(i) u(t), 1 = 0. (ii) d dt J(u(t)) ≤ 0 for almost every t. Moreover, t → J(u(t)) is non in- creasing. If α = 0, we have for almost every t that d dt J(u(t)) = 0 and t → J(u(t)) is constant. (iii) d dt ||u(t)|| 2 ≥ 0 and d dt ||u(t)|| 2 = 0 for α = 1.
(iv) If the flow converge to u * , we have

p * = J 2α-1 (u * ) ||p * || 2(1-α) 2 ||u * || 2α 2 u * ∈ ∂J(u * ) so that u * is an eigenfunction.
Proof. Property (iii) is obtained as follows:

d dt 1 2 ||u(t)|| 2 2 = u, u t = u, J(u) ||u|| 2 2 α u - J(u) ||p|| 2 2 1-α p = J α (u) ||u|| 2-2α 2 - J 2-2α (u) ||p|| 2-2α 2 ≥ 0.
For property (ii), we use once again Lemma 3.3 of [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]. For almost every t, it holds:

d dt J(u(t)) = p, u t = p, J(u) ||u|| 2 2 α u - J(u) ||p|| 2 2 1-α p = J 1-α (u) J 2α (u) ||u|| 2α 2 -||p|| 2α 2 ≤ 0 Since t → J(u(t)
) is absolutely continuous (thanks to Lemma 3.3 of [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]), we deduce that it is non increasing.

Properties of a semi-explicit scheme

We can look at the following semi implicit numerical scheme:

u k+1 -u k δt = J(u k ) ||u k || 2 2 α u k+1 - J(u k ) ||p k || 2 2 1-α p k+1 (14)
It is easier to analyse this scheme than the previous continuous equation.

Moreover, the properties that we prove on this scheme will be usefull in the next section.

Proposition 5. For u 0 of zero mean and δt such that

1 δt > J(u k ) ||u k || 2 2 α
, then the sequence (u k ) is defined for all k ≥ 0, and the trajectory u k given by the numerical scheme (14) satisfies:

1 u k , 1 = 0. 2 J(u k+1 ) ||u k+1 || 2 ≤ J(u k ) ||u k || 2 , 3 ||u k+1 || 2 2 ≥ u k+1 , u k ≥ ||u k || 2 2 4 ∀p k ∈ ∂J(u k ), ||p k+1 || 2 2 ≤ p k+1 , p k ≤ ||p k || 2
2 and p k+1 , u k ≥ 0 Proof. Let us rewrite the scheme [START_REF] Bresson | Total variation, cheeger cuts[END_REF] as

u k+1 -u k δt = β k u k+1 -γ k p k+1 (15)
where

β k = J(u k ) ||u k || 2 2 α and γ k = J(u k ) ||p k || 2 2 1-α
for the sake of clarity. We define

F (u, u k ) = 1 2γ k δt u -u k 2 2 - β k 2γ k u 2 2 + J(u), (16)
u k+1 is the unique minimizer of F (., u k ), as soon as 1 δt > β k , i.e. 1-δtβ k > 0. 1 Let us underline that if u 0 is of zero mean, since p is always of zeros mean, then u k is also of zero mean, so that property (i) of Proposition 4 is satisfied numerically.

2 If ||u k || 2 = 0, then u k+1 = u k = 0. Otherwise if ||u k || 2 > 0, we have: F (u k+1 , u k ) ≤ F ( ||u k+1 || 2 ||u k || 2 u k , u k ). Hence: u k+1 -u k 2 2 2γ k δt - β k 2γ k u k+1 2 2 + J(u k+1 ) ≤ ||u k+1 || 2 ||u k || 2 u k -u k 2 2γ k δt - β k 2γ k u k+1 2 2 + ||u k+1 || 2 ||u k || 2 J(u k ) 1 2γ k δt u k+1 -u k 2 2 + J(u k+1 ) ≤ 1 2γ k δt (||u k+1 || 2 -||u k || 2 ) 2 + ||u k+1 || 2 ||u k || 2 J(u k ). As u k+1 -u k 2 2 = ||u k+1 || 2 2 + ||u k || 2 2 -2 u k , u k+1 ≥ ||u k+1 || 2 2 + ||u k || 2 2 - 2||u k+1 || 2 .||u k || 2 = (||u k+1 || 2 -||u k || 2 ) 2 then we deduce that J(u k+1 ) ≤ ||u k+1 || 2 ||u k || 2 J(u k ). ( 17
)
3 We assume that 1 -δtβ k > 0. First notice that from [START_REF] Bresson | Multi-class transductive learning based on 1 relaxations of cheeger cut and mumford-shahpotts model[END_REF]:

u k+1 (1 -δtβ k ) = u k -γ k p k+1 u k+1 = 1 1 -δtβ k (u k -γ k δtp k+1 ) (18)
Next, as u k+1 = u k + δt(β k u k+1 -γ k δtp k+1 ), then:

||u k+1 || 2 2 = ||u k || 2 2 + 2δt u k , β k u k+1 -γ k p k+1 + (δt) 2 ||β k u k+1 -γ k p k+1 || 2 2 ≥ ||u k || 2 2 + 2δt u k , β k 1 -δtβ k (u k -γ k δtp k+1 ) -γ k p k+1 ≥ ||u k || 2 2 + 2δt β k 1 -δtβ k ||u k || 2 2 -γ k β k δt 1 -δtβ k + 1 u k , p k+1 ≥ ||u k || 2 2 + 2δt 1 -δtβ k β k ||u k || 2 2 -γ k J(u k ) . (19)
We now recall that

β k = J(u k ) ||u k || 2 2 α and γ k = J(u k ) ||p k || 2 2 1-α
, hence:

β k ||u k || 2 2 -γ k J(u k ) = J(u k ) ||u k || 2 2 α ||u k || 2 2 - J(u k ) ||p k || 2 2 1-α J(u k ) = (J(u k )) α ||u k || 2-2α 2 - (J(u k )) 2-2α ||p|| 2-2α 2 ≥ 0, (20)
since 19) and ( 20), we get

J(u k ) ≤ ||u k || 2 .||p k || 2 . From (
||u k+1 || 2 ≥ ||u k || 2 . ( 21
)
Notice that we can also deduce from relations ( 15), ( 19) and ( 20) that

||u k+1 || 2 2 ≥ ||u k || 2 2 + (δt) 2 ||β k u k+1 -γ k p k+1 || 2 2 ||u k+1 || 2 2 ≥ ||u k || 2 2 + ||u k+1 -u k || 2 2 2 u k+1 , u k ≥ 2||u k || 2 2 (22) so that u k+1 , u k ≤ ||u k+1 || 2 ||u k || 2 ≤ ||u k+1 || 2 2 .
4 The optimality conditions of the minimizer of ( 16) state that there exists p k+1 ∈ ∂J(u k+1 ) such that 1 γ k δt (u k+1 -u k ) -β k γ k u k+1 + p k+1 = 0, which gives

p k+1 = 1 γ k δt u k - 1 γ k 1 δt -β k u k+1 := µu k -νu k+1 , (23)
with µ ≥ ν ≥ 0. Taking the scalar product of ( 23) with p k+1 , we have:

p k+1 , p k+1 = µ u k , p k+1 -ν u k+1 , p k+1 ||p k+1 || 2 2 = µ u k , p k+1 -νJ(u k+1 ) ||p k+1 || 2 2 + νJ(u k+1 ) ≤ µJ(u k ), (24)
where we observe that u k , p k+1 ≥ 0. Next, by taking the scalar product of ( 23) with any p k ∈ ∂J(u k ), we have

p k+1 , p k = µ u k , p k -ν u k+1 , p k µJ(u k ) ≤ p k+1 , p k + ν J(u k+1 ). (25)
By ( 24) and ( 25) we get

||p k+1 || 2 2 ≤ p k+1 , p k so that ||p k+1 || 2 ≤ ||p k || 2 . Corollary 1. If 1 δt > J(u 0 ) ||u 0 || 2 2 α , then the assumption 1 δt > J(u k ) ||u k || 2 2 α
of Proposition 5 is valid ∀k ≥ 0.

Proof. Let us assume that 1 δt > J(u 0 )

||u 0 || 2 2 α
. Make the induction hypothesis

that 1 δt > J(u k ) ||u k || 2 2 α
, for all k ≤ N . Then to prove that it still holds for N + 1, one just need to notice that from relations ( 17) and ( 21), we have

J(u N +1 ) ||u N +1 || 2 ≤ J(u 0 ) ||u 0 || 2 J(u N +1 ) ||u N +1 || 2 2 ≤ J(u 0 ) ||u 0 || 2 2 , ( 26
) so that J(u N +1 ) ||u N +1 || 2 2 α ≤ J(u 0 ) ||u 0 || 2 2 α < 1 δt .
4

The case α = 1

From now, we will assume that ||u 0 || 2 = 1 and we will restrict our attention to the case when α = 1, where the flow ( 13) is:

u(0) = u 0 , u t = J(u) ||u|| 2 2 u -p, p ∈ ∂J(u). ( 27
)
This flow may be easier to analyze since we get rid of ||p|| 2 , while keeping constant ||u(t)|| 2 . Observing that J(u) = p, u so that u t = p, u ||u|| 2 u ||u|| 2 -p, the behaviour of this flow is illustrated in Figure 2. The PDE makes u evolve on the boundary of an 2 ball of radius ||x 0 || 2 (we assume that ||x 0 || 2 > C defined in (4)) until there exists a subgradient p ∈ ∂J(u) ⊂ K such that p is the orthogonal projection of u onto K. As characterized in Proposition 1, an eigenfunction is thus obtained as soon as p ∈ ∂J(u) and p = Proj K (u).

t = 0 t = 1 t = 2

Uniqueness of a solution of (27)

We start our analysis of the flow by stating a comparison result.

Proposition 6. Let u and v be two solutions of [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF] with respective initial condition u 0 and v 0 such that J(u 0 ) < +∞ and J(v 0 ) < +∞, with u 0 2 = v 0 2 = 1. Then we have:

d dt 1 2 u -v 2 2 ≤ J(u) + J(v) 2 u -v 2 2 ( 28
)
Uniqueness of a solution for the flow is then a direct consequence, as stated in the next corollary.

Corollary 2. Let u and v be two solutions of [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF] with respective initial condition u 0 and v 0 , such that J(u 0 ) < +∞ and J(v 0 ) < +∞, with u 0 2 = v 0 2 = 1. Then we have:

d dt 1 2 u -v 2 2 ≤ J(u 0 ) + J(v 0 ) 2 u -v 2 2 ( 29
)
Moreover, we have:

u -v 2 2 ≤ u 0 -v 0 2 2 exp ((J(u 0 ) + J(v 0 )(t -t 0 )) ( 30
)
The corollary is a direct consequence of the previous proposition, the fact that J(u) is decreasing, and Gronwall lemma. Let us now prove the proposition.

Proof. From the properties of the flow, we have u 2 = v 2 = 1. Moreover, we have J(u) ≤ J(u 0 ) and J(v) ≤ J(v 0 ) for all t. Let us compute:

d dt 1 2 u -v 2 2 = u -v, u t -v t = J(u)u -p -J(v)v + q, u -v , p ∈ ∂J(u) q ∈ ∂J(v) = J(u)u -J(v)v, u -v + q -p, u -v , p ∈ ∂J(u) q ∈ ∂J(v)
But q -p, u -v ≤ 0 since p ∈ ∂J(u) and q ∈ ∂J(v) (as J is convex). Hence:

d dt 1 2 u -v 2 ≤ J(u)u -J(v)v, u -v (31)
But:

J(u)u -J(v)v, u -v = -J(u) u, v + J(u) + J(v) -J(v) u, v = (J(u) + J(v))(1 -u, v ) = J(u) + J(v) 2 u -v 2 2
So that we get [START_REF] Hein | An inverse power method for nonlinear eigenproblems with applications in 1-spectral clustering and sparse pca[END_REF].

Notice of course that uniqueness of the solution of ( 27) comes at once from [START_REF] Merkurjev | Global binary optimization on graphs for classification of high-dimensional data[END_REF].

Hence:

χ ||.|| 2 ≤1 (v) ≥ χ ||.|| 2 ≤1 (u k+1 ) + v -u k+1 , - u k+1 -u k δt + J(u k )u k+1 -p k+1
(48) We can let k → +∞ so that:

χ ||.|| 2 ≤1 (v) ≥ χ ||.|| 2 ≤1 (u) + v -u, J(u)u -p . ( 49
)
We thus deduce that (44) holds.

Infinite dimensional case

In this case, we consider that J is defined on X := L 2 (Ω). We first state a preliminary convergence result.

Proposition 8. Let u 0 in L 2 (Ω), and the sequence u k defined by [START_REF] Moeller | Learning nonlinear spectral filters for color image reconstruction[END_REF]. There exists some u and p in X such that up to a subsequence, u k converges to u and p k converges to p in L 2 (Ω) weak.

Proof. Since u k 2 = 1 for all k, u k is a bounded sequence in L 2 (Ω). Hence (see e.g. [START_REF] Brezis | Analyse fonctionnelle. Théorie et applications[END_REF]) there exists u in L 2 (Ω) such that up to a subsequence, u k u in L 2 (Ω) weak. With the same reasoning, we can show that there exists p in L 2 (Ω) such that up to a subsequence, p k p in L 2 (Ω) weak.

To state a full convergence result, we need to add a technical hypothesis.

Theorem 3. Let u 0 in L 2 (Ω), and the sequence u k defined by [START_REF] Moeller | Learning nonlinear spectral filters for color image reconstruction[END_REF]. Assume that the sequence u k lives in a compact subset of L 2 (Ω) for the strong topology. There exists some u and p in L 2 (Ω) such that up to a subsequence, u k converges to u in L 2 (Ω) strong and p k converges to p in L 2 (Ω) weak, with p ∈ ∂J(u), and J(u k ) converges to J(u). Moreover, u satisfies the differential inclusion:

J(u)u -u -p ∈ ∂χ ||.|| 2 ≤1 (u) (50)
Proof. From the previous proposition, there exists u in L 2 (Ω) such that up to a subsequence, u k u in L 2 (Ω) weak. There exists also p in L 2 (Ω) such that up to a subsequence, p k p in L 2 (Ω) weak. Since u k lives in a compact subset of L 2 (Ω) for the strong topology, we have that

u k → u in L 2 (Ω) strong.
The rest of the proof is identical to the one of Theorem 2, since

X := L 2 (Ω).
where the gradient magnitude is based on 2 , |∇u| = (u x ) 2 + (u y ) 2 . On the top right and middle the results of u and p, respectively, are shown after 100 iterations. The initialization can be noise or some image to produce different eigenfunctions, here we chose the cameraman image. It can be seen that p is very similar to u in its shape, which is expected for eigenfunctions. In the ideal case, we should expect p = λu pointwise, therefore for every pixel the ratio p/u = λ should yield spatially a constant image. In Fig. 3 top right this ratio is shown, where most of the image is of constant value, but there are some deviations near the boundaries of the shape. There is still no definitive theory of eigenfunctions of discrete isotropic TV. Our experiments indicate that numerically one reaches in general only approximations of eigenfunctions of the continuous case. Convergence to precise eigenfunctions are reached in trivial cases, such as partitions by straight lines of the space. Consequently, the process is very stable when p and u are very similar, but full convergence is not attained numerically, as can be seen in Fig. 3 For segmentation and clustering purposes, TV on graphs is used as the regularizer of choice. One constructs a graph based on the input data and computes the iterative flow. We use J as defined in (3) with q = 1 where w ij is the weight of the graph between node i and node j. We give examples of graph based on an image, for segmentation, and one based on point cloud, for clustering. In Fig. 5 the graph is constructed from the image based on Euclidean patch-distances, as for instance in [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF]. We use a 5 × 5 search window for similar patches, so the pixel proximity relation is essentially very local. Our initialization of the flow is the input image f (left). We show the result of u after 50 iterations. We see that the process naturally converges to a segmentation of the data (see thresholded result on the right).

Conclusion

In this work we have presented a class of nonlinear flows for which their steady-states are eigenfunctions with respect to the subgradient of a desired one-homogeneous regularization functional, such as any flavor of totalvariation on grids or graphs. The flows were analyzed in finite dimensions both in the continuous time setting and in the discrete setting. The discrete setting is realized as a series of convex optimization iterations. Its properties and stability characteristics were shown. For a specific case of the proposed α-flow, a comprehensive theory was derived. It was shown that the discrete iterations converge to a steady state, which is an eigenfunction. Moreover, we have shown that the time continuous flow exists and has a unique solution.

These algorithms can be used for several applications related to segmentation and clustering, where graph total variation and eigenfunctions of the 1-Laplacian operator are used (see e.g. [START_REF] Bresson | Total variation, cheeger cuts[END_REF][START_REF] Bresson | Multi-class transductive learning based on 1 relaxations of cheeger cut and mumford-shahpotts model[END_REF][START_REF] Bühler | Spectral clustering based on the graph plaplacian[END_REF]). The flows are continuously evolving towards an eigenfunction and are very convenient to use when one has a rough initial estimate (for instance, using linear eigenfunctions as approximations). Future research will include further examination of the general continuous case. In addition, we will examine the simultaneous evolution of several flows to compute several eigenfunctions, e.g. for multi-class clustering applications.

A.2 A priori estimates

We first need to show some a priori estimates.

Proposition 10. Let u 0 in L 2 (Ω). Then t → J(ǔ δt (t, .)), and t → pδt (t, .) 2 are non increasing, ǔδt (t, .) 2 = u 0 2 for all t, and ûδt (t, .) 2 ≤ 3 u 0 2 .

Proof. This is a direct consequence of the previous section and equation (55).

Proposition 11. Let T > 0 be fixed. There exists a constant C > 0, which does not depend on δt, such that:

T 0 ∂ ûδt ∂t 2 L 2 (Ω) ≤ C (59)
Proof. We have:

t k+1 t k ∂ ûδt ∂t 2 L 2 (Ω) = δt Ω u k+1 (x) -u k (x) δt 2 dx. (60)
By using (42), we get:

t k+1 t k ∂ ûδt ∂t 2 L 2 (Ω) ≤ 2 (J(u k ) -J(u k+1 )) . (61)
Let us denote by K = [T /δt], then

K-1 n=0 t k+1 t k ∂ ûδt ∂t 2 L 2 (Ω) ≤ 2 (J(u 0 ) -J(u K )) ≤ 2J(u 0 ).
We thus deduce that:

T 0 ∂ ûδt ∂t 2 L 2 (Ω) dt ≤ 2T J(u 0 ) + T t K ∂ ûδt ∂t 2 L 2 (Ω) dt. (62)
But, by using (42), we have:

T t K ∂ ûδt ∂t 2 L 2 (Ω) dt ≤ 2 T -t K δt (J(u K ) -J(u K+1 )) ≤ 2J(u 0 ).
We then get from Proposition 10 that there exists B > 0 which does not depend on K and δt such that:

T t K ∂ ûδt ∂t 2 L 2 (Ω)
dt ≤ B. We then conclude thanks to (62).

Corollary 3. Let T > 0 be fixed. Then:

lim δt→0 T 0 ûδt -ǔδt 2 L 2 (Ω) dt = 0 (63)
Proof. Let us denote by K = [t/δt]. We have:

T 0 ûδt -ǔδt 2 L 2 (Ω) dt = K-1 k=0 t k+1 t k ûδt -ǔδt 2 L 2 (Ω) dt+ T t K ûδt -ǔδt 2 L 2 (Ω) dt, (64) but:
K-1 k=0 t k+1 t k ûδt -ǔδt 2 L 2 (Ω) dt = K-1 k=0 t k+1 t k (t -t k -δt)(u k+1 -u k ) 2 L 2 (Ω) dt (65)
We then deduce from (56) that:

K-1 k=0 t k+1 t k ûδt -ǔδt 2 L 2 (Ω) dt ≤ K-1 k=0 t k+1 t k δt ∂ ûδt ∂t 2 L 2 (Ω) dt ≤ (δt) 2 T 0 ∂ ûδt ∂t 2 L 2 (Ω)
dt →0 as δt → 0 , and:

T t K ûδt -ǔδt 2 L 2 (Ω) dt ≤ (δt) 3 u K+1 -u K δt 2 L 2 (Ω) dt →0 as δt → 0 . (66)

A.3 Convergence

Theorem 4. Uniqueness of Theorem4 comes from Corollary 2 (see also e.g. [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]). Let us now prove existence of a solution. We first remark that, from Propositon 10 and 11, ûδt is uniformly bounded in W 1,2 ((0, T ); L 2 (Ω)), since

u 2 W 1,2 ((0,T );L 2 (Ω)) = T 0 ∂ ûδt ∂t 2 L 2 (Ω) + T 0 ûδt 2 L 2 (Ω) .
Thus, up to a subsequence, there exists u in W 1,2 ((0, T ); L 2 (Ω)) such that ûδt u in W 1,2 ((0, T ); L 2 (Ω)) weak. Since W 1,2 ((0, T ); L 2 (Ω)) is compactly embedded in L 2 ((0, T ); L 2 (Ω)) (see [START_REF] Temam | Navier Stokes equations[END_REF], Theorem 2.1, chapter 3), ûδt → u strongly in L 2 ((0, T ); L 2 (Ω)).

Since ǔδt 2 = 1 for all t ∈ (0, T ), we have ǔδt Proposition 10). Hence ǔδt is a bounded sequence in L 2 ((0, T ); L 2 (Ω)). Thus there exists ǔ in L 2 ((0, T ); L 2 (Ω)) such that up to a subsequence, ǔδt ǔ in L 2 ((0, T ); L 2 (Ω)). Since we assume that ǔδt lives in some compact set for the strong topology of L 2 ((0, T ); L 2 (Ω)), we deduce that ǔδt → ǔ in L 2 ((0, T ); L 2 (Ω)) strong. We can also show that there exists p in L 2 ((0, T ); L 2 (Ω)) such that up to a subsequence, pδt p in L 2 ((0, T ); L 2 (Ω)) weak. From Proposition 3, we deduce that ǔδt → u strongly in L 2 ((0, T ); L 2 (Ω)), so that ǔ = u.

2 L 2 ((0,T );L 2 (Ω)) = T 0 ǔδt 2 L 2 (Ω) = T (thanks to
Let v in L 2 (Ω). We have:

J(v) ≥ J(ǔ δt ) + v -ǔδt , pδt . (67)
Let φ in C 0 c (0, T ) a test function, φ ≥ 0. We multiply (67) by φ and integrate on (0, T ):

This inequality holds for all φ ≥ 0, we deduce that for a.e. t in (0, T ):

J(v) ≥ J(u) + Ω (v -u)p dx. (72)
Hence p ∈ ∂J(u). Moreover, we have J(ǔ δt ) = ǔδt , pδt . Letting again δt → 0, we see that J(u δt ) → J(u) = u, p . The semi-discrete implicit scheme writes for pδt in ∂J(u ǔδt ) and for a.e. t ∈ (0, T ): (77)

Now, since ǔδt → u strongly in L 2 ((0, T ); L 2 (Ω)), ∂ ûδt ∂t ∂u ∂t in L 2 ((0, T ); L 2 (Ω)) weak, J(ǔ δt )ǔ δt → J(u)u strongly in L 2 ((0, T ); L 2 (Ω)), and pδt → p strongly in L 2 ((0, T ); L 2 (Ω)), the second term on the right hand-side of (76) tends to (78) This inequality holds for all φ ≥ 0, we deduce that for a.e. t in (0, T ):

χ ||.|| 2 ≤1 (v) ≥ χ ||.|| 2 ≤1 (u) + Ω (v -u) - ∂u ∂t + J(u)u -p dx, (79)
i.e.: -∂u ∂t + J(u) -p ∈ ∂χ ||.|| 2 ≤1 (u). Hence we deduce that u is a solution of (51) in the distributional sense.

Figure 1 :

 1 Figure 1: Illustration of an eigenfunction u where λu ∈ ∂J(u) ⊂ K. Observe that λu is the orthogonal projection of u onto K.

Figure 2 :

 2 Figure 2: Illustration of the evolution of u(t) with an explicit discretization of the flow (27). The vector p ∈ ∂J(u) and its projection onto u gives the direction u t of the flow.

) u k+1 -u k 2 Figure 3 :

 23 Figure 3: Results of the flow for isotropic TV. Top: u and p and the ratio u/p after 100 iterations (the ratio should be close to a constant function). Bottom, the values of J(u k ) and u k+1 -u k are plotted as a function of iterations k.

Fig. 4 2 Figure 4 :

 424 Fig. 4 shows the case of local anisotropic TV, where the gradient magnitude is based on 1 , |∇u| = |u x | + |u y |. In this case, u and p have exactly the same shape, and the ratio p/u (top-right) is constant, up to numerical precision. As we reach a precise eigenfunction the algorithm fully converges to a steady state, as seen (bottom) on the values of J(u k ) and u k+1 -u k 2 as a function of the iteration k. These experiments are useful to examine the algorithm and to compute local discrete TV eigenfunctions.

Figure 5 :

 5 Figure 5: Results of the flow for TV defined on graphs constructed from the image, based on patch distances. Initializing with the input image, the flow yields a segmentation of the image.

Figure 6 :

 6 Figure 6: Results of the flow for TV defined on graphs based on point cloud distances. The processes converges to natural clustering of the data.

  δt )φ dt + T 0 v -ǔδt , pδt φ(t) dt, δt )φ(t) dt + T 0 Ω (v -ǔδt )p δt φ(t) dtdx. (69)We want to let δt → 0 in (69). By convexity, we have:lim inf T 0 J(ǔ δt )φ(t) dt ≥ T 0 J(u)φ(t) dt.(70)Now, since ǔδt → u strongly in L 2 ((0, T ); L 2 (Ω)) and pδt → p in L 2 ((0, T ); L 2 (Ω)) strong, the second term on the right hand-side of (69) tends to u)pφ(t) dtdx.

0 χ 0 χχχ 0 χ

 000 x) + J(ǔ δt (t, x))ǔ δt (t, x)) -pδt (t, x) ∈ ∂χ ||.|| 2 ≤1 (ǔ δt (t, x))), (73)We thus have for all v in L 2 (Ω), and a.e. t ∈ (0, T ):χ ||.|| 2 ≤1 (v) ≥ χ ||.|| 2 ≤1 (ǔ δt ) + v -ǔδt , -∂ ûδt ∂t + J(ǔ δt )ǔ δt -pδt . (74)Let φ in C 0 c (0, T ) a test function, φ ≥ 0. We multiply (74) by φ and integrate on (0, T ):T ||.|| 2 ≤1 (v)φ(t) dt ≥ T ||.|| 2 ≤1 (ǔ δt )φ(t) dt+ T 0 v -ǔδt , -∂ ûδt ∂t + J(ǔ δt )ǔ δt -pδt φ(t) dt, ||.|| 2 ≤1 (v)φ(t) dt ≥ T 0 χ ||.|| 2 ≤1 (ǔ δt )φ(t) dt+ T 0 Ω (v-ǔ δt ) -∂ ûδt ∂t + J(ǔ δt )ǔ δt -pδt φ(t)dtdx. ||.|| 2 ≤1 (ǔ δt )φ(t) dt ≥ T ||.|| 2 ≤1 (u)φ(t) dt.

  u)u -p φ(t) dtdx.

Acknowledgements The authors would like to thank Antonin Chambolle for many discussions related to this work. This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR GOTMI) (ANR-16-CE33-0010-01). Jean-Franois Aujol is a member of Institut Universitaire de France (IUF). Guy Gilboa acknowledges support by the Israel Science Foundation (grant No. 718/15). This study as been carried out with financial support from EU RISE project NoMADS.

Properties of a semi-explicit scheme

For α = 1, the numerical scheme [START_REF] Bresson | Total variation, cheeger cuts[END_REF] becomes:

-p k+1 [START_REF] Merriman | Diffusion generated motion by mean curvature[END_REF] associated with the minimization of

u k+1 is a minimizer of F (., u k ), as soon as 1 δt > J(u 0 )

. From Proposition 4, the continuous flow keeps ||u|| 2 constant for α = 1, but the discrete properties studied in Proposition 5 just ensure that ||u k || 2 is non decreasing. As a consequence, instead of dealing with [START_REF] Meyer | Oscillating patterns in image processing and in some nonlinear evolution equations[END_REF], we consider the following renormalization to ensure that for

This scheme is associated with the minimization of

Proposition 7. u k+1 defined in [START_REF] Moeller | Learning nonlinear spectral filters for color image reconstruction[END_REF] is the minimizer of F .

Proof. We define u k+1 as the minimizer of [START_REF] Müller | Advanced image reconstruction and denoising: Bregmanized (higher order) total variation and application in pet[END_REF]. Then:

So that:

Hence

We deduce that u k+1 is the L 2 projection on the ball of radius 1 of

But from [START_REF] Merriman | Diffusion generated motion by mean curvature[END_REF], since u k 2 = 1, we know that:

So we see that u k+1 is the L 2 projection on the ball of radius 1 of u k+1/2 . Moreover, since the scheme defined by [START_REF] Merriman | Diffusion generated motion by mean curvature[END_REF] is such that u k+1/2 2 ≥ u k 2 , we deduce that u k+1/2 2 ≥ 1, and thus

Hence u k+1 is also solution of [START_REF] Moeller | Learning nonlinear spectral filters for color image reconstruction[END_REF].

Thanks to Proposition [START_REF] Bellettini | The total variation flow in R N[END_REF], we are now in position to analyse the sequence u k defined by [START_REF] Moeller | Learning nonlinear spectral filters for color image reconstruction[END_REF].

Theorem 1. Let u 0 in X, and the sequence u k defined by [START_REF] Moeller | Learning nonlinear spectral filters for color image reconstruction[END_REF]. Then the sequences J(u k ) and p k 2 are non increasing, u k 2 = u 0 2 for all k, and u k+1 -u k → 0.

Proof. ||u k || 2 constant. Let us proceed by contradiction and assume that ||u k || 2 = 1 and ||u k+1 || 2 < 1, where u k+1 is the minimizer of [START_REF] Müller | Advanced image reconstruction and denoising: Bregmanized (higher order) total variation and application in pet[END_REF]. As the constraint χ ||.|| 2 ≤1 (u k+1 ) is not saturated then

, we can use Proposition 5 with α = 1 and therefore conclude.

Convergence of u k+1 -u k . Assume that J(u 0 ) < 1 δt to have F convex, and ||u k || 2 = 1, then:

since

Summing on k from 0 to N -1 relation [START_REF] Zeune | Multiscale segmentation via bregman distances and nonlinear spectral analysis[END_REF], we deduce that:

In particular k u k+1 -u k 2 2 converges, and u k+1 -u k → 0.

Convergence of the semi-implicit scheme

We are now in position to state some convergence results. We first consider the case when X is a finite dimensional space. We will then consider the general case when X is an infinite dimensional case, for which we need to add a technical hypothesis to get a convergence result.

Finite dimensional case

Theorem 2. Let u 0 in X, with X of finite dimension, and the sequence u k defined by [START_REF] Moeller | Learning nonlinear spectral filters for color image reconstruction[END_REF]. There exists some u and p in X such that up to a subsequence, u k converges to u in X and p k converges to p in X, with p ∈ ∂J(u), and J(u k ) converges to J(u). Moreover, u satisfies the differential inclusion:

Proof. From Theorem 1, there exists u in X such that up to a subsequence, u k → u in X. There exists also p in X such that up to a subsequence, p k → p in X.

Let v in X. We have:

We can let k → +∞ and using the lower semi continuity of J we get:

Hence p ∈ ∂J(u). Moreover, we have J(u k ) = u k , p k . Letting again k → +∞, we se that J(u k ) → J(u) = u, p .

Let again v in X. We have, for p k+1 in ∂J(u k+1 ):

Eigenfunction

The following result shows that the limit u of the semi-implicit scheme is indeed an eigenfunction.

Proposition 9. If u satisfies either Equation 44 or Equation 50, then u is an eigenfunction.

Proof. From either Equation 44 or Equation 50 it follows that there exists p ∈ ∂J(u) and q ∈ ∂χ ||.|| 2 ≤1 (u) such that: -J(u)u+p+q = 0. Since ||u|| 2 = 1, we have ∂χ ||.|| 2 ≤1 (u) = {γu, γ ≥ 0}, and therefore:

and u is an eigenfunction. Moreover, as J(u) = p, u , we get γ = 0. Now that we have analyzed a semi-explicit scheme for computing an eigenfunction, we turn our attention to the time continuous problem (evolution equation) in the next section.

Existence of solution for (27)

In this section the time continuous flow is analyzed. Let us rewrite here [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF]:

We have the following existence and uniqueness result. As in the case of Theorem 3, we need to add a technical hypothesis. Theorem 4. Let u 0 in L 2 (Ω) with J(u 0 ) < +∞. Assume that ǔδt lives in some compact set K T for the strong topology of L 2 ((0, T ); L 2 (Ω)). Then problem(51) admits exactly one solution in W 1,2 ((0, T ); L 2 (Ω)).

The proof of this theorem is detailed in Appendix A

Numerical Results

Here we give a few examples of running the flow in several settings. We first examine local TV regularizers. Fig. 3 shows the results using isotropic TV, A Existence of a solution for [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF] To show that problem [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF] has a solution, we start from the semi-discrete problem [START_REF] Moeller | Learning nonlinear spectral filters for color image reconstruction[END_REF]. Thanks to Proposition [START_REF] Bellettini | The total variation flow in R N[END_REF], we therefore consider a sequence (u k ) satisfying [START_REF] Müller | Advanced image reconstruction and denoising: Bregmanized (higher order) total variation and application in pet[END_REF]. We know that (u k ) satisfies:

with p k+1 in ∂J(u k+1 ). From the result of Section (4),we know that the sequence (u k ) exists and is unique provided δt small enough.

A.1 Definitions of interpolate functions

For t 0 = 0 and t k = kδt, we classically introduce two piecewise constant functions defined on Ω × R + (see e.g. [START_REF] Belahmidi | Time-delay regularization of anisotropic diffusion and image processing[END_REF][START_REF] Aubert | A variational approach to removing multiplicative noise[END_REF]):


Shalika S W Walker
email: s.w.walker@tue.nl

Warody Lombardi

Suzanne Lesecq

Samira Roshany-Yamchi

Application of Distributed Model Predictive Approaches to Temperature and CO 2 Concentration Control in Buildings

Keywords: Model Predictive Control, centralized, distributed, RC model, Energy management, temperature, CO 2 concentration

In the context of energy consumption reduction, this paper focuses on the application of Model Predictive Control to occupants' thermal comfort together with indoor air quality control while improving the whole building energy efficiency. First, an open-space office split in three zones, located in Cork Institute of Technology, is modeled. A centralized MPC is designed to control the temperature and CO 2 concentration in the three zones. Then, a distributed version of the MPC, with three separate local controllers, is considered. Finally, simulation results show that the distributed MPC solution achieves control performance quite close to the centralized version with less computing effort.

INTRODUCTION

Buildings consume more than 40% of the total primary energy resources throughout the world [START_REF] Shaikh | A review on optimized control systems for building energy and comfort management of smart sustainable buildings[END_REF], [START_REF] Cigler | Model Predictive Control for Buildings[END_REF]]. Moreover, inefficiencies of the deployed sensing and control strategies cause energy waste that should be avoid by a better coordination among Building Automation Systems (BMS) and appropriate control approaches. The minimization of the energy consumed by buildings is essential for their sustainability. However, this minimization may badly affect the occupants' comfort, e.g. by reducing (resp. increasing) the temperature in the building when the outside temperature is low (resp. high). Basically, a higher degree of indoor comfort is expected by the occupants along with the increased time they spend inside buildings. While performing their daily activities in buildings, energy savings should not negatively impact occupants' health or decrease their welfare [START_REF] Castilla | A multivariable nonlinear mpc control strategy for thermal comfort and indoor-air quality[END_REF]], thus leading to contradictory objectives [START_REF] Wang | Multi-objective optimal analysis of comfort and energy management for intelligent buildings[END_REF]] at control level. Energy and Comfort Management Systems try to fulfill occupants' comfort expectations while reducing energy consumption.

Indoor Environmental Quality (IEQ) is related to thermal aspects [START_REF] Sarbu | Aspects of indoor environmental quality assessment in buildings[END_REF]], Indoor Air Quality (IAQ), acoustic and visual (lighting) levels [START_REF] Castilla | A multivariable nonlinear mpc control strategy for thermal comfort and indoor-air quality[END_REF]], while humidity level also affects the This work is supported by TOPAS H2020 project, GA nb 676760 (https://www.topas-eeb.eu/). This work has been performed while S. Walker and W. lombardy were resp. internship and research-engineer in CEA-LETI. comfort feeling. IEQ/IAQ is an active research area, from clinical and medical viewpoints, to control, but also building construction and retrofit, communication, etc. These research works are often related to energy efficiency and HVAC control. IAQ regulations are regularly reinforced (e.g. French statutory-orders n2011-1728 (air quality in public buildings) and n2011-1727 (benzene, formaldehyde and CO 2 levels)) while building energy efficiency is encouraged via directives, e.g. [Union (2010)]. This situation advocates for more multi-disciplinary researches to better optimize buildings, taking into account occupants' comfort and advanced control of BMS [START_REF] Shaikh | A review on optimized control systems for building energy and comfort management of smart sustainable buildings[END_REF]].

Air quality, either indoor or outdoor, is one of the major health concerns [START_REF] Gurjar | Air pollution: health and environmental impacts[END_REF]]. Moreover, quality of the indoor environment has strong health impact because of the close structure of the buildings. Indoor air pollutants are mostly emitted from sources inside the buildings but they can also enter from outside. CO 2 is usually considered the main IAQ indicator. Indoor CO 2 is mainly produced by the occupants exhalation. Poor indoor air quality can lead for the occupants to suffer from Sick Building Syndrome and building related illnesses. Thermal comfort is defined as "the condition of mind which expresses satisfaction with the thermal environment" in the international regulations ISO-7730 and ASHRAE-55 [ASHRAE (1992)]. Therefore, it is related to conscious intellectual activity influenced by the physical, psychological and physiological factors of the occupants.

Buildings have been equipped with BMS to manage the various systems installed, from lighting, ventilation, heating, but also fire alarm, and security. When properly tuned, they can offer better energy efficiency for the building, for instance by switching off the light when nobody is in a given area. With the integration of actuator/sensor networks, these systems can monitor and control the interior conditions so as to fulfill IEQ while keeping the energy consumption at a minimum. Unfortunately, this ideal situation is seldom reality because of the contradictory influence of building energy efficiency and occupants' comfort on the energy bill. Thus, building controllers should take into account multiple objectives, or at least, take a Multi-Input Multi-Output (MIMO) viewpoint to deal with such constrained systems.

The present paper presents an application of Model Predictive Control (MPC) to deal with thermal and CO 2 concentration control in indoor environment. Temperature T and CO 2 concentration are contradictory by nature. Basically, the CO 2 concentration can be easily decreased by the injection of fresh outdoor air inside the building using for instance forced or natural ventilation. The side effect of this injection is the decrease of the indoor temperature when the outdoor one is colder, leading to more heating. Both aspects are physically coupled and can be modeled using a MIMO model. The proposed controller is applied to an open-space office, split in three zones, and located in Cork Institute of Technology (CIT), Ireland. A centralized MPC approach is first designed. It requires all data to be collected in a unique control point. For large buildings, and in presence of communication issues, this can lead to bad functioning of the control law. Thus, a distributed version of the controller, with loosely coupled areas, might be of interest to overcome these situations. Such a distributed approach is therefore proposed. This Distributed MPC implements a separate local controller for each zone.

The rest of the paper is organized as follows. Section 2 recalls basics and notations for MPC. Section 3 presents the control-oriented thermal model developed here for the setup considered. Sections 4 and 5 are dedicated to the design of the centralized and distributed MPC approaches. Section 6 shows some simulation results and discusses implementation cost. A conclusion gives future work directions.

BASICS ON MODEL PREDICTIVE CONTROL

As a control methodology, Model Predictive Control (MPC) can naturally deal with MIMO systems subject to physical constraints. The control law is computed in real-time and on-line via the solution of an optimization problem. Contradictory objectives, e.g. building energy efficiency and occupants' comfort, can be simultaneously considered.

Summary of MPC for Building Management

MPC has gained popularity among researchers and industry in various engineering fields as an effective approach to deal with multi-variable constrained control problems. In the building sector, MPC is used for building climate control, cooling, heating, ventilation, etc. The main aim of MPC when applied to building management is to reduce the energy consumption, taking into account occupants' comfort expectation [START_REF] Cigler | Model Predictive Control for Buildings[END_REF]]. MPC design is traditionally a three-step process [START_REF] Camacho | Model predictive control[END_REF]]:

(1) Modeling of the system. The efficiency and accuracy of the control actions highly depend on the system model. This latter is expressed with a statespace representation, which for a Linear Time Invariant (LTI) model is given by:

x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) (1)
where y(k) ∈ IR m is the measurement (output), x(k) ∈ IR n is the system state and the control (input) vector is u. Matrices A, B, C are real of appropriate dimension.

(2) Cost objective function. Here, the cost function is chosen to minimize the energy consumption and ensure IEQ. Energy consumption and IEQ must be satisfied under constraints on the system outputs and inputs. In the present context, the square of the control vector u is related to energy flow. Therefore, minimizing u 2 is equivalent to the energy consumption minimization. This can be achieved through the minimization of the cost function:

J = N j=1 { (y -y ref )(k + j) 2 Q Comfort + ∆u(k + j) 2 R Energy } (2)
under the constraints Eq.( 1)

y min ≤ y ≤ y max , u min ≤ u ≤ u max ∆u min ≤ ∆u ≤ ∆u max
where ∆u(k) = u(k) -u(k -1). y ref is the comfort level (temperature and CO 2 concentration) that should be maintained. Q > 0 and R > 0 are weighting matrices. The first part of the cost function is related to the comfort objective while the second part deals with the energy consumption. Note that the ∆ formulation implies the usage of an integrator. (3) Prediction horizon N and weighting matrices Q, R. These tuning parameters influence the closed loop response time, and the relative weight of both objectives, namely comfort and energy efficiency.

MPC architectures

Three main Model Predictive Control architectures can be found in the literature, namely Centralized MPC, Decentralized MPC and Distributed MPC [START_REF] Christofides | Distributed model predictive control: A tutorial review and future research directions[END_REF], [START_REF] Scattolini | Architectures for distributed and hierarchical model predictive control-a review[END_REF]]. Here, the centralized version is used as reference while the distributed one is considered so as to reduce the computational and communication burdens of the centralized version. Basically, the distributed MPC implements local controllers with smaller size optimization problems to be solved. The centralized and distributed versions are sum up as follows:

• Centralized MPC (CMPC). MPC is classically implemented in a centralized scheme where all the control signals are computed using a single objective function. When the system becomes more complex

BUILDING MODELING FOR MPC

MPC requires a model of the controlled system to predict the system behavior, as briefly recalled in section 2. In this paper, the building thermal behavior is described using an equivalent Resistance-Capacitance (RC) network. The RC parameters are calculated using basic knowledge of the building geometry and construction materials. Mass balance for the CO 2 concentration is used to model air quality. Both models are merged in a state-space model that will be used in the MPC formulation.

Description of the Experimental Building

The practical study is conducted on a two-storey building from CIT, Ireland. Figure 1 shows the area under study, located at the first floor. This area is divided into three zones that are not separated by walls. Moreover, the BMS installed controls the space under study with three separate units. Hereafter, only the model for zone z1 is presented, the models for the other zones being deduced straightforwardly.

Thermal Model

From the heat transfer and heat storage equations, and applying the equivalent thermal Resistance-Capacitance (RC) model, zone z1 is modeled for temperature control. Beside these equations, there is also an effect of the ventilation and of the occupants on the temperature rate • the air in each zone has a unique temperature across its whole volume (lumped model); • the specific heat of the air c p is constant and equal to 1.005KJ/kg.K; • the radiative coupling between the inner walls of the building is ignored since the temperature differences between the walls are negligible; • the metabolic heat gain per occupant is 100W .

The rate of change of temperature T z1 in z1 is given by:

C z1 dT z1 dt = T ww,z1 -T z1 R ww,z1 /2 A ww,z1 + T sw,z1 -T z1 R sw,z1 /2 A sw,z1 + T f,z1 -T z1 R f,z1 /2 A f,z1 + T c,z1 -T z1 R c,z1 /2 A c,z1 + ϕ + T r,z1 -T z1 R r,z1 A r,z1 + εσA r,z1 (T 4 r,z1 -T 4 z1 ) + Q z1 • cp • (Tout -T z1 ) (3)
where ϕ is the heat flux. C z1 is the thermal capacitance. ww is the west wall, sw the south wall, f the floor, and c the ceiling. R a,z1 , T a,z1 and A a,z1 are respectively the thermal resistance, temperature and area related to the zone limit a ∈ {r, c, f, sw, ww}. Q z1 is the ventilation rate in z1. In the sequel, denote f 1 = dT z1 /dt. Similarly to (3), f a = dT a,z1 /dt, a ∈ {r, c, f, sw, ww} are defined. Due to lack of space, they are not reported here. Define:

f z1 = [f 1 , f r , f c , f f , f sw , f ww ] T (
4) Moreover, the dynamic thermal equations reflect the absence of wall between zones by adding a natural convective heat transfer term (not reported here).

Air Quality Model

The mass balance for the CO 2 concentration in z1 is:

g 1 = dC0 2,z1(t) dt = (G(t) * n o,z1 + Q z1 {C0 2,out (t) -C0 2,z1 (t)})/V z1 (5)
where G(t) represents the CO 2 generated per person [START_REF] Emmerich | State-of-the-art review of CO2 demand controlled ventilation technology and application[END_REF]] in z1, n o,z1 is the number of occupants in z1 and V z1 is the total volume of z1. Note that the diffusion between zones that are not separated by walls is added using the Fick's law to obtain a more realistic mass balance equation (not reported here).

Model Linearization

The state-space representation is obtained via the linearization of the non linear dynamic equations using the Jacobian linearization around an equilibrium point.

Thermal Model Linearization Define the state vector for the thermal model in z1 as:

x T,z1 = [T z1 T r,z1 T c,z1 T f,z1 T sw,z1 T ww,z1 ] T (6) From ( 4), the state matrix for z1 is given by:

A T,z1 = ∂f z1 ∂x T,z1 (7)
The experimental building uses natural ventilation only.

The control variable is u T,z1 = [N z1 P z1 ] T (8) where N z1 corresponds to the windows opened. For simplicity, N z1 = no w=1 N z1,j where n o is the number of windows opened. N z1,j ∈]0, 1] is the window j position where 0 (resp. 1) means "window j is closed (resp. fully open)". P z1 ∈ [0, P z1,max ] is the power consumption of the heaters. The natural ventilation rate for a single sided ventilation system can be computed with:

Q = 0.025 • W S • A w • N w1 (9)
where A w is the effective area of a window, N w1 is the number of windows in z1 and W S is the wind speed. Then:

B T,Z1 = ∂f z1 ∂u T,z1 (10)
Air Quality Model Linearization Define the state and control vectors for the air quality model: 5), it comes:

x CO2,z1 = [CO 2,z1 ] u CO2,z1 = [N z1 ] (11) From (
A CO2,z1 = ∂g 1 ∂x CO2,z1 B CO2,z1 = ∂g 1 ∂u CO2,z1 (12)
3.5 Linearized model for z1

From ( 6), ( 8) and ( 11), it comes:

x z1 = [x T T,z1 x CO2,z1 ] T u z1 = [N z1 P z1 ] T (13)
The state and control matrices are derived from ( 7), ( 10) and ( 12):

A z1 = A T,z1 0 0 A CO2,z1 B z1 = A T,z1 B CO2,z1 0 (14)

CENTRALISED MPC

The whole system under study is modeled with a state space representation derived from section 3:

x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) (15) where x = [x T z1 x T z2 x T z3 ] T u = [u T z1 u T z2 u T z3 ] T (16)
and

A = A z1 A 12 A 13 A 21 A z2 A 23 A 31 A 32 A z3 B = B z1 0 0 0 B z2 0 0 0 B z3 (17)
A ij , i = j, accounts for the coupling between zones i and j. Matrix C is defined according to the measured outputs:

y = [T z1 CO 2,z1 T z2 CO 2,z2 T z3 CO 2,z3
] T (18) Since the state vector x cannot be fully measured, a Luenberger observer has been designed in order to estimate the different temperatures in the three zones.

First, a centralized Model Predictive Controller is designed. The three zones are controlled thanks to the minimization of a unique cost function:

min ∆u(k) J = min ∆u(k) N j=1 { y -y ref 2 Q Comfort + ∆u(k) 2 R Energy } (19)
subject to

Eq.( 15)

y min ≤ y ≤ y max , u min ≤ u ≤ u max ∆u min ≤ ∆u ≤ ∆u max (20)
where ∆u(k) = u(k) -u(k -1). y ref contains the comfort references that should be maintained, with thermal and CO 2 parts for each zone:

y T,ref = 20 • C, y CO2,ref = 800ppm (21)
Q > 0 and R > 0 are weighting diagonal matrices. Depending on the numerical value of their elements, one can emphasize more on the energy consumption minimization or on the thermal comfort and/or on the air quality comfort. Their choice is not discussed in the paper.

DISTRIBUTED MPC

An independent, non-iterative partially connected topology is considered for the distributed control scheme. Therefore, only neighboring controllers communicate. For each zone i, i = 1 : 3, a separate controller is designed. The state-space model is defined by:

x i (k + 1) = A ii x i (k) + B ii u i (k) + j =i A ij x j (k) (22) y i = Cx i (k) (23)
where x i = x zi is the state vector defined similarly to (13).

The output vector is given by:

y i = [T zi CO 2,zi ] T (24
) from which matrix C can be derived straightforwardly. The state and control matrices are given by:

A i = A zi B i = B zi (25
) A ij in 22 models the coupling of neighboring zones i and j. As a consequence, all the zones in the considered space do not (necessarily) appear in (22). Each zone is controlled thanks to the minimization of a dedicated cost function:

min ∆ui(k) J i = min ∆ui(k) N j=1 { y i -y i,ref 2 Qi Comfort + ∆u i (k)
y i,ref = [20 • C 800ppm] T (27)
Q i > 0 and R i > 0 are defined per zone, in a similar way as for the centralized version of section 4.

SIMULATION RESULTS

Both MPC strategies are now validated. The first validation step is performed in simulation using Matlab 2015b. Figure 3 illustrates the simulation environment. Here, the control of the ventilation system relies on the monitoring of CO 2 concentration in each zone. This latter is estimated from occupancy count during the time period [6am, 9pm], see Figure 6 where it appears in orange. The numerical values have been extracted from real-life occupancy count in the open-space office considered. The comfort level for CO 2 concentration is fixed at 800 ppm. Thermal comfort is evaluated with the comparison of the indoor temperature in each zone with a reference one equal to 20 • C. Only heating of the zones is considered along with natural ventilation. The outdoor temperature variation is also considered (see Fig. 4). It is extracted from real-life measurements. The wind speed is supposed constant at 5 m/s, which is the average wind speed in Cork area, Ireland, where the experimental building is located.

6.1 Centralized MPC for thermal and CO 2 comfort regulation

Figure 5 shows simulation results for the temperature variation in the three zones with the centralized MPC approach from section 4. The indoor temperature is controlled nearly as expected, with a maximum error of 2 • C. Note that the temperature constraints in (20)have been set to ±1 • C around the reference temperature. When more people enter the controlled space, temperature exceeding appears because of extra heating provided by these people, and that cannot be not anticipated. Moreover, outdoor temperature also provides extra heating. This phenomenon can be decreased at the cost of extra energy consumption for the ventilation that will help extract this extra heat.

Figure 6 shows the CO 2 concentration in the 3 zones together with the occupancy count and C0 2 concentration reference. As can be seen, the controller performs as expected, with the CO 2 concentration following the reference despite disturbances introduced by the occupants. Note that the constraints on CO 2 are set to ±100ppm around the reference CO 2 concentration.

6.2 Distributed MPC for thermal and CO 2 concentration regulation and comparison with the centralized scheme

The same scenario (i.e. external temperature and wind speed, occupancy count, references) is considered for Distributed MPC evaluation. Figures 7 and8 show the evolution of the temperature and CO 2 concentration in the three zones.

Results of both control schemes are reported on these graphs for comparison purpose. As can be seen, they provide close results, with slight under performance for the distributed version. However, this latter presents some nice features. First, it does not require the collection of all the information in a unique controller, decreasing therefore the communication burden, especially when large buildings are considered. Moreover, the distributed scheme is less sensitive to communication breakdown as it processes information locally (the coupling with neighboring zones can be temporarily neglected). Second, the distributed version obviously requires less computational capabilities for each individual controller.

Comparison of energy consumption and computational effort for centralized and distributed schemes

The total energy consumption for the conventional control (i.e. on/off switch of the actuators) and for both MPC strategies is reported in table 6.3. As expected, the energy gain, when compared to the conventional control is slightly smaller for the distributed scheme. Table 6.3 gives an estimate of the computing effort. There results have been obtained using Matlab 2015b. The distributed MPC is taken as reference. Even if these results might be optimized, the computation involved by the centralized controller, when compared to one controller of the distributed scheme, is clearly more demanding, both in terms of memory and computational burden.

CONCLUSION

This paper presents the application of Model Predictive Control to cope with energy efficiency and occupants' thermal and C0 2 concentration comfort in buildings. An open-space office split in three zones, located in Cork Institute of Technology, Ireland, is considered as experimental setup. This space is first modeled. Then, Centralized and Distributed Model Predictive Controllers are designed. Simulation results show that the distributed MPC solution achieves control performance close to the centralized version with less computing effort and communication burden. Both approaches are currently under implementation on the real testbed in the context of the H2020 TOPAs project. Real-life results will be reported during the oral presentation.

Fig. 1 .

 1 Fig.1. Selected space of the building divided into 3 zones and spatially distributed, CMPC requires a significant amount of computing time and data transfer, leading to communication overhead. In case of communication breakdown, the controller is no more able to compute the control signals. As a consequence, for large (scale) systems, as buildings may be, it is highly advisable to divide the system into several (possibly loosely coupled) subsystems.• Distributed MPC (DMPC). In this control architecture, the system is divided into subsystems. Each one has its own local controller that optimizes its own objective function. Information from other subsystems is exchanged between local controllers to achieve better performances. Here, partially connected subsystems are considered. Information about the state in each local system is transmitted to the neighboring local controllers, once in every sampling period[START_REF] Rawlings | Coordinating multiple optimization-based controllers: New opportunities and challenges[END_REF],[START_REF] Farina | Distributed predictive control: a non-cooperative algorithm with neighborto-neighbor communication for linear systems[END_REF],[START_REF] Farina | Distributed noncooperative mpc with neighbor-to-neighbor communication[END_REF]].

Fig. 2 .

 2 Fig. 2. Equivalent RC model for zone z1 of change inside the space. Here, only natural ventilation is used. All these aspects are included in the RC model. The resulting equivalent RC model is given in figure 2 [ Širokỳ et al. (2011),[START_REF] Nagy | Low complexity room thermal models for energy optimization in industrial buildings[END_REF]]. The thermal model z1 is obtained via nodal analysis of the equivalent RC circuit, following the approach in[START_REF] Maasoumy Haghighi | Modeling and Optimal Control Algorithm Design for HVAC Systems in Energy Efficient Buildings[END_REF]]. Consider the assumptions hereafter:

Fig. 3 .

 3 Fig. 3. Simulation setup

Fig. 5 .

 5 Fig. 5. Indoor temperature in the three zones, Centralized MPC

Fig. 7 .

 7 Fig. 7. Indoor temperature in the three zones, centralized (black) and distributed (red) control schemes

  Eq. (22, 23)y i,min ≤ y i ≤ y i,max , u i,min ≤ u ≤ u i,max ∆u i,min ≤ ∆u i ≤ ∆u i,max

	2 Ri	} (26)
	Energy
	subject to

Table 1 .

 1 Energy consump. vs. control strategy

	Energy (kWh)	Conventional	Centralised Distributed
	Total	181	127.83	138.59
	Gain	0 %	30 %	23.5 %

Table 2 .

 2 Computational cost comparison

		Pb. size	Nb. constraints Computational time
				(normalized)
	Centralised	21	660	4.7
	Distributed	7	220	1