Pierre Girardeau
email: pierre.girardeau@cermics.enpc.fr

A COMPARISON OF SAMPLE-BASED STOCHASTIC OPTIMAL CONTROL METHODS

Keywords: 2010 Mathematics Subject Classification. Primary 93E20; Secondary 49M05, 49L20 Stochastic optimal control, Stochastic programming, Particle methods

 and has been widely noticed by practitioners. However, in the second part, we show using the same example that, by mixing Stochastic Programming and Dynamic Programming ideas, the particle method described by Carpentier et al (2009) copes with this numerical difficulty: the number of scenarios needed to obtain a given precision now does not depend on the time horizon. Unfortunately, we also observe that serious obstacles still arise from the system state space dimension.

Introduction

Consider a controlled dynamical system, affected by some exogenous noises, say uncertain parameters for instance. Stochastic optimal control consists in driving this system, having at each time step partial or total observations of those noises, so as to minimize some expected cost integrated through time and while satisfying a number of constraints. Many applications can be handled through such a model. For instance, think of a power producer that has to plan the use of a set of power plants (the control would be the production of the plants at each time step) in order to supply an uncertain power demand (the noises) while minimizing a production cost over a certain time period. We here consider discrete and finite time horizon. We are hence looking for feedback functions which, at each time step and for each possible observation of the system, provide a decision to be taken. We consider the laws of the random variables involved to be continuous; consequently this is an infinite-dimensional optimization problem. In most situations, no analytic solution can be found and one has to consider approximations of the original problem to be able to solve it numerically.

In the manner of the Monte Carlo approach for computing the expectation of a random variable, we seek a tractable approximation of the probabilistic structure of the original problem based on samples of some noise scenarios. A scenario consists in a sample of the noises along the whole time horizon. Note that such a resolution method is stochastic in itself. Indeed, the solution given by such a method is a random variable in the sense that it depends on the scenarios sampled and used in the resolution. Therefore, when studying the performance of an algorithm based on such an approach, a variance term naturally appears, representing the sensitivity of the solution regarding these scenarios. This variance term is added to the squared bias term which arises from the necessity to approximate the solution in feedback by a function in a predefined space. Note that it is also present in the performance evaluation of deterministic techniques, like Dynamic Programming [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF] for instance.

For stochastic optimal control problems, it is common to represent the diffusion of "likely futures" using a scenario tree structure, leading to so-called multistage stochastic programs. This kind of representation goes back to [START_REF] Dantzig | Linear programming under uncertainty[END_REF] and has been widely studied within the Stochastic Programming community (see Prékopa, 1995;Shapiro and Ruszczynski, 2003, for a broad overview of this approach). This methodology consists in discretizing the probabilistic structure of the problem, then rewriting the constraints and objective function of the original problem on this discrete structure and finally solving it using a well-suited mathematical programming technique. In other words, the problem is solved for a particular sample of the random variables (a scenario tree) and the solution is hence a random variable in itself. We have to keep this fact in mind while evaluating the error. The main interest of such a methodology is that it leads back to a deterministic problem, on which we can apply classical tools of numerical optimization [START_REF] Bonnans | Numerical Optimization, Theoretical and Practical Aspects[END_REF]. For instance, one may then try to solve large-scale problems using decomposition techniques (see [START_REF] Carpentier | Stochastic optimal control and decomposition-coordination methods[END_REF][START_REF] Higle | Stochastic Decomposition: A Statistical Method for Large Scale Stochastic Linear Programming[END_REF]Shapiro and Ruszczynski, 2003, Ch. 3).

Despite the benefits of this methodology, it has been already observed that, in order to obtain a given precision, the number of scenarios needed to build the scenario tree has to grow exponentially with the time horizon of the problem. One of the aims of this paper is to highlight and to quantify this practical observation on an example. Thus, we start in Section 1 by introducing the Mean Squared Error, that is the distance between the strategy obtained using scenario trees and the (supposedly unique) optimal strategy. This is the indicator we use to estimate the performance of the method. Then, in Section 2, we study on a numerical example the relation linking the Mean Squared Error to the number of scenarios used for the discretization. [START_REF] Shapiro | On complexity of multistage stochastic programs[END_REF] obtains similar conclusions by working on the performance regarding the objective function and using large deviations tools.

We show that this negative observation is to be attributed to the scenario tree approach rather than to be considered as a feature of the original problem. Indeed, when using the particle method [START_REF] Dallagi | Méthodes particulaires en commande optimale stochastique[END_REF][START_REF] Carpentier | Particle Methods for Stochastic Optimal Control Problems[END_REF], in Section 3 on the same example, we show that the number of scenarios needed to obtain a given precision does not grow when the time horizon gets longer. This methodology, which mixes ideas from Stochastic Programming and Dynamic Programming, is a variational technique which consists in writing the optimality conditions of the optimization problem first, and then solving them using sampling techniques. This is a gradient-like method that builds an adaptive mesh over the state space as gradient iterations progress towards the solution. This mesh aims at discretizing the space in the regions mostly visited by the state vector at the optimum. According to the results obtained on the example we present here, the algorithm seems to be well-suited for multi-stage problems. However, with this approach, it seems that the difficulties arising from the dimension of the state space remain; this point is discussed in Section 4.

Mathematical formulation

Consider a controlled dynamic system, affected by random variables called noises. We aim to find strategies that minimize some expected cost over a finite time horizon, while satisfying a number of constraints. We suppose that the problem has a unique solution and propose to compare the approximate strategies using the Mean Squared Error (MSE).

1.1. The problem. We denote by T the finite time horizon of the problem and consider discrete time {0, . . . , T }. Let (Ω, A, P) be a probability space. Three kinds of random variables on this space are involved in the problem 1 :

• the state 2 X = (X t ) t=0,...,T , the values of which lie in a finite-dimensional vector space X; • the control U = (U t ) t=0,...,T -1 , the values of which lie in a finite-dimensional vector space U; • the noise W = (W t ) t=0,...,T , the values of which lie in a finite-dimensional vector space W, equipped with the σ-field W and the probability P W , which is the transport of the probability P by W .

Since some of the numerical methods we use in the following are gradient-based, it is natural to assume that the state and control lie in a Hilbert space. Hence we assume all three random variables are in L 2 (Ω, A, P). The noise variable W is given, i.e. its probability law is known, while the state and control variables are optimization variables. We here suppose that the information structure has perfect memory: at each time step t, the noise W t is observed and kept in memory, in such a way that the decision at time t is based on the knowledge of (W 0 , . . . , W t ). Hence we introduce the σ-field F t = σ{W 0 , . . . , W t } that represents the information available at time t, the associated filtration (F t ) t=0,...,T , and require the control at time t to be measurable with respect to F t . This constraint will be written: U t F t . At time 0, we assume that the value of the state is the realization of some random variable W 0 . Then, at each time step t, based on the available information, a decision U t is taken, a cost C t (X t , U t ) is incurred and the state evolves according to X t+1 = f t (X t , U t , W t+1 ). At the end of the time period, a final cost V (X T ) is incurred.

The aim is to find a pair (X, U ) that minimizes the expected sum of the costs over the time horizon while satisfying the dynamics and measurability constraints.

1 Throughout this paper, random variables will be denoted by bold letters (e.g. W ∈ L 2 (Ω, A, P; R p ).

2 In the Stochastic Optimal Control framework, the terminology "state" has a rather precise sense: this is the minimal information that completely sums up the past of the system in order to compute its future behaviour and cost function knowing all future inputs (see Powell, 2007, Def. 5.4.1). In that sense, the variable X can only be considered as the state variable when assumption 1 is made later on in §1.2. Until then, our use of this terminology is somewhat abusive.

Mathematically speaking, the problem we consider may be formulated as follows:

min

X,U E T -1 t=0 C t (X t , U t ) + V (X T ) , (1a) s.t. X t+1 = f t (X t , U t , W t+1 ) , ∀t = 0, . . . , T -1, (1b) X 0 = W 0 , (1c) U t F t , ∀t = 0, . . . , T. (1d)
Equations ( 1b) and (1c) describe the dynamics of the state variable and are P-almost sure relations. Equation (1d) states the information structure of the problem: it is the measurability constraint that forces the decision to be a function of all the past noises, and to be independent of future noise values. Hence we refer to this relation as the non-anticipativity constraint.

1.2. Mean Squared Error (MSE) for strategies. Now suppose that Problem (1) has a unique solution denoted by (X * , U * ). Suppose that we also have a numerical method, based on a scenario discretization, that gives us an approximate solution of (1) that we denote by (X ♯ , U ♯ ). Since this approximate solution given by the numerical method depends on, say, N samples, the pair (X ♯ , U ♯ ) lies in the probability space associated with these samples, that is (W, W, P W ) ⊗N . For instance, scenario tree-based methods compute values of the optimal decision for a finite number of state values that depend on the samples used for building the tree structure.

We make the following assumption:

Assumption 1. Random variables W 0 , . . . , W T are independent.
Under this assumption, the problem lies in the framework of Dynamic Programming (DP). In particular, we know that optimal controls can be expressed as feedback functions depending only on the state variable x. In other words, there exists some sequence of functions γ * such that U * t = γ * t (X * t ). In order to be able to compare the optimal and approximate solutions, we need to produce a solution of the same nature as the optimal solution, that is a feedback function such as γ * out of (X ♯ , U ♯ ). We hence suppose that we have such an interpolation/regression operator that gives us a feedback function Γ ♯ = (Γ ♯ t ) t=0,...,T out of the approximate solution (X ♯ , U ♯ ). More details on the interpolation/regression operators we use are given in §2.2. Note that the approximate strategy Γ ♯ = (Γ ♯ t ) t=0,...,T -1 is a random variable that lies in the samples probability space introduced earlier, that is (W, W, P W ) ⊗N .

We are interested in evaluating the efficiency of our numerical method for solving Problem (1) by computing the distance between the approximate strategy Γ ♯ and the optimal one γ * . Since those functions take values over the state space, we need to define some measure on this space. In our context, the measure that seems the most natural is the one corresponding to the density of the optimal state X * . Indeed, using this measure, we allocate more weight to regions which are often "visited" by the optimal state and we do not take into account decisions in regions that are never "visited" by the optimal state, since those decisions are unlikely to be used in practice. We denote by µ * t the density of the optimal state at time t. We can now introduce the MSE associated with strategy Γ ♯ , which is our performance indicator in this study:

MSE E T -1 t=0 X γ * t (x) -Γ ♯ t (x) 2 µ * t (x) dx .
Note that the expectation in the MSE is taken over all possible samples; in other words, it lies on the samples probability space (W, W, P W ) ⊗N . The MSE is thus the expected distance between the optimal strategy and the approximate strategy.

Recall that this expectation is taken with respect to the drawings used by the algorithm to produce strategy Γ ♯ . It is common to decompose the MSE using the socalled variance and squared bias. For this purpose, we introduce the expectation γ ♯ of the approximate feedback function Γ ♯ . Then

MSE = T -1 t=0 X γ * t (x) -γ ♯ t (x) 2 µ * t (x) dx Squared bias + E T -1 t=0 X γ ♯ t (x) -Γ ♯ t (x) 2 µ * t (x) dx Variance . (2)
The squared bias is the square of the distance between the optimal strategy and the expected approximate strategy, while the variance is the expected square of the distance between the approximate strategy and the expected approximate strategy. Both are real values.

Remark 1. Another relevant choice as a performance indicator would have been the following:

(1) Having the approximate strategy Γ ♯ , build the"true" state and control obtained while simulating this strategy, that is the pair (X † , U † ) that satisfies the dynamics (1b) and (1c):

X † 0 W 0 , X † t+1 f t X † t , U † t , W t+1 , with U † t Γ ♯ t (X † t ), ∀t = 0, . . . , T -1.
(2) Since the optimal decision U * is a random variable on (Ω, A, P) and the decision U † is a random variable that lies on the tensor product of (Ω, A, P) and the samples probability space (W, W, P W ) ⊗N , compute the expectation:

E U † -U * 2 U ,
where the expectation lies on the tensor product of the probability spaces (Ω, A, P) and (W, W, P W ) ⊗N . Note that this quantity may also be written

E( U † -U * 2 L 2 (Ω,A,P) )
, where the expectation now lies on the samples probability space only. This quantity is another relevant performance indicator since it is positive and it equals zero if and only if U † equals U * almost surely. Moreover, we note that the choice of U † for evaluating the performance of the numerical method is the one made by [START_REF] Shapiro | On complexity of multistage stochastic programs[END_REF], even though the comparison is performed using the objective function rather than the strategies. 1.3. Computing the MSE. For computing the squared bias and the variance as defined in Equation (2), we must face two main issues.

(1) At each time step t, we must compute integrals over the whole state space with respect to the density µ * t . In the following example, we are able to explicitly compute the optimal control strategy. We can then numerically integrate the Fokker-Planck equation to derive the density of the optimal state on a sufficiently dense grid and perform the integration by quadrature techniques. However, this technique suffers for the same curse of dimensionality as DP. For this reason we choose to perform the integration using a Quasi-Monte Carlo technique. This sampling method aims at distributing the samples associated with density µ * t equally over the state space using so-called low-discrepancy quasi-random sequences. This provides a very efficient convergence speed for the numerical computation of expectations. There exist numerous methods of Quasi-Monte Carlo sampling (see [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF], for a survey of these methods).

(2) We must be able to compute the expectation in the variance term. This expectation is taken over the samples used during the algorithm to produce the approximate strategy, e.g. the scenarios in the case of a scenario tree technique. We approximate this expectation by Monte Carlo, by performing the experiment a great number of times independently one from another.

In our case, 10 4 experiments were enough to precisely evaluate the variance. We are now able to estimate the efficiency of numerical methods using the MSE. We point out the relation between the error and the parameters of the methods, namely the number of scenarios used to discretize randomness.

Scenario tree-based methods

In the context of discrete time and finite horizon stochastic optimal control problems, a popular and widely used resolution technique consists in approximating the information structure of Problem (1) using scenario trees. This methodology has been widely studied by the Stochastic Programming community. We do not detail this method and refer to Shapiro and Ruszczynski (2003) for a more in depth presentation of Stochastic Programming.

2.1. Brief presentation. Because of the measurability constraints in Problem (1), and since there is no reason for the noise probability densities to have finite support, decision variables X and U are in general infinite-dimensional. Hence, except for some very particular cases where an explicit solution can be found, solving the underlying optimization problem directly is intractable. Scenario trees provide a way to approximate filtrations using noise samples to produce finite support approximations.

Let N be the set of all nodes in the tree, R be the set of the root nodes and L the set of the leaves. Building a scenario tree consists in defining the following functions:

• the time function θ : N → {0, . . . , T } which, with every node, associates the corresponding time step and the corresponding multi-application θ -1 which, to every time step t ∈ {0, . . . , T }, associates the set of all nodes at time t; • the function ν : N \R → N \L which, with each non-root node, associates the preceding node in the tree; • the weight function π : N → [0, 1] which, with every node in the tree, associates the probability to go through that node (the sum of all π(i), with i ∈ θ -1 (t) is equal to 1, for each t); • the multi-application F = ν -1 which, with every node i in the tree, associates the set of its successors (we impose the convention that F (i) = ∅, ∀i ∈ L); • the function F + which, with each node i, associates the whole sub-tree starting from i:

F + (i) = F (i) ∪ F 2 (i) ∪ • • • ∪ F T -θ(i) (i).
We consider a regularly branching tree, i.e. with a constant branching factor denoted by n b , so that card(F (i)) = n b , for every node i ∈ N \L. At time t = 0, we draw a n b -sample of the random variable W 0 , that we denote by (W ′i ) i∈θ -1 (0) . For each root node i, we draw at time t = 1 a n b -sample of the random variable W 1 3 , that we denote by (W ′j ) j∈F (i) . Samples of the random variable W 1 are drawn independently root node per root node. This procedure continues until final time t = T , so that we have n T +1 b leaves to the tree. Hence one can note that this can reveal a heavy procedure when T becomes large, since we need to draw n T +1 b samples of W T .

With this material we can now define a discretized version of Problem (1): min

X ′ ,U ′ i∈N \L π(i) • C θ(i) X ′i , U ′i + i∈L V X ′i , (3a) s.c. X ′i = f θ(ν(i)) X ′ν(i) , U ′ν(i) , W ′i , ∀i ∈ N \R, (3b) X ′i = W ′i , ∀i ∈ R. (3c)
The expectation in the objective function (1a) has been replaced in (3a) by a discrete weighted sum of costs on the tree nodes. Equations ( 3b) and (3c) correspond to Equations ( 1b) and (1c). The only constraint that seems to be missing is Equation (1d). Actually, it is now reflected in the structure of the tree. Indeed, starting from any node in the tree, there is only one possibility for going back to the root (we know exactly the values of all past noises), but there is a whole sub-tree that goes to the set of leaves, corresponding to the final time step (we only know the laws of future noises). In other words, the non-anticipativity constraint is coded in the structure of the tree.

In most cases, the scenario tree is not drawn directly. One draws a number of scenarios independently from one another, and then one builds a tree structure. Building such a tree in an efficient way (regarding the implied discretization error) is a challenging task. There exists a large literature on the way one can build scenario trees in a reasonable way [START_REF] Pflug | Scenario tree generation for multiperiod financial optimization by optimal discretization[END_REF][START_REF] Heitsch | Scenario Reduction Algorithms in Stochastic Programming[END_REF], as well as on stability of the underlying optimization problem regarding this discretization [START_REF] Heitsch | Stability of multistage stochastic programs[END_REF]. We will not focus on this point here. As it has already been mentioned, we suppose that the tree has been built using a branching factor n b , which does not depend on time, leading to n T +1 b leaves. Recall that this procedure requires the drawing of N n T +1 b different scenarios.

2.2. Case study. We apply this methodology to a simple instance where Problems (1) and (3) can be solved analytically. Then, we apply the methodology described in §1.3 in order to compute the MSE and to point out its relation with the number of scenarios N used to discretize the original problem. Let ε be some positive real number, the problem we consider is the following: min

X,U E ε T -1 t=0 U 2 t + X 2 T , (4a) s.c. X t+1 = X t + U t + W t+1 , ∀t = 0, . . . , T -1, (4b) X 0 = W 0 , (4c) U t F t . (4d)
The state and control variables lie in L 2 (Ω, A, P; R). Noise variables W 0 , . . . , W T are i.i.d. random variables with uniform law on [-1, 1]. One easily shows that the optimal control for Problem (4) reads U * t = γ * t (X * t ) with:

(5)

γ * t (x) = - x T -t + ε , ∀x ∈ R.
We now solve the problem on the tree. On each node i, one has a noise sample W ′i and solving the problem on the tree leads to values X ′i and U ′i for the state and control, respectively. One easily shows that:

U ′i = - X ′i + 1 n T -θ(i) b j∈F + (i) W ′j T -θ(i) + ε , ∀i ∈ N \L.
We now need to build a feedback function from these state and control values. This can be performed using a interpolation/regression operator. We here use the simplest one, which is the nearest neighbour interpolation operator (see Gersho and Gray, 1992, §10.4, for a rigorous study of this operator and its possible implementations). At each time step t, we build the Voronoi diagram of the set of state values {X ′i } i∈θ -1 (t) . We denote by C ′i the cell for which X ′i is the center. This cell is a random variable in itself, because it depends on the random variables drawn to build the scenario tree. On this cell, the feedback function will be constant, equal to U ′i . Thus we obtain the following strategy at time t:

(6)

Γ ♯ t (x) = i∈θ -1 (t) U ′i 1 C ′i (x) .
Remark 2. The fact that the Voronoi cells C ′i are random variables, depending on the scenarios, makes the theoretical study of the error associated with strategy Γ ♯ intricate. Indeed, in order to compute the expected approximate strategy γ ♯ , as well as to evaluate the variance term in Equation ( 2), we have to compute expectations over all possible scenario drawings, leading through Equation ( 6) to expectations over all possible Voronoi cells. This is generally a challenging task. This is the reason why we here appeal to numerical experiments to study the behaviour of the error with respect to the number of scenarios.

The results for this approach are presented in Figure 1 for T = 4 and a branching factor of n b = 3. The approximate strategy Γ ♯ 0 for the first time step is hence a piecewise constant function with 3 pieces. On the right-hand side of Figure 1, we draw 4 samples of this strategy (dashed curves) to emphasize that the scenario-tree method is stochastic: each sample Γ ♯ 0 is derived from a different scenario tree. We also draw the averaged strategy γ ♯ 0 over the 10 4 scenarios (dotted curve). One can remark that even if sample strategies are not continuous functions, the average strategy seems to be continuous, even smooth. On the left-hand side of figure 1, we draw the exact strategy (solid curve) obtained using equation ( 5) and the same average approximate strategy as on the right-hand side (dotted curve).

We can observe both the variance and the bias terms of the MSE on Figure 1: the variance term is the average L 2 -distance between the dashed curves (of course if there were 10 4 of them) and the dotted curve on the right hand-side, whereas the bias term is the L 2 -distance between the solid curve and the dotted curve on the left hand-side. From Equation (4c), the density of the optimal state, under which we compute the L 2 -distances cited above, equals the density of W 0 , which follows a uniform law on [-1, 1].

We are now able to compute the MSE using the protocol described in §1.3, for several values of the branching factor. We draw in Figure 2 the evolution of both the squared bias and the variance with respect to the branching factor, for the strategy at each time step (here T = 4). For computational time reasons, we were not able to perform this experiment for strategies at the third and fourth time steps

Γ ♯ 0 Γ ♯ 1 Γ ♯ 2 Γ ♯ 3 10 0 10 1 10 2 10 -4 10 -2 10 0 n b variance Γ ♯ 0 Γ ♯ 1 Γ ♯ 2 Γ ♯ 3 Figure 2.
Squared bias and variance for scenario trees On this figure which uses a logarithmic scaling on both axes, one can observe that the variance term on the right-hand side obviously dominates the squared bias term on the left-hand side. Moreover, while the convergence rate for the variance seems to be constant with respect to time, the bias decreases with a highest rate as time grows. This is natural since the tree has more nodes in the last time steps than in the first ones to represent the approximate strategy.

A careful inspection of the variance shows that the MSE appears experimentally to have a convergence rate of n -1 b for every time step. In other words, the number of scenarios needed to obtain a given precision varies exponentially with respect to the time horizon: it is proportional to n T b . Still, recall that in order to have a branching factor n b equal to 10 2 , since we have a 4-period problem, we need to draw 10 8 scenarios!

Particle methods

We here propose a similar study using another numerical method based on scenarios: the particle method. This variational technique consists in first writing the optimality conditions for the problem. Then, scenarios are used to compute the expectations that lie in the optimality conditions. We briefly present the idea of this method and then study its error with respect to the number of scenarios used for sampling. A detailed description can be found in the work by [START_REF] Carpentier | Particle Methods for Stochastic Optimal Control Problems[END_REF].

3.1. Brief presentation. The particle method consists in solving first-order optimality conditions that have been approximated through sampling. We first present the optimality conditions and then describe the sampling procedure. The resulting approximate optimality conditions are solved numerically using a gradient technique.

We suppose that the cost function and the system dynamics are differentiable with respect to both the state and the control, and that their derivatives are Lipschitz-continuous.

We suppose that Assumption 1 is still in force. In this context, Carpentier et al (2009, Theorem 2.6) state that, if a solution (X * , U * ) to Problem 1 exists, then, for every time step t = 0, . . . , T , there exists a random variable Λ * t ∈ L 2 (Ω, A, P; X), called the adjoint state, such that the following conditions hold:

X * 0 = W 0 , (7a) X * t+1 = f t (X * t , U * t , W t+1 ) , (7b) Λ * T = ∂V ∂x (X * T ) ⊤ , (7c) Λ * t = ∂C t ∂x (X * t , U * t ) ⊤ + E ∂f t ∂x (X * t , U * t , W t+1 ) ⊤ Λ * t+1 X * t , (7d) 0 = ∂C t ∂u (X * t , U * t ) ⊤ + E ∂f t ∂u (X * t , U * t , W t+1 ) ⊤ Λ * t+1 X * t . (7e)
Note that these optimality conditions are specific to the Markovian case and are not obtained in a straightforward manner. The authors first write first-order optimality conditions conditionally to the filtration F t . Then Assumption 1 is used at each time step to recursively replace the conditional expectations with respect to F t by conditional expectations with respect to the optimal state X * t . Equations (7a) and (7b) simply recall the state dynamics of the system. Equations (7c) and (7d) introduce a backwards relation on the adjoint states. Actually, the adjoint state at time t may be seen as the sensitivity of the optimal cost from time t to the end of the time horizon with respect to the state variable at time t. Hence there exists a relation between Equations (7c) and (7d) and the classical DP equation. Finally, Equation (7e) states that the derivative of the optimal future cost at a given time step with respect to the control at the same time step equals zero at the optimum. This is the classical first-order optimality condition when no additional constraints (apart for the non-anticipativity constraints) are imposed to the control.

Note that the adjoint state Λ * t+1 in Equations ( 7d) and ( 7e) is, by construction, measurable with respect to the state variable X * t+1 = f t (X * t , U * t , W t+1 ). Because of Assumption 1, the conditional expectations in Equations ( 7d) and (7e) are actually expectations that are only supported by the random variable W t+1 . One has to keep this in mind when discretizing these conditions.

We now want solve these conditions numerically using sampling techniques and gradient methods. Suppose we have N samples W ′1 , . . . , W ′N , called scenarios or noise particles, for the random variable W . After the k-th iteration, suppose we have N samples for the current state X (k) , control U (k) , and adjoint state Λ (k) . In order to compute conditions (7d) and (7e), we need to evaluate expectations of functions involving Λ (k) t+1 , knowing the value X (k) t of the state at time t, for every particle k = 1, . . . , N . We have two ingredients to perform this operation:

The algorithm stops when all gradient particles are sufficiently small, ensuring that condition (8e) is almost fulfilled, for every i = 1, . . . , N .

Remark 3. Note that, at each iteration, we have N samples (or particles) for the state and control at every time step. Using a regression operator, we are able to compute a feedback function such as the one introduced in Section 1.2, denoted by Γ ♯ . This feedback function associates a decision with every possible state of the system and this is the quantity on which we base our error analysis.

3.2. Case study. using the same example as in §2.2, we represent in the righthand side of Figure 3 some samples of the approximate strategy (dashed curves) obtained by the particle method using 3 4 scenarios for discretization. Observe that these are piecewise constant functions with 3 4 pieces, instead of only 3 pieces when using scenario trees. Indeed, using the particle method one has the same number of nodes at every time step. In other words, all the scenarios are used at each time step. In the left-hand side of the same figure we draw the average approximate

-0.5 0 0.5 -0.2 -0.1 0 0.1 0.2 state control γ * 0 γ ♯ 0 -0.5 0 0.5 -0.2 -0.1 0 0.1 0.2 state control γ ♯ 0 Γ ♯ 0 Figure 3.
Strategies for particle methods strategy (dotted curve) and the optimal strategy (solid curve). The bias term in the MSE is still the L 2 -distance between the dotted curve and the solid curve, while the variance term is the average of the L 2 -distances between the dashed curves and the dotted curve. Using the same number of scenarios as we did for for scenario trees in §2.2, we observe that both the bias and variance look much smaller to those observed in Figure 1.

We can now compute the MSE for several values of the number N of scenarios used for discretization. We observe in figure 4 how the squared bias and the variance associated with each strategy, i.e. for each time step, behave when the number of scenarios grows. Note that the squared bias, on the left-hand side is dominated by the variance, on the right-hand side. A careful inspection of the latter shows that the term leads experimentally to a convergence rate that is slightly less than N -1 . When comparing with the results obtained using scenario trees, one should notice that the x-axes do not refer to the same quantities in Figures 2 and4. Hence, we here observe a convergence rate in N -1 , N being the number of scenarios, while we observed a convergence rate of n -1 b = N -1 T in the case of scenario trees, which is much smaller.

Unlike scenario trees, the error associated with particle methods does not depend on the time horizon on this example: it is the same at every time step.

The state space dimension issue

Most numerical methods that aim at solving stochastic optimal control problems encounter difficulties when the dimension of the state space becomes large. Since we

10 0 10 1 10 2 10 -3 10 -2 10 -1 10 0 N squared bias Γ ♯ 0 Γ ♯ 1 Γ ♯ 2 Γ ♯ 3 10 0 10 1 10 2 10 -3 10 -2 10 -1 10 0 N variance Γ ♯ 0 Γ ♯ 1 Γ ♯ 2 Γ ♯ 3 Figure 4
. Squared bias and variance for particles.

are looking for strategies, that is functions that map a decision from every possible state of the system, even just storing such functions becomes challenging when the dimension of the state space is large. In the stochastic optimal control community framework, this difficulty is known as the curse of dimensionality.

Note that particle methods are not prone to the same curse of dimensionality as DP, for instance. While DP encounters difficulties because it requires to build a somewhat regular grid to cover the state space, particle methods are designed to construct an adaptive mesh over the state space: the state particles are built so as to concentrate in regions that are often visited by the optimal state. As such, regardless of the dimension of the state space, whenever the dispersion of the optimal state is small, the particle method shall produce satisfying results.

In the example we here study, the density of the optimal state does not seem to have this "tightness property". On a simple extension of Problem (4) to a twodimensional state space case, we observe in Figure 5 how the squared bias and the variance associated with each strategy, i.e. for each time step, behave when the number of scenarios grows. Note that the squared bias, on the left-hand side, now prevails over the variance term, on the right-hand side. A careful inspection of the latter shows that the term leads experimentally to a convergence that is slightly less than N -0.5 . The variance itself also seems to decrease at a lower rate than in the one-dimensional case. Both of these observations indicate that the performance of the method could benefit from a better regression operator than the nearest neighbour that we used here. Note that, for computational time reasons, we cannot produce such experimental results for higher dimensions. Hence we are not able to make a thorough study of the relation between the dimension of the state space and the convergence rate of both the squared bias and the variance.

On the other hand, when dealing with large state space dimensions, a clever idea is to make use of decomposition schemes in order to replace the solving of a high dimensional problem by the iterative solving of several smaller problems. This is indeed a common and powerful technique when using scenario trees (see [START_REF] Carpentier | Stochastic optimal control and decomposition-coordination methods[END_REF][START_REF] Higle | Stochastic Decomposition: A Statistical Method for Large Scale Stochastic Linear Programming[END_REF]Shapiro and Ruszczynski, 2003, Ch. 3). Both the improvement of the regression operator and the use of a decomposition scheme are the subject of ongoing research.

Conclusion

We presented a numerical study concerning the comparison of two scenariobased approaches for solving stochastic optimal control problems with a discrete finite horizon. In the first section, we introduced the performance indicator that was used to compare both methods, namely the Mean Squared Error (MSE).

The first approach we considered was scenario tree modeling. We observed that the number of scenarios needed to obtain a given accuracy grew exponentially with the time horizon, making the implementation for multi-stage problems hardly tractable.

The second approach we considered was particle methods. In this case, we observed that the associated MSE does not depend on time. Thus, this recent approach seems to be well-suited for multi-stage problems. Unfortunately, it seems that its performance deteriorates rapidly when the dimension of the state space increases. However, this method is quite new and as such has received little attention. Its implementation and hence performance would certainly benefit from computational studies. Future works will be concerned with this topic.

Figure 1 .

 1 Figure 1. Strategies for scenario trees

Figure 5 .

 5 Figure 5. Squared bias and variance for particles in dimension 2.

This tree building procedure, often called conditional sampling, usually requires to draw samples of W 1 knowing the value of the preceding noises, that is W 0 . Because of Assumption 1, this is not needed here.

(1) at the optimum, we know that for every time step t = 0, . . . , T -1, the adjoint state Λ * t+1 is a function of the state X * t+1 = f t (X * t , U * t , W t+1 ), where U * t is measurable with respect to X * t , and W t+1 is independent of X * t ;

(2) we have N samples of all the random variables involved, namely X (k) , U (k) , Λ (k) and W . All we need is to define a regression operator, based on the current state and adjoint state samples, that associates with every value of the current state an estimate for the corresponding value of the adjoint state. At time t, we denote this regression operator by Λt . Using this material, we are now able to write the discretized optimality conditions, namely the initial condition on the state for every particle i = 1, . . . , N :

the state dynamics for every time step t = 0, . . . , T -1, and for every particle i = 1, . . . , N :

the final condition on the adjoint state, for every particle i = 1, . . . , N :

the backwards adjoint state dynamics, for every time step t = T -1, . . . , 0 and every particle i = 1, . . . , N :

and the stationarity condition, for every time step t = T -1, . . . , 0 and every particle i = 1, . . . , N :

Note that in both equations (8d) and (8e), the sums only involve the index j, and hence the noise variable only. This corresponds to the fact that the conditional expectations in (7d) and (7e) with respect to X * t are in fact expectations. The algorithm then consists in solving these conditions using a gradient method. Thus, each iteration k of the particle method consists in three steps:

(1) integrate the N state dynamics (8a) and (8b) using the current control particles U ′1,(k) , . . . , U ′N,(k) ; (2) integrate the N backwards adjoint state dynamics (8c) and (8d) using the current state and noise particles; (3) compute the N "gradient particles" G ′1,(k) , . . . , G ′N,(k) (the right-hand side of equation ( 8e)) and update every control particle i using a gradient step:

U ′i,(k+1) = U ′i,(k) -ρG ′i,(k) .


Axel Freyn

Serge Florens

Optimal broadening of finite energy spectra in the numerical renormalization group: application to dissipative dynamics in two-level systems

Keywords:

Numerical renormalization group (NRG) calculations of quantum impurity models, based on a logarithmic discretization in energy of electronic or bosonic Hamiltonians, provide a powerful tool to describe physics involving widely separated energy scales, as typically encountered in nanostructures and strongly correlated materials. This main advantage of the NRG was however considered a drawback for resolving sharp spectral features at finite energy, such as dissipative atomic peaks. Surprisingly, we find a bunching of many-body levels in NRG spectra near dissipative resonances, and exploit this by combining the widely-used Oliveira's z-trick, using an averaging over few discrete NRG spectra, with an optimized frequency-dependent broadening parameter b(ω). This strategy offers a tremendous gain in computational power and extracts all the needed information from the raw NRG data without a priori knowledge of the various energy scales at play. As an application we investigate with high precision the crossover from coherent to incoherent dynamics in the spin boson model.

A general hallmark of many-particle interaction, as found in a variety of condensed matter systems such as nanostructures and strongly correlated materials, lies in the presence of several energy scales, possibly widely separated from each other due to renormalization and dynamical effects. Two well-studied examples found in electronic systems concern the Kondo effect for magnetic impurities in metals and Fermi liquids in proximity to a Mott insulating phase, two instances where low-energy quasiparticles emerge below a typical temperature which is quite reduced from the bare Fermi energy [START_REF] Hewson | The Kondo Problem to Heavy Fermions[END_REF] . These lowlying excitations do however coexist with higher energy atomic levels, also broadened and displaced in a strong manner from their bare atomistic values due to the dissipation brought by the electronic environment. Such complex physical effects, taking place on a broad range of energies, entail great practical difficulties for most direct numerical approaches. These are partially lifted using Wilson's original idea of the logarithmic discretization 2,3 , as implemented in numerical renormalization group calculations (see Ref. 4 for a recent review). This technique has been improved in the last twenty years to calculate static and dynamic quantities both for fermionic 5 and bosonic quantum impurity models 6 . Important practical applications until now involve the calculation of transport in the Kondo regime for Kondo alloys and artificial quantum dots 7 , as well as the accurate description of the zero-temperature Mott transition by combining NRG 8 with Dynamical Mean Field Theory (DMFT) 9 . More generally, exponentially small energy scales are also found in the vicinity of quantum critical points 10 , so that impurity models provide a simplified testbed for the theory of quantum critical phenomena [START_REF] Vojta | [END_REF] . Again the NRG is the ideal technique for studying such impurity quantum phase transitions 4 , with potential implications for artificial nanostructures 12,13 . Despite these successes, the foundation of NRG on a logarithmic discretization mesh implies a loss in accuracy for the high energy spectral features, which has plagued most calculations so far. Not only are atomic-like excitations physically observable, using photoemission or tunneling spectroscopies, but they may also be intimately tied to the low-energy excitations. Such interesting behavior occurs in the vicinity of the Mott metal-insulator transition 9 , where the self-consistently determined electronic environment in the DMFT picture shows electronic states violently redistributed over all energy scales. The numerical cost of converging the DMFT equations certainly requires efficient and accurate NRG calculations for the spectral functions, without a priori knowledge of the excitations involved. For this reason, the idea of averaging over N z realization of the Wilson chain, using the so-called z-trick 14 to fill in the missing spectral information (to be discussed later on), seems prohibitive for most practical calculations, and offers a limited gain for very narrow spectral structures (see Ref. 15 for a detailed study).

In this Rapid Communication, we show that the broadening procedure used to smoothen the discrete NRG data is one critical step for obtaining an optimal resolution at all energies. We henceforth develop a simple adaptive procedure where the standard broadening parameter b is taken to be frequency dependent. This choice is dictated by our surprising observation that the density of many-body NRG levels increases sharply as soon as narrow atomic resonances are encountered. Together with the usual z-trick, combining several but reasonably few NRG calculations, this extra b(ω)-trick allows to compute finite frequency spectra using as little as N z = 10 NRG calculations, in situations where large scale NRG z-averaging would be prohibitive. An improved broadening procedure leading to errors in spectral functions limited to few percents could constitute a further step in applying the NRG to a wider class of problems, such as DMFT+NRG calculations.

In order to explicitly demonstrate these ideas, we focus on the simplest quantum impurity model, namely the spin boson Hamiltonian:

H = ∆ 2 σ x + λ 2 σ z i (a † i + a i ) + i ω i a † i a i , (1)
that involves a two-level system, described by a quantum spin 1/2, and a bosonic bath a † i with continuous spectrum ω i of energies. ∆ is a transverse magnetic field, while the coupling constant λ controls the strength of longitudinal dissipation. The bosonic spectrum is assumed to be subohmic with bath exponent 0 < s ≤ 1 (this includes the well-studied ohmic case s = 1):

J(ω) ≡ i πλ 2 δ(ω-ω i ) = 2παω 1-s c ω s θ(ω)θ(ω c -ω) (2)
where ω c is a high-energy cutoff. For small values of the dissipation α, this model is known to exhibit coherent precession of the spin around the x axis at zero temperature. By increasing the coupling to the bath for values of s 1, Rabi oscillations are progressively damped 16 , before the two-level system localizes in one potential minima via a quantum phase transition 6 . When s 0 however, the phenomena of localization and decoherence occur in reverse order 17 . Despite its simplicity, this model embodies all the effects typical of strong correlations: low-energy scales are indeed dynamically generated near the localization transition, while discrete spin excitations due to the transverse field ∆ are deeply affected by the bosonic environment, that leads to broadening and frequency shift in the magnetic response

C(ω) = 1 2π +∞ -∞ dt C(t)e iωt (3)
with C(t) = 1 2 [σ z (t), σ z ] the spin autocorrelation function.

The implementation of the NRG procedure follows the standard route 2 as discussed in Ref. 6 for the extension to bosonic Hamiltonians. The bosonic band (2) is logarithmically discretized with the Wilson parameter Λ > 1, first on the highest energy interval [Λ -z ω c , ω c ], and then iteratively on successive decreasing energy windows [Λ -n-z ω c , Λ -n-z+1 ω c ] for n strictly positive integer. This also introduces the crucial z parameter (0 < z < 1) that is used to average over N z Wilson chains 14 , allowing to obtain better resolution on the finite energy states. The rest of the NRG follows Ref. 6, coupling iteratively the kept energy levels up to iteration n to the states living in the shell n + 1, and truncating the successive Hamiltonians to keep up with a manageable number of eigenstates. All subsequent calculations were performed with Λ = 2, N b = 8 kept bosonic states on the added bosonic "site", and N = 160 kept NRG states, ensuring good convergence. The resulting discrete spectra at successive NRG iterations are combined using the interpolation scheme proposed in Ref. 18 (see however Ref. 19 for a more rigorous implementation), leading to a set of z-dependent many-body energy levels ǫ a,z labelled by quantum number a. The spin-spin correlation function is thus readily obtained at zero temperature as:

C(ω) = 1 2N z a,z | 0, z|σ z |a, z | 2 δ (|ω| + ǫ 0,z -ǫ a,z ) (4)
where ǫ 0,z is the ground state energy, and these raw NRG data are displayed on Fig. 1, to be discussed below. Let us now discuss the general structure of the NRG spectra. For a single NRG calculation with a given value of z, the a priori energy resolution at the scale ω N = Λ -N -z is δω N = Λ -N -z -Λ -N -z-1 = (1 -Λ -1 )ω N , so that the resolution degrades at increasing energy. For obtaining smooth NRG spectra, the delta-peaks in (4) are usually broadened 4 at energy ω N on the same scale:

δ (|ω| -ω N ) → e -b 2 /4 ω N b √ π e -[log(|ω|/ωN ) 2 /b] 2 (5)
with b ≈ 0.7 typically. Combining N z NRG runs using the z-averaging in (1) allows in principle to improve by a factor N z the accuracy at high energy, since the broadening parameter may now be decreased down to b = 0.7/N z . This procedure actually faces two problems: i) for very sharp resonances (typically several orders of magnitude narrower that the natural high-energy cutoff), parallelizing N z ≫ 10 NRG calculations becomes too prohibitive, especially with the aim of DMFT simulations; ii) states far from the resonances can then be too much underbroadened, so that oscillations of period Λ/N z can become quite prominent. In this view, quantitative spectra can be accurately extracted from the NRG data only in the continuum limit, either with Λ → 1 or N z → +∞, as assumed in the previous literature 15,20,21 . We wish to show however that NRG provides much more information close to dissipative resonances than usually expected. For this purpose, we closely examine the raw NRG spectra for resonance widths as small as 10 -4 ω c , using reasonably few interleaved averaging N z = 20, much smaller that the naively needed 10 4 NRG runs. Fig. 1 shows that for frequencies far from the resonance located at ∆, the highest weight NRG states come in packets of N z peaks for a given Wilson shell. Surprisingly, the density of NRG eigenstates demonstrates a huge increase precisely at the peak value, so that information about the resonance width and height seems really encoded in the discrete results! The standard broadening b = 0.7/N z = 0.035 is however way too large to benefit from this effect, and the corresponding smoothed spectra are indeed completely inaccurate. For definite comparison, we have plotted the analytical result obtained from the expansion at weak dissipation 22 :

C(ω) = - 1 π Im ω -Σ(ω) ω[ω -Σ(ω)] -∆ 2 (6) Σ(ω) = dǫ π J(ǫ) ω -ǫ + i0 + = ωc 0 dǫ 2αω 1-s c ǫ s ω -ǫ + i0 + .
In order to judiciously exploit this unexpected finding, the straightforward idea we follow here is to adapt the broadening parameter b in (5) to the frequency dependence of the local density of NRG peaks. A very natural approach adapted to the NRG logarithmic discretization is to extract b(ω) from the logarithmic derivative of the integrated spectrum up to frequency ω:

b(ω) = b 0 2   q + d log ω 0 C d log ω -1 + q + d log ω +∞ C d log ω -1   (7)
where q ≃ 1 is regularization parameter, whose precise value is not sensitive to the final results, and b 0 sets the typical broadening at very low and very high frequencies (far from the atomic resonances). Note indeed that we average in this expression two frequency sweeps from ω = 0 and ω = +∞ respectively, in order to treat on an equal basis low and high frequency tails. Because the actual NRG data is fully discrete (4), we extract b(ω) recursively using equation ( 7) on the broadened NRG spectra. This procedure converges after few iterations to the results displayed on Figs. 1 and2. Given the small nu- merical effort devoted to generate these NRG data, the quality of final result is quite astonishing. The reason for this success is given on Fig. 2, where the b(ω) parameter takes values as small as 10 -4 at the resonance, while it increases drastically far away from it, naturally cancelling out a great part of the NRG oscillations due to the discretization of the Wilson chain. We emphasize that the only free parameter is the typical low-frequency broadening b 0 in (7), which is easily adjusted: we always find a range of b 0 values where the self-consistently broadened spectrum is relatively independent of the chosen b 0 . More generally, one can also avoid the use of Ansatz (7) by determining b(ω) from estimates of the local density of highest weight states in the raw NRG data, giving similar results. Large scale NRG simulations with N z ≫ 10 and fixed b, as performed e.g. in 15 , may thus benefit greatly from an adaptive broadening procedure. We also note that the high frequency tails are also faithfully reproduced using the b(ω)-trick. Finally, in situations where relatively broader resonances are present, we show that quantitative results can be obtained at lower computational cost, by decreasing significantly the number of z-averaging even as low as N z < 5, see Fig. 3.

In order to check that our broadening procedure is completely robust for the whole range of parameters in the model (1), we investigate the effect of increasing dissipation for several bath exponent values s. The quantitative comparison of the Rabi resonance to the analytical formula (6) in its domain of validity, namely α ≪ 1, gives us indeed confidence in the adaptative method. This is proved on Fig. 4 and5, which consider the same parameters as in Ref. 17, with important quantitative improvement. For the small value of s = 0.1 taken on Fig. 4 perfect matching of the NRG data on the resonance to the lowest order calculation in α shows that the dissipation mechanism does not care for the low-energy behavior of the spin dynamics, even at the quantum critical point 6,17 where the spin localizes for the values α c = 0.0071 and 0.105 with s = 0.1 and 0.5 respectively (∆/ω c = 0.1 and Λ = 2 here). Interestingly, the low-energy tails show significant deviations from the analytical result (6) even far from the quantum phase transition, so that lowest-order calculations 22 do not apply for the long-time dynamics. Improved resummation of the perturbation theory to all orders in α will be considered in a forthcoming work 23 . For the intermediate value s = 0.5 shown on Fig. 5, clear non-perturbative effects on the resonance are seen inbetween the weak coupling regime and the quantum critical point, so that perturbation theory in α breaks down. To conclude, we have investigated the properties of discrete NRG spectra near atomic states in the spin boson model, and found that a bunching of the many-body levels occurs whenever sharp resonances are encountered. An adaptive broadening scheme was proposed, showing a drastic improvement in computation power for the calculation of accurate spectral functions over the whole energy range. This procedure will certainly allow to take further advantage of the potentialities of the NRG in a wide range of physical situations, from ab-initio calculations for magnetic impurities in metals 24 to the difficult problem of simulating strongly correlated materials in the framework of the Dynamical Mean Field Theory 9 .

  FIG. 1: (Color online) Spin susceptibility C(ω) of the subohmic spin boson model at s = 0.1, ∆ = 0.1ωc and α = 0.000125. Raw data | 0, z|σz|a, z | 2 /[2Nz(ǫa,z -ǫ0,z)] are given as circles for Nz = 20 combined NRG calculations, solid line is the perturbative result (6), and the three dashed lines are the various NRG broadenings discussed in the text.

perturbative

  FIG. 2: (Color online) Comparison of various broadening procedures near the resonance for the parameters of Fig. 1: standard interleaved averaging with Nz = 20 and b = 0.7/Nz = 0.035 (dash-dotted line), b(ω)-trick with b0 = 0.12 (dashed line), and analytical result (full line). The adaptive broadening parameter b(ω) is also given (note the logarithmic scale on the right).

perturbative

  FIG. 3: (Color online) Similar comparison as done on Fig. 2, with stronger dissipation α = 0.02. Optimized broadenings with b(ω)-trick were done with Nz = 2, 20 and b0 = 0.6, 0.12 respectively.

FIG. 4 :

 4 FIG. 4: (Color online) Comparison of optimized NRG spectra (solid lines) and perturbative calculations (dashed lines) with bath exponent s = 0.1, ∆ = 0.1ωc and increasing values (arrow) of dissipation α = 0.002, 0.004, 0.007.

FIG. 5 :

 5 FIG.5:(Color online) Similar plot as on Fig.4, with bath exponent s = 0.5, ∆ = 0.1ωc and α = 0.02, 0.06, 0.10.


Valentin Gavan
email: valentin.gavan@insa-lyon.fr

Monika Woloszyn

Frédéric Kuznik

Jean-Jacques Roux

Experimental study of a mechanically ventilated double-skin façade with venetian blinds: A full-scale investigation in controlled environment

Keywords: double-skin façade, full-scale experiment, mechanical ventilation, venetian blinds

The aim of this article is to present results of an experimental campaign performed on a full scale facility provided with a double-skin façade. The behaviour of this architectural concept is tested under controlled climatic conditions. A summer case is scrutinised under different configurations: variation of the airflow through the double-skin façade and different angle of the solar shading device. This paper describes the experimental conditions, as well the test facility and the tested façade element. The results show the temperatures of the test cell and the façade and how they depend on the climatic conditions and the blind blade angles. One objective of this research was to measure and provide extensive data set detailing air and surface temperatures on the double-skin façade, together with airflow rates and air velocities. The experiments are fully described so that the results can be used for the validation of numerical models dealing with ventilated double-skin façades with venetian blinds.

Introduction

Double-skin façades are highly technological building components which are deployed especially for office buildings in order to satisfy special requirements. Although the double-skin façade (DSF) concept is not new, there is a growing tendency from the architects to use it. In new architectural projects, the DSF are designed to fulfil several envelope functions, such as thermal and acoustic insulation, optimization of natural lighting and improvement of ventilation system. Given that the study of DSF constructions has appeared during recent years, a satisfying level has not been yet achieved and the lack of knowledge is still present.

The DSFs were studied in the past years to improve their capacities headed for the energy performance simultaneously with preservation of the thermal and visual comfort. In addition to the envelope itself (façade), DSF have a second glazed layer (generally without a structural function) placed at a certain distance from the inner layer [START_REF] Streicher | Bestfaçade -Best Practice for Double Skin Façades[END_REF]. These two panes of glass act as an insulation between outside and inside enabling the air to circulate the cavities of the façade. The zone positioned between these two layers is named buffer zone, channel or air gap and is generally ventilated. Since the distance between the glass layers goes from a few centimetres to over one meter, this significant air volume can be used various ways according to climatic conditions, hours of occupation, orientation, type of construction and HVAC system. Its temperature is influenced by solar radiation, outside temperature, presence of shading systems.

Although numerous papers describe different configurations of DSF, experimental results under controlled climatic environment are rarely given. Measurements are very often carried out under real weather conditions where there is a little or no control of the key environmental factors [START_REF] Pasquay | Natural ventilation in high-rise buildings with double façades, saving or waste of energy[END_REF], [START_REF] Park | Real-time optimization of a doubleskin façade based on lumped modelling and occupant preference[END_REF], [START_REF] Stec | Symbiosis of the double skin façade with the HVAC system[END_REF], [START_REF] Corgnati | Experimental assessment of the performance of an active façade during actual operating conditions[END_REF]. All these experiments treat different façade configurations. Results of measurements are used for validation and enhancement of developed models. Generally, these results are employed to improve DSF performance at the design stage.

In [START_REF] Pasquay | Natural ventilation in high-rise buildings with double façades, saving or waste of energy[END_REF] several high-rise naturally-ventilated DSF on real buildings are monitored. In [START_REF] Stec | Symbiosis of the double skin façade with the HVAC system[END_REF] three main experiments are evaluated, one in a small-scale laboratory facility and two under real weather conditions. The small-scale test cell measurements are performed to validate a simplified DSF model without an influence of wind and solar radiation. Measurements under real weather conditions were used to validate the building model coupled with the DSF. In [START_REF] Corgnati | Experimental assessment of the performance of an active façade during actual operating conditions[END_REF] the authors evaluate extensively the operation of a DSF integrated with HVAC strategies. The air is drawn from the room and flows through the façade to an air treatment unit. The main inconvenient of this strategy is that during the winter, due to the occupants, condensation may appear on the outer glazing of the façade.

The major inconvenient of this large variety of experimental situations reside generally in the different climatic conditions, as well as building or DSF characteristics. In this case, aiming inter-building and/or inter-façade comparisons became very complicated.

Experimental set-up

The first part of this paper is dedicated to an extensive description of the experimental set-up. The second chapter deals with the results obtained during the experimental campaign. Finally, conclusions are given in the last part of the paper.

Minibat test cell

The study of thermal characteristics or influence of a system connected to a building structure, requires reliable and detailed experimental data. However, the consistency of such measurements can only be acquired through controllable experimental environment, situation constituting the essence of the climatic test cell called Minibat. Used intensively for the last years to validate zone and CFD models or to measure pollutant, Minibat cell was upgraded recently to a new potential, the assessment of a double-skin façade system.

The Minibat experimental facility is composed by two identical test cells (rooms), each one with the dimensions L × W × H = 3.1 m x 3.1 m x 2.5 m (Figure 1). Throughout measurements, the door between the tests cells was sealed and only the data from test cell 1 were considered, henceforth called test cell.

Five walls of the test cell are in contact with a thermal buffer zone with controlled temperature. The composition of the walls is given in Table 1. The thermophysical properties of the materials are summarized in Table 2. To create a homogeneous and regular volume of air inside the thermal buffer zone, a distribution air network is installed. The air diffusers are installed in the upper part of the thermal buffer zone and air outlets in the lower part. The air temperature inside the thermal buffer zone is entirely controlled with an accuracy of 0.5°C. The sixth face of the test cell, where the DSF is installed, is in contact with a weather generator capable of simulating the outdoor conditions (temperature and solar radiation). An air treatment unit controls the air temperature in the weather generator. The average air temperature can vary between 0°C and 40°C with an accuracy of 1.5°C.

Large scale solar simulator

The large scale solar simulator installed in the Minibat test cell provides a source of artificial solar radiation for irradiating the test façade surface (see Figure 1). The solar simulator is designed to test elements W × H = 3.1 m x 3.1 m. The lamps are physically separated from weather generator by a glass pane. Thus, the separation between solar and weather generator allows measurements with different scenarios (summer, winter, middle season configuration), without deterioration of the lamps characteristics. This environment that adjoins the front irradiated surface of the façade and the solar simulator corresponds to outdoor environment.

The solar radiation is simulated using twelve 1000W CSI lamp (a gaseous discharge lamp of metal halide type, CSI meaning compact source iodide). Figure 2 shows the comparison between the CSI spectrum and the solar spectrum [START_REF] Allard | Thermal Experiments of Full-Scale Dwelling Cells in Artificial Climatic Conditions[END_REF]. Then, the lamps are capable of producing a radiative spectrum close to the solar spectrum. The lamps are installed in three rows of four lamps in front of the weather generator (Figure 3). Every lamp is electrically controlled by a supervising unit which allows the functioning of different combination of lamps. During our tests all the lamps were switched on.

Double-skin façade

As a result of a literature review, the more conventional façade element comprising a single storey "boxwindow" type is selected in our investigation [START_REF] Loncour | Les double façades ventilées. Classification & Illustration des concepts de façades[END_REF]. In this respect, a full sized double-skin façade incorporating a sun-shading blind was constructed and installed in the Minibat test cell.

The double-skin façade, presented in Figure 4, is made of two aluminium frames: external and internal.

The internal is composed of one fixed part (W × H = 2.8 m x 2.3 m) while the external one is divided into two parts that can be separately opened. The aluminium frame (6 cm thick by 10 cm of width) is a specially designed 3 cell formation limiting thermal bridges. The glass panels, separated by a gap of 20 cm, are two simple clear glazing without a surface treatment. The glass area is W × H = 1.3 m × 1.93 m (2 panes) for the external pane and W × H = 2.8 m × 1.93 m (one pane) for the internal pane. The outer skin was designed to be opened and closed for easy installing the measurement devices and directing the solar blinds.

The DSF includes four openings for ventilation (W × H = 2.6 m x 0.04 m). They cover the whole width of the façade, at top and bottom of each of the two glass panels. Thus, various ventilation configurations of double-skin can be tested by simple obstruction of these openings.

The sun-shading devices, commonly used in practice and installed for these experiments, are venetian blinds placed in the middle of the channel. Since the sun-shading device is motor driven, the blinds can be oriented automatically during the tests, from the exterior of the façade. The orientation is representing an imposed angle (β) between 0° and 90°. The 0° angle corresponds to horizontal completely open blinds.

The 90° angle corresponds to completely closed blinds facing the solar simulator. Made of aluminium, tested sun-shading blades are 2.5 cm of width (total sun-shading system dimension is W × H = 2.8 m × 2.3 m), 0.21 mm thick and with a weight of 699 g/m² of weight. The chosen colour of the blinds is yellow gloss. Table 3 summarizes the DSF properties.

The DSF is exposed to the controlled climatic conditions (temperature and solar radiation), provided by the weather generator and the solar simulator. In order to represent precisely a mechanically ventilated façade and ensure accurate airflow control a special air supply system was designed (Figure 5). This system is connected to the lower ventilation opening of the DSF. It is composed of a ventilator which blows the air into a convergent diffuser by the intermediate of a duct (Figure 5a). A differential pressure flowmeter is installed inside this duct to measure accurately the airflow blown by the ventilator into the air channel of the façade (Figure 5c). Preliminary tests show that the velocity at the diffuser outlet is nearly uniform.

Metrology

The metrology installed on the experimental facility is composed in two systems, one for the test cell and the other for the DSF. Both measurement units are composed of 70 thermocouples and 7 RTD sensors, one inline flowmeter and one pyranometer. Table 4 summarises the sensors properties. The entire test cell is equipped with type-K thermocouples on its interior and external surfaces (facing inside and facing thermal buffer zone). 9 thermocouples are located on each face of five walls in order to obtain the mean surfaces temperatures.

The DSF airflow measurement system is composed of two differential pressure devices (two blades, thus limiting the differential heads and turbulences) and a sensor box with the differential pressure elements (Figure 5c). These two pressure elements have a distribution of apertures on the airfoil permitting the control of the average airflow (by averaging the differential pressures). The pressure transmitters working with a differential probe are configured with a square root function. Via this function, and from the differential pressure, the transmitter calculates air velocity and/or airflow in the duct, and thus in the DSF.

The DSF itself is instrumented using thermocouples, distributed on all the elements of the façade (Figure 6). Type-T thermocouples with flattened joint are used here for a better contact with the glass surface. Among these sensors (Table 4), 15 thermocouples are affected for the glass surface temperatures and 10 for air temperatures (Figure 6). All the thermocouples installed in the DSF are shielded to the solar simulator radiation using a special paint. This paint is Labsphere's Spectraflect.

Air velocity transducers are also installed. The sensors are TSI 8475 Air velocity Transducer. They sensors are selected to measure air velocities within the mechanically ventilated cavity. Thus, four air velocity transducers were used for the measurements. The transducers were located at the middle height (+1.15m) and in the upper part (+2.1m) of the façade cavities (external and internal air channels). Each air velocity transducer is suspended vertically in the cavity with the sensors facing the upward airflow. The uncertainties concerning velocity measurements are given in Table 4.

The test cell sensors are connected to a multiplexing station composed of 10 acquisition cards, each one having 39 inputs. Measurements are then carried out using a multimeter Keithley 2700. DSF sensors are connected to a dedicated Multimeter/Data Acquisition System. Acquisition and recording of the data are done by a computer station using Labview software.

Experimental procedure and preliminary results

Experimental procedure

Every test consists in setting up parameters of the test cell facility (boundary conditions) and DSF.

Simultaneously, the solar simulator is switched on and synchronized with the acquisition system. The time delay between two measurements is 300 s (i.e. a measurement frequency of 0.003 Hz). The test cell reaches the steady state condition after 48 hours when the DSF is ventilated, respectively 72 hours with no airflow in the DSF. All the values presented in this paper were obtained averaging 100 values recorded at steady state. Comprehensive and more detailed information about the configurations are given in Figure 7.

Thermal buffer zone

Given that the thermal buffer zone temperature around all the five boundaries of the test cell is accurately controlled, it is possible to test several cases. Among these, a main scenario was selected for the tests. This scenario is represented by a summer situation which is representative of the hot period. This configuration of the thermal buffer zone simulates the cell test in contact with adjacent rooms at constant temperature. Therefore, the temperature of the thermal buffer zone is fixed according to the studied case.

Thus, the air temperature of the thermal buffer zone was fixed at 26°C (Figure 7).

Weather generator

During tested scenario, the air temperature inside the weather generator is fixed at 32°C (Figure 7). The weather generator oscillation around the temperature setpoint is reduced to approximately ±0.2-0.3°C for all configurations, which confirms the stability of the weather generator regulation system.

Solar simulator radiation results

During all the tests, the solar simulator is unchanged constant. The orientation angle of these lamps is chosen following preliminary tests, to uniformly irradiate the DSF structure.

To obtain uniform radiation repartition, the lamps are redirected until reaching a homogenous radiative solar flux on the external DSF glass pane. The measurements are realised simultaneously on a virtual grid (Figure 8) on the external and the internal glass pane. With these measurements, is possible to evaluate the values of solar radiation in different location of the DSF. Figure 9 presents the solar irradiation isovalues on the external glazing after the reorientation of the lamps. External glass sheet receives the most significant part of the radiation of solar simulator spots. The reorientation of lamps allows to improve and to normalize the initial regions with a high irradiation (in real case, the sun is sufficiently far to act as punctual source and to radiate uniformly). Figure 10 presents the same direction as Figure 9, but the measures were taken on the internal DSF glass pane. Figure 10 shows as well, that a significant part of the solar radiation is absorbed or reflected by the glass panes or by the venetian blinds of the DSF. Figure 11 summarises all the internal solar radiation measurements for blinds angles (β) varying from 0 to 90 degrees. For convenience, the measured values on different location of the façade glass (exterior and interior, Figure 9 For each change of the blinds angle the average irradiation is decreasing by approximately 50% comparing to previous angle (for example, between horizontal blinds and 30° blinds). Moreover, for completely closed blinds (90°), since the theoretical transmittance is equal to 0, the first consequence can be that no radiation will enter inside test cell. But with multiple reflections and solar radiation leakage regions of the façade, the test cell receives solar radiation.

DSF airflow and air velocities results

In this part the airflow and the air velocities inside the DSF are analysed. The air velocities are analysed by the intermediate of a coefficient. This coefficient (see equation 2) is a function of measured air velocities inside the DSF gaps, at middle height. As well, this coefficient depends of the ventilation opening surface and of the imposed airflow at the entering of the DSF.

v vent H Q S v     15 . 1
(2) Figure 12 3 Results

In this part, results obtained from the measurements are detailed. Various scenarios are considered.

Temperature results

The scenario is characterized mainly by the weather generator air temperature fixed at 32°C (a hot summer day) and the buffer zone temperature fixed at 26°C (the test cell is supposed to be in contact with a air-conditioned zone). Under these aspects, several configurations are scrutinised: different airflows combined with different orientation of the shading device. In our case, the presence of venetian blinds allows to vary the angle between 0° (completely opened blinds) and 90° (completely closed blinds). The airflow range varies between 0 and 600 m 3 /h corresponding to 0-193 m 3 /h for one meter of façade. All the tested configurations are presented in Figure 7.

Table 5 presents the totality of the boundary conditions of the test cell facility; the temperatures of the thermal buffer zone air volumes (North and South face, Ceiling and Floor). The fluctuations between the setpoint and the temperature in the thermal buffer zone, on the North, South and top air volumes are around 1.3°C. For the bottom face of the thermal buffer zone, in contact with the floor, the temperature recorded during the measurements is in the region of 22°C. Table 6 presents the temperatures recorded on interior surface of the walls of the test cell. For convenience, the 9 thermocouples recordings are time averaged and only one value is presented for every wall. Similarly to thermal buffer zone, the walls temperatures are averaged from 100 values recorded at steady state. For example, during the measurements with the completely opened blinds, increasing the airflow in the façade is decreasing the walls surface temperatures (North, South, Ceiling and Floor). For the East wall (the back wall of the test cell) the temperature is relatively constant for all airflows. This fact is due to the solar radiation passing through the façade, irradiating the wall surface. Modifying the blades angles, the solar direct flux through the DSF is reduced and the consequence is observed on the back wall surface temperature. This wall surface temperature decreases simultaneously with the increasing of the blinds angle.

Installing thermocouples in different points of the DSF (Figure 6) and inside test cell (6 Pt100 sensors are installed to measure the vertical temperature profile), a temperature evaluation is realised for different location. Consequently, every sensor will record the temperature of the near environment (surface or air).

Thus, for measured configurations (airflow and blades angle) top and bottom recordings of test cell and DSF layers were averaged at steady state. Table 7 shows these temperatures for all the summer configurations. In the light of this situation is possible to calculate a temperature difference top-bottom for each element. This vertical temperature difference of the DSF elements, particularly for air temperatures, is important to be evaluated. This evaluation shows if the DSF can increase sufficiently the air temperature, air that can be provided for ventilation. This situation is supposed advantageous during middle or moderate cold season.

Figure 13 shows that blinds angle are highly influencing the temperature regimes DSF and test cell. The test cell air temperature is more influenced by the blinds angles than by the airflow. Thus, closing progressively the blades, a more consequent decrease of temperature is produced inside the test cell than by doubling the airflow.

For both mechanical ventilation and zero ventilation through the DSF, blinds temperature is superior to all other temperatures (air and surfaces of the DSF). This temperature difference vary proportionally with the blinds angles, from approximately 1.5°C for blinds at 0°, respectively to 3°C for blinds at 90°(closed sun-shading system). The blinds increased temperature is due to the progressive interception of the solar radiation. Thus, more the blinds are closed more their temperature is increasing.

Figure 14 shows an outlook of the entire set of temperature measurements for the test cell. Measured temperature inside test cell varies few for angles superior to 60. Subsequently to our measurements and strictly thermo-energetic, hot season strategies with angles superior to 60° (even with important airflow in the DSF) cannot be considerate giving a thermally sufficient impact. Thus, since a 60° angle gives a similar decrease of temperature as a 90° angle (closed blinds), a lower angle is more suitable due to his higher daylighting potential.

DSF energy efficiency results

This consequence presented in Figure 14 is reinforced calculating the energy efficiency of the DSF (see equation 3). This coefficient is a measure of the energy "removed" from the test cell by the airflow imposed in the DSF. This coefficient depends of the measured air temperatures at the lower and upper ventilation openings of the DSF. Moreover, this coefficient is correlated with the solar simulator radiative flux and the glass surface of the DSF.

g ext S in out p m S E T T c Q E       ) ( (3)
Therefore, Figure 15 presents the energy efficiency function of different airflow. Increasing more the blinds angle is not traduced by removing more heat from the façade. For a small airflow (200 m 3 /h) the energy efficiency is nearly constant for all angles excepting complete opened blinds. The maximum of energy efficiency is for blinds angle of 45°. After this value, the trend is directed downward which makes that a 60° angle to be equivalent thermo-energetically with a 30° angle. From a natural lighting perspective, small angles can lead to better use of daylight inside allowing more solar radiation to pass towards the room. This limited blinds inclination due to thermal aspects can simplify further the DSF control (high glazed surface with bad consequences on thermal comfort, but with good effects for visual comfort) and thus giving more access to natural light.

Conclusions

In order to investigate glazed building envelope, an experimental research has been carried out with double-skin façades. This study allows measurements of double-skin façades, under controlled thermal and radiative environment. The tests concern the behaviour of the test cell provided with a double-skin façade during the hot season characterised by different set-ups (blinds angle and ventilation).

These experimental campaigns enable us to evaluate precisely the thermal behaviour of a room provided with a double-skin façade, comparing to real building monitoring. However, this entire set of experimental results enables to carry out validation of DSF models. Thus, this extensive data set (boundary conditions, walls surface temperature, DSF temperatures, etc.) can provide a practical resource for further research and for use by building designers and modellers.

In order to validate the use of double-skin façades, further investigations are needed. Various external conditions must be tested to evaluate the gain during spring, autumn. As well, numerical simulations must be carried out, in order to investigate the behaviour for real buildings. The thermal and visual aspects are to take into account simultaneously.

  and Figure 10) are averaged. Thus, Equation (1) defines the shading coefficient ζ, function of the radiative flux on the DSF external and the internal glass sheet. This shading coefficient describes the amount of solar radiation which is stopped by closing progressively the blinds.

S

  Specific heat capacity [J/kg K] out T Air temperature in the upper part of the DSF [°C] in T Air temperature in the lower part of the DSF [°C] g Glass surface of the DSF [m²]

Figure 2 .

 2 Figure 2. Comparison of a CSI simulator spectrum with solar spectrum Figure 3. Large scale solar simulator. (Metric units) Figure 4. DSF section view: 1, 6 mm glass layer; 2, motorized solar protection; 3, metallic support; 4, solar protection action box; 5, concrete beam; 6, aluminium frame; 7, pane opening articulation; 8, ventilation openings. (Metric units) Figure 5. Double-skin façade air supply system. (Metric units) a) Section view. b) Plan view. c) Flowmeter Figure 6. Thermocouples position on double-skin façade elements. (Metric units) Figure 7. Measurements configurations Figure 8. Solar insulation measure grid (external and internal glass). (Metric units) Figure 9. Solar simulator insolation distribution on the external glass pane Figure 10. Solar simulator insolation distribution on the internal glass pane (β=0°) Figure 11. Shading coefficient of the DSF function of different blinds angle Figure 12. Middle height air velocities inside the façade cavities a) External cavity; b) Internal cavity; Figure 13. Double-skin façade temperature profile. a) Blinds at 0° angle (completely opened); b) Blinds at 30° angle; c) Blinds at 45° angle; d) Blinds at 60° angle; e) Blinds at 90° angle (completely closed); Figure 14. Test cell air temperature Figure 15. Energy removed from the test cell

Table 1 .

 1 Wall composition (from the interior of the test cell to the thermal buffer zone)

	Wall	Material	Thickness [mm]
	Floor	Cellular concrete	200
		Plasterboard	10
	Vertical wall	Polystyrene Plasterboard	50 10
		Agglomerated wood	50
		Plasterboard	10
	Ceiling	Plywood Mineral wool	8 55
		Wood	25

Table 2 .

 2 Thermophysical properties of the envelope Figure 1. Test cell facility scheme: 1, cooling unit; 2, weather generator; 3, double-skin façade; 4, protection glass of the solar simulator; 5, test cell; 6, concrete; 7, air blowing plenum; 8, solar simulator's heat removal ventilators; 9, air extraction plenum; 10, HVAC unit of the buffer zone; 11, solar simulator; 12, controlled buffer zone. (Metric units)

	Material	Conductivity [W/m 2 °C]	Capacity [J/kg°C]	Density [kg/m 3 ]
	Plasterboard	0.35	1620	817
	Agglomerated wood	0.136	1640	544
	Plywood	0.11	1600	417
	Polystyrene	0.04	1380	25
	Mineral wool	0.06	1000	72
	Siporex concrete	0.16	1000	368