Julien Bruneau
email: julien.bruneau@inria.fr

Charles Consel
email: charles.consel@inria.fr

Marcia O'malley
email: omalleym@rice.edu

Walid Taha
email: walid.taha@hh.se

Wail Masry Hannourah
email: wailhannourah@yahoo.com

Virtual Testing for Smart Buildings

Keywords: Virtual testing, Smart buildings, HVAC, Pervasive computing

Smart buildings promise to revolutionize the way we live. Applications ranging from climate control to fire management can have significant impact on the quality and cost of these services. However, smart buildings and any technology with direct effect on human safety and life must undergo extensive testing. Virtual testing by means of computer simulation can significantly reduce the cost of testing and, as a result, accelerate the development of novel applications. Unfortunately, building physically-accurate simulation codes can be labor intensive.

To address this problem, we propose a framework for rapid, physically-accurate virtual testing. The proposed framework supports analytical modeling of both a discrete distributed system as well as the physical environment that hosts it. The discrete models supported are accurate enough to allow the automatic generation of a dedicated programming framework that will help the developer in the implementation of these systems. The physical environment models supported are equational specifications that are accurate enough to produce running simulation codes. Combined, these two frameworks enable simulating both active systems and physical environments. These simulations can be used to monitor the behavior and gather statistics about the performance of an application in the context of precise virtual experiments. To illustrate the approach, we present models of Heating, Ventilating and Air-Conditioning (HVAC) systems. Using these models, we construct virtual experiments that illustrate how the approach can be used to optimize energy and cost of climate control for a building.

I. INTRODUCTION

Buildings are designed to achieve a wide range of goals. For example, ensuring occupant safety and health requires specialized appliances such as fire alarm systems [1] and proper ventilation [START_REF]Refrigerating and Air-Conditioning Engineers website[END_REF]. These concerns are very real: Each year, hundreds of people die from carbon monoxide poisoning due to poor ventilation design. Similarly, ensuring occupant comfort (which is generally essential for ensuring their productivity [START_REF] Hedge | Linking environmental conditions to productivity[END_REF]) also requires specialized appliances. The stringent uptime requirements placed on building appliances also give rise to concerns about energy [START_REF]Energy Star website[END_REF], cost, and environmental impact [5]. Traditional technologies such as brick and mortar, standard heating and airconditioning units, have all been carefully scrutinized with respect to their impact on all of these goals, and have passed the test of extended timelines. Often, the design of such types of systems is intensionally kept simple and passive in order to meet some of these goals, such as safety and maintainability.

1) Smart Buildings.: Recently, there has been an increasing interest in smart building technologies. Such technologies bear a promise to revolutionize the way we live [START_REF] Harper | Inside the Smart Home[END_REF]. Applications ranging from climate control to fire alarm systems, to lighting management, to security systems can be found in smart buildings [START_REF] Lin | Internet-based monitoring and controls for hvac applications[END_REF], [START_REF] Mueller | Demo: A generic platform for sensor network applications[END_REF], [START_REF] Lee | Network-based fire-detection system via controller area network for smart home automation[END_REF].

Intrinsic to the idea of smart buildings is the introduction of higher levels of active control into the traditional components of a building. Higher levels of active control are achieved by the use of more sophisticated control algorithms, more extensive sensing of the physical environment, more actuation of various physical subsystems, and more communication between different components of the system. However, because of the direct impact that buildings have on humans and on their resources, smart building technologies must be held to stringent standards for correctness, availability, and a wide range of other software quality considerations.

2) Virtual Testing.: Testing is a crucial tool for eliminating poor designs, and developing a degree of confidence in promising designs. Testing is equally important for traditional physical design and design involving active control. But testing smart building technologies in physical buildings can be slow and prohibitively expensive. Addressing these two problems can accelerate the rate of innovation and deployment of successful smart building applications. In particular, using computer simulations to carry out some part of the testing virtually can help achieve this goal. The effectiveness of this approach, however, depends heavily on the accuracy with which we model both the control and the physical components of the system. Furthermore, building accurate simulation codes, especially for physical systems, can be labor intensive and can slow down the whole development process for smart building applications.

The goal of this paper is to address three technical challenges that must be overcome in order to enable effective virtual testing of smart buildings. The first is to accurately capture the distributed and networked nature of the active devices in the system. The second is to accurately capture the physical properties of the building. The third is to auto-matically map such models directly to executable simulation codes.

3) State of the Art.: Existing approaches only cope at most with one of the three challenges raised by the virtual testing of smart buildings. For instance, several projects simulated building active devices using MAT-LAB/Simulink [START_REF] Riederer | MATLAB/Simulink for Building and HVAC Simulation -State of the Art[END_REF]. These projects focus on the fine-grained modeling of these devices. They provide libraries of digital components that can be used to model devices. However, they do not attempt to use analytically sound models of the physical environment surrounding such devices. COMSOL allows to accurately simulate the surrounding physical environment. For instance, it provides a heat transfer module and an acoustics one. However, these simulations are based on the Finite Element Method (FEM) and are prohibitively expensive for modeling the physical environment of a whole smart building. Other tools allow faster simulation of the physical properties of a building. Modelica is one of these tools. Modelica is an equation-oriented modeling language. The main draw back is that modeling systems in Modelica that combine discrete and continuous behaviors can be somewhat challenging.

4) Contributions.:

In this paper, we present a new framework that allows an effective virtual testing of smart buildings where

• Physical aspects of smart homes are modeled explicitly in Acumen [START_REF] Zhu | Mathematical Equations as Executable Models of Mechnical Systems[END_REF], a domain-specific language with specialized support for describing continuous systems (Section II). • The active components of smart homes are modeled explicitly in DiaSpec [START_REF] Cassou | A generative programming approach to developing pervasive computing systems[END_REF], a domain-specific language with specialized support for modeling distributed, digital systems (Section III). • The complete models, containing both Acumen and DiaSpec components, are mapped to executable codes (Section IV). This is achieved by combining Acumen's simulation capability and the DiaSim simulator [START_REF] Bruneau | Diasim, a parameterized simulator for pervasive computing applications[END_REF].

• A virtual testing platform to monitor and gather statistics about smart building appliances can be built (Section V). • Virtual experiments using different building models and control strategies can be analyzed (Section VI).

II. MODELING THE PHYSICAL ENVIRONMENT

Virtual testing of smart building designs requires accurately capturing its physical properties. A building is a multiphysics system involving multiple interrelated physical characteristics. Physicists have already defined these phenomena with mathematical definitions (e.g., with ordinary/partial differential equations). These analytical descriptions are an ideal material to reuse in order to capture the physical properties of a building. In this section, we first explain the approach that we adopt to modeling heat transfer and temperature change in a building. We then show how the differential equations that capture this model are expressed in Acumen. For reasons of space, we only present our temperature model. However, we have also modeled relative humidity and carbon dioxide density in Acumen.

A. Modeling Temperature and Heat Transfer in a Building

Human comfort and safety are highly sensitive to the temperature of the surrounding air. As a result, it is critically important to accurately model the factors that impact temperature, including appliances that can be used to control it, as well as the processes by which the temperature of the air in a given room is changed.

We present the models used in our study as a series of successive refinements of a basic model. All models are compartment models, in that they treat each room as one state variable. The basic model, as well as the refinements which include additional terms, are all differential equations.

1) Heat Transfer.: Heat transfer is the rate at which heat moves through a medium or from one medium to another, and is a topic studied extensively in thermodynamics (e.g., [START_REF] Kuehn | Fundamentals: 2005 Ashrae Handbook[END_REF]). Heat transfer between two media is linearly proportional to the difference in temperature between the two media. In addition, it is also affected by the thermal resistance of the boundary between the two media. For simplicity, we assign each room in a building one temperature value. Reasonable values for the thermal resistance of building walls, windows and doors can be determined using a reference book in the heating and cooling domain [START_REF] Kuehn | Fundamentals: 2005 Ashrae Handbook[END_REF].

Let us assume that all rooms are numbered. We will use the subscripts i and k to refer to room numbers. Let N eighbors(i) be the set of numbers representing the rooms neighboring room i. Let T i denote the temperature of room i.

To help introduce the reader to the notation and the equations that define heat transfer in a building, we begin by assuming that the only factor affecting temperature in a particular room is heat from neighboring rooms. To express even this simple process, we need some additional notation. In particular, we will also use the following convention:

dTi dt Rate of temperature change in room i ( • C.h -1 ), C i Thermal capacitance of room i (J. • C -1 ),
R ik Thermal resistance of the boundary between rooms i and k ( • C.h.J -1 ). It takes into account the heterogenous elements of this boundary (e.g., walls, windows, doors). The equation constraining the rate of change for each and every room i is given by the equation:

dT i dt = 1 C i k∈N eighbors(i) T k -T i R ik (1)
Because this equation is instantiated for each room, the whole building is modeled by a set of such equations.

2) Air-Conditioning Unit.: We now consider adding airconditioning (AC) units to our model. An AC unit consists of a heater and a cooler. We need only to introduce four additional parameters:

B h(i) The heater in room i is active (1 if present and active, 0 o.w.), P i

Heater power (W), B c(i) The cooler in room i is active

(1 if present and active, 0 o.w.), Q i

Cooler power (W), The equation above only needs to be extended as follows:

dT i dt = 1 C i ( k∈N eighbors(i) T k -T i R ik +(B h(i) * P i -B c(i) * Q i ))
(2) 3) Occupants.: Occupants can be modeled as heat sources. The set of occupants of room i is denoted Occupants(i). We only need one additional parameter to incorporate this aspect of building:

H j Heat dissipation of occupant number j (W), Thus, the final equation can be expressed as:

dT i dt = 1 C i k∈N eighbors(i) T k -T i R ik + 1 C i * (B h(i) * P i -B c(i) * Q i ) (3)
+ 1 C i * j∈Occupants(i) H j
Other heat sources, such as equipment, appliances, and lights, were neglected for simplicity but will be included at a later stage. Now, we are ready to present the actual Acumen code used in the experiments we report on in the rest of the paper.

B. The Heat Model in Acumen

By design, the Acumen modeling language enables direct specification and simulation of continuous and discrete systems. Our virtual testing framework only uses its continuous system modeling and simulation capability, and not its support for modeling discrete behaviors.

Figure 1 presents the Acumen specification of the temperature defined in Equation 3. The continuous section in Figure 1 specifies the temperature rate of change for each room of the building. Equations in Acumen can refer to derivatives of variables. For example, T' refers to the first derivative of T. Finally, the boundary conditions subsection allows one to define the initial state of the physical environment. This initial state can be easily changed by setting new boundary conditions. As can be noticed, defining physical characteristics in Acumen is straightforward and has a direct correspondance to equational definitions. Thus, Acumen leverages standard mathematical notation used to define physical phenomena.

( * Building topology, Room 0 corresponds to the outside * ) building=((0),(0,2),(0,1,3),(0,2));

( * Temperature in each building room, T0 is the outside temperature * ) T=(T0,T1,T2,T3);

( * Other variable definitions * ) ... continuous foreach room in length(building) begin T'[room] = 1/C[room] * ((sum n < length(building[room]) in ((T[building[room][n]]-T[room]) / R th[room][n])) + Bh[room] * P[room]-Bc[room] * Q[room] + (sum p < length(occupants[room]) in H[occupants[room][p]]

III. MODELING ACTIVE COMPONENTS

Active systems deployed in smart buildings involve a wide range of devices and software components, communicate using a variety of protocols, and rely on intricate distributed systems technologies. Accurate virtual testing for smart buildings requires accurately capturing the distributed and networked nature of these active systems.

We rely on the DiaSpec language for modeling active systems. DiaSpec first consists of a taxonomy language dedicated to describing classes of devices that are relevant to the target active systems. These devices can be hardware or software. DiaSpec also provides an Architecture Description Language (ADL) layer. This ADL layer enables to model active controllers.

In this section, we explain how a taxonomy of devices is created in DiaSpec, present the ADL layer, and briefly summarize the benefits of using DiaSpec for modeling active systems in smart buildings.

A. Modeling Active Devices in DiaSpec

DiaSpec provides a dedicated language for specifying the taxonomy of devices relevant to a target application domain. As shown in Figure 2, a device consists of sensing capabilities, producing data, and actuating capabilities, providing actions. Following this structure, a device description declares a data source for each one of its sensing capabilities. For example, the TemperatureSensor class of devices in Figure 2 provides a temperature data source. An actuating capability corresponds to a set of method declarations. Thus, the Heater class of devices provides an OnOff action. This action allows to turn on and off a heater instance. A device declaration also includes attributes, characterizing properties of device instances (e.g., the LocatedDevice instances are characterized by the room in which they are deployed). Finally, device declarations are organized hierarchically allowing device classes to inherit attributes, sources and actions. In our example, the Heater and TemperatureSensor classes of devices inherit the room attribute.

B. Modeling Active Controllers in DiaSpec

DiaSpec provides an ADL layer for specifying the architecture of active controllers. This layer is dedicated to an architectural pattern commonly used in the pervasive computing domain [START_REF] Dey | A conceptual framework and a toolkit for supporting the rapid prototyping of context-aware applications[END_REF]. Active controllers are decomposed in two types of components: context and controller. Context components are fueled by sensing devices. These components filter, interpret and aggregate these data to make them amenable to the applications needs. Controller components receive application-level data from context components and determine the actions to be triggered on devices.

Figure 3 illustrates the definition of an active controller in DiaSpec with a temperature regulator. It is composed of a TemperatureRegulation context component and of a HeaterController controller component. First, we define a TemperatureRegulation context component. This context component receives low-level temperature data from temperature sensors. This context component is responsible for determining whether a temperature regulation is needed. The produced Regulation information is coupled to a room index so that the regulation is realized room per room. Finally, we define a HeaterController component. This controller component is responsible for applying the temperature regulation. It receives regulation information from the TemperatureRegulation context component. It acts on heaters to apply this regulation.

C. DiaSpec and Virtual Testing

Using DiaSpec brings multiple benefits to our virtual testing framework regarding development and simulation support, going beyond a contemplative approach. Not only is a DiaSpec description a high-level model of an active system, it is also processed by a compiler that generates a dedicated Java programming framework. This programming framework enables the developer to easily implement the active components that he modeled using DiaSpec. Indeed, the provided support is dedicated to what was described with DiaSpec. For instance, a remote procedure call from an active controller to a device is realized by simply calling a Java method provided in the generated programming framework. This framework abstracts over the distributed and heterogeneous nature of this type of system. Further information on the generated programming support can be found elsewhere [START_REF] Cassou | A generative programming approach to developing pervasive computing systems[END_REF].

Using DiaSpec for modeling the active system allows us to test active controllers, without code modification, in a 2D simulator, named DiaSim [START_REF] Bruneau | Diasim, a parameterized simulator for pervasive computing applications[END_REF]. The simulation support provided by DiaSim is described next.

IV. MAPPING THE MODELS INTO EXECUTABLE SIMULATION CODES

The final technical challenge to enabling effective virtual testing of smart building is to automatically map the models of both physical environment and active components to executable simulation codes. This section explains how this mapping is done, and explains how the simulations of the Acumen and DiaSpec models interact.

1) Execution of Acumen Models.: Acumen models are directly mapped to executable code. The process of mapping Acumen models to executable code is described in more details elsewhere [START_REF] Zhu | Mathematical Equations as Executable Models of Mechnical Systems[END_REF]. The user only needs to specify the length of the simulation and the step size of each simulation iteration. Continuous variables in an Acumen model are discretized with respect to this step size. The simulation accuracy depends on the chosen step size. However, having a small step size requires larger computation resources. The user needs to properly set the step size value so that an acceptable level of accuracy is obtained in the simulation.

2) Execution of Active Controllers.: From the DiaSpec taxonomy and architecture, the DiaSpec compiler generates a dedicated Java programming framework. The developer uses this framework to implement both the active devices and the active controllers. Once the active controllers are developed, we can execute and test them with the DiaSim simulator. DiaSim allows active controllers to interact with simulated devices. The active controller implementation does not have to be modified in order to be simulated in DiaSim. This allows to test and assess the application logic of the active controllers. The simulated device localizations in the physical environment are graphically edited in the DiaSim editor. Once the physical environment is edited, we generate a simulation support that is used for executing the active system simulation. DiaSim also allows to test active controllers in hybrid environments, mixing virtual and real active devices. This feature allows an incremental deployment of a smart building. DiaSim is more extensively presented elsewhere [START_REF] Bruneau | Diasim, a parameterized simulator for pervasive computing applications[END_REF].

3) Smart Building Simulation.: The Acumen and DiaSim simulations interact in two ways. First, simulated sensors need to retrieve their sensed information from the Acumen model. Second, simulated actuators may modify the state of variables in the Acumen model. These two interactions are achieved using Acumen's socket-based Java API to interact with an Acumen simulation. 1 While editing the simulated active system in the DiaSim editor, the user specifies the Acumen variables of interest for each simulated sensor. Thus, the data sensed by a simulated sensor correspond to a variable in Acumen. Continuous variables in Acumen are discretized with respect to a given step size defined by the user so that the discrete components can use them. Finally, the user specifies which variable is modified by an actuator action and how it is modified.

Finally, we generate Java implementations for executing simulated devices. The necessary glue code to interface with Acumen is generated in the simulated device implementation with respect to what the user has specified during the edition.

V. MONITORING VIRTUAL SMART BUILDINGS

Virtual testing of smart building active components involves modeling and simulation of numerous interacting components. In addition, it is difficult for the user to determine if an active controller behaves correctly without observing the physical environment. To cope with these difficulties, we first provide support for monitoring the virtual testing of smart buildings. We also provide tools for analyzing the results after the completion of virtual testing.

4) 2D Graphical Rendering.: Simulations can be monitored using a 2D renderer provided by DiaSim. This renderer is based on the Siafu simulator [START_REF] Martin | A generic large scale simulator for ubiquitous computing[END_REF]. The image rendered for an example virtual house is presented in Figure 4. This rendering facilitates monitoring the state of the active devices. Information is displayed when a sensor publishes an event and when an action is called on an actuator. Moreover, graphical rendering helps the user visualize the state of the physical environment. Indeed, it is possible to display a greyscale overlay that represents the values of the physical environment characteristics (e.g., temperature and luminosity).

5) Virtual Testing Analysis.: During a simulation, we log the values of the continuous and discrete variables. This allows to save the physical environmnent state and the active device states at each simulation iteration. At the end of a simulation, all logged variables are automatically plotted. This plot can serve as a basis for analyzing events that occured during the simulation. If the user is interested in a more specific analysis (e.g., energy consumption of the 1 Thanks to Cherif Salama for providing this API. active devices), it is easy to create dedicated plots with tools such as GNUplot. For example, we created GNUplot scripts for reading the logged variables and displaying the energy use and comfort level achieved in Section VI.

VI. A VIRTUAL EXPERIMENT

This section presents the active system analyses that enable our virtual testing framework. Active systems can first be evaluated using algorithmic variations. Indeed, multiple algorithms for the same active controller can be compared. Virtually testing these algorithms can help us find the most efficient one with respect to the functionalities of the active controllers. A second variation allowed by virtual testing is structural variation. Multiple configurations of sensors and actuators can be tried to choose the most efficient one. These variations are facilitated by the DiaSpec approach, confining changes to the application logic.

To illustrate both types of variation, we evaluate a Heating, Ventilating and Air-Conditioning (HVAC) system deployed in a house. The plan of the house is illustrated in Figure 4. We apply algorithmic and structural variations to this system for illustrating the usefulness of our virtual testing framework.

A. Heating, Ventilating and Air-Conditioning Systems

HVAC systems are responsible for providing thermal comfort and acceptable indoor air quality to the building occupants. These systems are a standard part of mechanical engineering curricula, see for example [START_REF] Bobenhausen | Simplified design of HVAC systems[END_REF], [START_REF] Levenhagen | HVAC control system design diagrams[END_REF]. HVAC systems regulate multiple physical properties of a building. For example, temperature and humidity must be regulated so that the occupants feel comfortable. As well, carbon dioxide density needs to be regulated for keeping an acceptable indoor air quality. Finally, the amount of airflow introduced to an air-conditioned zone must be controlled to remain pleasant for the occupants. To regulate these multiple physical characteristics, HVAC systems interact with numerous active devices. It retrieves information from sensors such as temperature, humidity and carbon dioxide sensors. It also controls heaters, humidifiers and ventilators to provide comfort and acceptable indoor air quality to occupants. Developing an HVAC system that regulates these multiple physical characteristics is complex. Virtual testing is useful for testing such system.

B. Global vs. Local Temperature Management

Multiple Regulating and Standards organizations define comfortable ranges of temperature depending on the indoor humidity and the outside temperature. For instance, ASHRAE defines the comfortable humidity and temperature ranges in its "Thermal Environmental Conditions for Human Occupancy" standard. We use these values for evaluating our HVAC control algorithms.

Originally, to keep the temperature in the comfort zone, HVAC systems had a single thermostat and regulated the temperature everywhere using this single thermostat. We call this strategy global temperature management. Global temperature management is obviously not optimal, because the temperature can be different in each room. To achieve finer grained, and more efficient regulation, EnergyStar recommends managing the temperature locally [START_REF]A Guide to Energy-Efficient Heating and Cooling by Energy Star[END_REF]. In this case, the building is divided in areas, each area with its own thermostat. We call this strategy local temperature management. To illustrate virtual testing, we set up an experiment to evaluate the effect of applying this recommendation to the HVAC system of a 3-room house. We then observe the comfort and the energy use differences between these two temperature management policies.

C. HVAC System Evaluation

Both global and local regulations use the same regulation logic. The only difference is the scope of this regulation. We name T c-min and T c-max , respectively the lowest comfortable temperature and the highest comfortable temperature. Our regulator needs to keep the temperature in a range [T min , T max ]. If the temperature is below T min , the HVAC system ventilates hot air. If the temperature is above T max , the hot air stops being ventilated. Obviously, the logic of our regulator is very simple and could be refined. However, it is enough for illustrating the variations we apply for testing this regulator.

In this virtual experiment, we use the temperature model presented in Section II. In this model, our virtually tested HVAC system has a heating power of 1500 Watts and a cooling power of 600 Watts. The boundaries between areas of our virtual physical environment are illustrated in Figure 4. The calculated thermal resistance coefficients depend on the windows, doors and walls that compose these boundaries.

1) Algorithmic Variation.: In this section, we test several algorithms for our temperature regulator and choose the most efficient one with respect to its energy use. T min and T max are defined in Equation 4. We test our HVAC system with different values of T min and T max over a period of one month. We decrease these two values as long as comfort is provided 100% of the time. We choose January because heating is critical during this month. The chosen outside temperatures correspond to Bordeaux average temperatures in January. We test three different algorithms. The percentage of time when the hot air is ventilated is presented in Figure 5. We see that these algorithmic variations bring differences in terms of heater time of use. A low range of temperature results in less heating time. Since the heaters have the same power, the differences in the time of use of these heaters allow to evaluate the energy consumption gain. The third algorithm in Figure 5 allows a 38% gain of energy use compared to the first one.

2) Structural Variation.: We also test two different configurations of active components for comparing global and local temperature regulation. The global configuration consists of one temperature sensor in the living-room and one heater in each room. The local configuration consists of one temperature sensor and one heater in each room. We use the most efficient algorithm from the previous section. We test our HVAC system over the month of January. The comfort provided by these two configurations is presented in Figure 6. Local temperature regulation provides perfect comfort to the three rooms of the house. In comparison, global temperature regulation is not as comfortable for the occupants. The bathroom temperature is comfortable for only 86% of the time. The percentage of time when hot air is ventilated for the two types of regulation is presented in Figure 7. We can see that global temperature management ventilates more hot air than local temperature management. Our virtual testing framework allowed us to see that following the EnergyStar recommandation enables to get a better comfort with less energy use from the HVAC system.

3) Virtual Testing Accuracy.: The user needs to carefully choose the accuracy of his virtual testing. The accuracy of virtual testing depends on two parameters: physical environment model accuracy and size of the time discretization step. The temperature model we present in this paper considers that each room has a single state. A smaller discretization of the space would allow the simulation of the physical environment to be more precise. However, a smaller space discretization requires more computation.

Likewise, the size of the time discretization step impacts the virtual testing accuracy. In our evaluation, we choose a discretization step of one minute. A one minute accuracy is enough for evaluating an HVAC system over a month. However, this criteria needs to be carefully chosen depending on the tested smart building application.

VII. RELATED WORK

Existing tools provide partial means to cope with the three challenges stated in the introduction of this paper.

Several projects focus on modeling and simulating building active devices using MATLAB/Simulink [START_REF] Riederer | MATLAB/Simulink for Building and HVAC Simulation -State of the Art[END_REF]. These projects allow a fine-grained modeling of these devices. They provide libraries of digital components that the tester can use to model his devices. However, they do not allow an analytically sound simulation of the physical environment that impacts such devices. Ptolemy goes one step further than MATLAB/Simulink. Ptolemy provides a library of computation models that the user can compose for modeling and simulating embedded systems. Its main contribution is to allow the simulation of systems that contain heterogeneous computation models. Using Ptolemy would be complementary with our approach. Ptolemy could be used for defining the computation models of active devices, whereas DiaSpec describes the outside interface of these active devices. However, the support provided by Ptolemy for modeling continuous systems is too restricted for modeling the physical environment of a building.

A few simulators exist in the pervasive computing domain [START_REF] Barton | Ubiwise, a ubiquitous wireless infrastructure simulation environment[END_REF], [START_REF] O'neill | A testbed for evaluating human interaction with ubiquitous computing environments[END_REF], [START_REF] Morla | Evaluating a location-based application: A hybrid test and simulation environment[END_REF], [START_REF] Reynolds | Requirements for an ubiquitous computing simulation and emulation environment[END_REF]. Ubiwise [START_REF] Barton | Ubiwise, a ubiquitous wireless infrastructure simulation environment[END_REF] and Tatus [START_REF] O'neill | A testbed for evaluating human interaction with ubiquitous computing environments[END_REF] are built upon 3D game graphics engines, respectively Quake III Arena and Half-Life. However the game graphics engine becomes a burden when it comes to define new scenarios; users can neither add their own actuators and sensors, nor simulate arbitrary physical characteristics. Lancaster [START_REF] Morla | Evaluating a location-based application: A hybrid test and simulation environment[END_REF] and PiCSE [START_REF] Reynolds | Requirements for an ubiquitous computing simulation and emulation environment[END_REF] allow to add new types of sensors and actuators for testing active systems. However, there is no support for modeling and simulating the physical environment that surrounds active systems.

Other projects focus on modeling the physical environment. COMSOL [START_REF]COMSOL website[END_REF] allows to accurately simulate the surrounding physical environment. For instance, it provides a heat transfer module and an acoustics one. However, these simulations are based on the Finite Element Method (FEM) and are too slow for modeling the physical environment of a whole building. Other tools allow a faster simulation of the physical properties of a building. Modelica [START_REF] Tiller | Introduction to Physical Modeling With Modelica[END_REF] is one of these tools, and is in fact a closely related language to Acumen. The differences between the two are primarily in Acumen's support for binding time separation and partial derivatives [START_REF] Zhu | Mathematical Equations as Executable Models of Mechnical Systems[END_REF], but these are in fact orthogonal to the models used here. For more sophisticated models of heat transfer, however, partial derivatives are needed.

Our study of the related works showed that testing virtual smart building applications with existing tools is a very complex task. Our approach is a step towards decreasing this complexity by using DSLs to model and execute the physical and digital aspects of smart buildings. Indeed, these DSLs provide the user with a declarative and high-level support that simplifies the virtual testing of these applications.

VIII. CONCLUSION

In this paper, we have presented a virtual testing framework for smart buildings. This framework first provides an analytically sound modeling of the physical environment of smart buildings using Acumen. Our framework also enables the modeling of active components of a building using DiaSpec. Finally, our virtual testing framework map directly both the physical environment model and the active component model into executable simulation codes.

This work can be expanded in several directions. First, we plan on testing other smart building applications in our virtual testing framework. One of these applications is a fire alarm system. Another future work is to couple our virtual testing framework for smart buildings to a human behavior model. Such coupling would allow to see how people behave depending on the active components present in a building, and depending on the physical environment. An example would be to observe people behavior in a building in case of fire.

  )); end boundary conditions T0 with T0(0) = 10; ( * We define the boundary conditions for all continuous variables * )...

Figure 1 .

 1 Figure 1. Temperature specification in Acumen.

Figure 2 .

 2 Figure 2. Device taxonomy in DiaSpec.

contextFigure 3 .

 3 Figure 3. Temperature regulator architecture in DiaSpec.

Figure 4 .

 4 Figure 4. 2D graphical rendering of a virtual house.

Figure 5 .

 5 Figure 5. Comparison of different algorithms for regulating the temperature in terms of energy use.

Figure 6 .

 6 Figure 6.Comparison of the comfort provided by global and local temperature management in the house.


Fiona Gillespie

Wayne K Metcalfe

Patricia Dirrenberger

Olivier Herbinet

Pierre- Alexandre Glaude

Frédérique Battin-Leclerc

Henry J Curran

Measurements of flat-flame velocities of diethyl ether in air

Keywords: Laminar burning velocity, Diethyl ether, DEE, High temperature, Combustion

This study presents new adiabatic laminar burning velocities of diethyl ether in air, measured on a flat-flame burner using the heat flux method. The experimental pressure was 1 atm and temperatures of the fresh gas mixture ranged from 298 to 398 K. Flame velocities were recorded at equivalence ratios from 0.55 to 1.60, for which stabilization of the flame was possible. The maximum laminar burning velocity was found at an equivalence ratio of 1.10 or 1.15 at different temperatures. These results are compared with experimental and computational data reported in the literature. The data reported in this study deviate significantly from previous experimental results and are well-predicted by a previously reported chemical kinetic mechanism.

Introduction

Oxygenated biofuels are important in today's world due to a decline in oil and fossil fuels and an increase in CO 2 emissions and atmospheric pollution. Furthermore, as the price of fossil-derived fuels rises, focus is being drawn towards bio-derived fuels, the most common of which is bio-ethanol. Bio-ethanol is used predominantly in the United States and Brazil, in mixtures with gasoline of 10 and 85 vol. %; however, bio-ethanol is produced primarily from food sources (corn, sugarcane, etc.) and its energy output in not ideal [START_REF] Kohse-Hoinghaus | Biofuel combustion chemistry: from ethanol to biodiesel[END_REF] and [START_REF] Pimentel | Ethanol fuels: energy balance, economics, and environmental impacts are negative[END_REF]. Diethyl ether (DEE) is produced by dehydration of ethanol and is considered to be an excellent compression ignition fuel [START_REF] Bailey | Diethyl ether (DEE) as a renewable fuel[END_REF]. It has a cetane number above 125 and a larger heating value than ethanol.

Davoud and Hinshelwood [START_REF] Davoud | The thermal decomposition of diethyl ether[END_REF] studied the thermal decomposition of DEE with focus on the amount of acetaldehyde present at each stage of the reaction. Freeman, Danby and Hinshelwood [START_REF] Freeman | The thermal decomposition of diethyl ether. I. Rate-pressure relations[END_REF], [START_REF] Danby | The thermal decomposition of diethyl ether. II. Analytical survey of the reaction products as a function of reaction conditions[END_REF], [START_REF] Freeman | The thermal decomposition of diethyl ether. III. The action of inhibitors and the mechanism of the reaction[END_REF], [START_REF] Freeman | The thermal decomposition of diethyl ether. IV. Production of cyanides and the mechanism of the nitric-oxide-induced reaction[END_REF] and [START_REF] Freeman | The thermal decomposition of diethyl ether. V. The production of ethanol from diethyl ether and the pyrolysis of ethanol[END_REF] also studied the thermal decomposition of diethyl ether and focused on rate-pressure relations, effect of reaction conditions, nitric acid as an inhibitor, production of cyanides and production of ethanol. Waddington [START_REF] Waddington | The gaseous oxidation of diethyl ether[END_REF] investigated the gaseous oxidation of DEE at low temperatures and Laidler and McKenney [START_REF] Laidler | Kinetics and mechanisms of the pyrolysis of diethyl ether. I. The uninhibited reaction[END_REF] and [START_REF] Laidler | Kinetics and mechanisms of the pyrolysis of diethyl ether. II. The reaction Inhibited by nitric oxide[END_REF] examined its pyrolysis from 833 to 893 K and 15-320 mmHg and the inhibition of pyrolysis by nitric acid from 833 to 913 K and 10-360 mmHg.

Diethyl ether has previously been investigated as an additive with common fuels in use in different engines. Mack et al. [START_REF] Mack | Investigation of HCCI combustion of diethyl ether and ethanol mixtures using carbon 14 tracing and numerical simulations[END_REF] experimentally investigated the addition of DEE to ethanol in a homogenous charge compression ignition (HCCI) engine in order to find a fuel blend that will elongate the combustion event and reduce the autoignition temperature relative to ethanol. However, they did not observe these desired results due to the similar nature of DEE and ethanol. Mosbach et al. [START_REF] Mosbach | Simulating a homogeneous charge compression ignition engine fuelled with a DEE/EtOH blend[END_REF] numerically simulated this data and predicted that the addition of diethyl ether accelerated combustion due to the greater ignitability of DEE. Miller Jothi et al. [START_REF] Miller Jothi | Experimental studies on homogeneous charge CI engine fueled with LPG using DEE as an ignition enhancer[END_REF] explored the use of DEE as a cetane number improver for use with liquefied petroleum gas (LPG) in an HCCI engine. Their findings show a reduction in nitric oxide (NO), particulate matter (PM) and smoke emissions and an increase in hydrocarbon (HC) and carbon monoxide (CO) emissions compared to diesel fuel. Sudheesh and Mallikarjuna [START_REF] Sudheesh | Diethyl ether as an ignition improver for biogas homogenous charge compression ignition (HCCI) operation -an experimental investigation[END_REF] added diethyl ether to biogas in an HCCI engine in order to improve ignition. This addition resulted in higher brake thermal efficiency and lower HC, NO and smoke emissions, however, an increase in CO emissions was observed. Cinar et al. [START_REF] Cinar | Effects of premixed diethyl ether (DEE) on combustion and exhaust emissions in a HCCI-DI diesel engine[END_REF] studied the effects of adding DEE to diesel fuel in an HCCI direct-injection engine and found that nitrogen oxide (NOx) and soot emissions decreased and CO and HC emissions increased simultaneously. Iranmanesh et al. added DEE to diesel [START_REF] Iranmanesh | Potential of diethyl ether as a supplementary fuel to improve combustion and emission characteristics of diesel engines[END_REF] and biodiesel [START_REF] Iranmanesh | Potential of diethyl ether as a blended supplementary oxygenated fuel with biodiesel to improve combustion and emission characteristics of diesel engines[END_REF] in a diesel engine and found that NOx emissions decreased with addition to both fuels. However, they found that HC emissions increased with addition to biodiesel and carbon dioxide (CO2) emissions increased with addition to diesel.

Recently, Yasunaga et al. [START_REF] Yasunaga | A multiple shock tube and chemical kinetic modeling study of diethyl ether pyrolysis and oxidation[END_REF] studied the oxidation and pyrolysis of DEE at high temperatures behind reflected shock waves. Their results showed that the most important reactions in the combustion of

Agnew and Agnew [START_REF] Agnew | Visible emission spectra of two-stage flames of diethyl ether produced in flat-flame burner[END_REF] investigated the cool and two-stage flames of diethyl ether formed in flatflame burners using infrared radiation techniques and observed emissions from excited formaldehyde in the range from 3200 to 5100 Å. Griffiths and Inomata [START_REF] Griffiths | Oscillatory cool flames in the combustion of diethyl ether[END_REF] studied the oscillatory cool flames of diethyl ether in air experimentally and numerically to investigate the low temperature oxidation of diethyl ether, which, they found, was associated with the hydroperoxide CH 3 CH(OOH)OC 2 H 5 .

Burning velocities of DEE have been measured by Gibbs and Calcote [START_REF] Gibbs | Effect of molecular structure on burning velocity[END_REF] using a Bunsen burner flame cone and the apex-cone method. Di et al. [START_REF] Di | Measurement of laminar burning velocities and Markstein lengths for diethyl ether-air mixtures at different initial pressure and temperature[END_REF] and Zhang et al. [START_REF] Zhang | Experimental study on combustion characteristics of N 2 -diluted diethyl ether-air mixtures[END_REF] and [START_REF] Zhang | Flame instability analysis of diethyl ether-air premixed mixtures at elevated pressures[END_REF] used a spherically propagating flame in a cylindrical vessel to obtain laminar burning velocities and study the combustion characteristics of diethyl ether. However no results have been reported for DEE flame velocities using the heat flux method.

This study presents new adiabatic laminar burning velocities of diethyl ether in air, measured on a flat-flame burner using the heat flux method [START_REF] De Goey | Stabilization of adiabatic premixed laminar flames on a flat flame burner[END_REF]. The experimental pressure was 1 atm and initial temperatures of the fresh gas mixture ranged from 298 to 398 K. These results are compared with the experimental data of Gibbs and Calcote [START_REF] Gibbs | Effect of molecular structure on burning velocity[END_REF] and Di et al. [START_REF] Di | Measurement of laminar burning velocities and Markstein lengths for diethyl ether-air mixtures at different initial pressure and temperature[END_REF], and simulations using the mechanism of Yasunaga et al. [START_REF] Yasunaga | A multiple shock tube and chemical kinetic modeling study of diethyl ether pyrolysis and oxidation[END_REF].

Experimental method

The heat flux method was first proposed by de Goey et al. [START_REF] De Goey | Stabilization of adiabatic premixed laminar flames on a flat flame burner[END_REF] in 1993. This method is beneficial because stretching of the flame does not occur and therefore stretch correction or extrapolation of the data is not necessary. This particular apparatus has previously been used to measure laminar flame velocities of components of surrogate mixtures of natural gases [START_REF] Dirrenberger | Measurements of laminar flame velocity for components of natural gas[END_REF].

In the flat-flame adiabatic burner, a brass burner plate, used to stabilize the unstretched flame, is mounted on a plenum mixing chamber, Fig. 1. The burner plate is 2 mm thick and 30 mm in diameter. The plate is perforated with small holes of diameter 0.5 mm, with the pitch between them being 0.7 mm. Eight Type K thermocouples, of diameter 0.5 mm, are soldered to the plate surface at different distances and angles from the centre to the edge of the plate. The plenum chamber is encompassed by a thermostatic oil jacket, the temperature of which is set to the desired initial temperature of the unburned gas mixture. The circumference of the burner plate is heated with thermostatic oil set to 50 K above the temperature of the unburned gas mixture. In practice, if an initial temperature of 298 K is desired, the temperature of the plenum chamber is set to 298 K and that of the burner plate to 348 K. The oil is applied to the burner plate in order to maintain the plate at a constant temperature and to heat up the cooler unburned gas mixture as it flows through the plate, which compensates for the heat loss necessary to stabilize the flame.

Fig. 1. Schematic of flat-flame adiabatic burner [START_REF] Dirrenberger | Measurements of laminar flame velocity for components of natural gas[END_REF].

Flow rates for liquid compounds are measured using a Bronkhorst mini-CORI-FLOW Mass Flow Controller connected to a Controlled Evaporator Mixer (CEM). The CEM is connected to the plenum chamber using a duct which is wrapped with heating cable in order to keep the compound in the gaseous state as it enters the chamber. Temperature controllers indicate the temperature of the experimental setup at two locations: one below the chamber after the Mass Flow Controllers and one inside the chamber below the burner plate.

Flow rates for gaseous compounds are measured using Bronkhorst High-Tech Mass Flow Controllers (MFCs). At most conditions, the MFC used for oxygen has a maximum flow of 3 nL/min (normal litres per minute) while that used for nitrogen has a maximum flow of 15 nL/min. When the adiabatic flame burning velocity of a gas mixture is above a certain limit, specific to each different temperature tested, two MFCs are used to allow a sufficient flow of oxygen through the burner plate. Bronkhorst conversion coefficient factors are used to convert an MFC designed for nitrogen (with a maximum flow of 3 nL/min) for use with a flow of oxygen.

Adiabatic flame burning velocities in a flat-flame burner are measured using the radial temperature distribution on the surface of the burner plate. If the sum of the heat loss and heat gain is larger than zero, the gas velocity is lower than the adiabatic flame burning velocity. Thus, the centre of the burner plate is hotter than the circumference and the flame is stabilized under sub-adiabatic conditions. Conversely, if the sum of the heat loss and heat gain is lower than zero, the flame is stabilized under super-adiabatic conditions and the centre of the plate is cooler than the circumference. When the flow rate of the gas mixture is changed, an appropriate value of the gas velocity is found that will cancel out the net heat flux. Therefore, the heat loss from the flame to the burner plate is compensated by the heat flux from the burner plate to the unburned gas mixture. The radial temperature distribution in the burner plate is uniform and equal to the temperature of the heating jacket. At this point, the flame is adiabatic, the sum of the heat loss and heat gain from the flame is zero and the adiabatic burning velocity is measured.

Uncertainties in determination of the flame burning velocity are due to errors in:

• Mass flow measurements, with a maximum error of 0.5% per MFC, as given by the manufacturer, Bronkhorst, and a maximum global error in laminar flame velocity of 1.5%.

• The thermocouples, with an error of 0.2 cm s -1 in the laminar flame velocity. This error was measured by changing the gas flow rate by small amounts that did not affect the uniform temperature profile of the flame on the burner plate noticeably.

• Flame distortions, for example edge effects, with an error of 0.2 cm s -1 in the laminar flame velocity in the range of equivalence ratios studied. This error was measured similar to the error in the thermocouples.

Errors of around 1% were present when determining the equivalence ratio, φ, due to mass flow effects, as given by Bronkhorst. The temperature of the unburned gas mixture before and after it flowed through the plenum chamber was measured with a thermocouple. An error of 2 K was estimated. Errors in gas and liquid fuel purity were not significant as high purity compounds were used.

Results

Flat-flame velocities of diethyl ether

The apparatus has previously been validated with methane, which is described in detail by Dirrenberger et al. [START_REF] Dirrenberger | Measurements of laminar flame velocity for components of natural gas[END_REF]. Diethyl ether was studied in air at a pressure of 1 atm, initial temperatures between 298 and 398 K and equivalence ratios from 0.55 to 1.60, Fig. 2. Error bars have been added to these data points. The adiabatic burning velocity increases sequentially with an increase in initial temperature. The maximum flame burning velocity is recorded at φ = 1.15 at 398 and 378 K and at φ = 1.10 at 388, 358, 328 and 298 K. It appears that the maximum burning velocities are shifting to the higher equivalence ratio of 1.15 at higher temperatures, but this is not the case for a temperature of 388 K. This may be due to uncertainties in the flame burning velocity. The adiabatic burning velocity measurements obtained are shown in Table 1.

Fig. 2. Adiabatic burning velocities of diethyl ether.

Table 1. Adiabatic burning velocities (in cm s -1 ) of diethyl ether obtained in the flat-flame burner.

Maximum burning velocities at each equivalence ratio are in bold.

Comparison with literature data

In 1959, Gibbs and Calcote [START_REF] Gibbs | Effect of molecular structure on burning velocity[END_REF] measured the burning velocities of 77 different compounds in order to show the effect of molecular structure on these velocities. Burning velocities of diethyl ether were determined at 298 K and atmospheric pressure, Fig. 3. They found the maximum burning velocity to be at φ = 1.0 with a value of 48.0 cm s -1 , which is in contrast with the results obtained in this study which show a maximum burning velocity of 41.81 cm s -1 at φ = 1.10. This is due to the fact that Gibbs and Calcote did not take into account stretch correction of their flames, as most modern studies do.

In addition, Davis and Law [START_REF] Davis | Determination of and fuel structure effects on laminar flame Speeds of C 1 to C 8 hydrocarbons[END_REF] have reported that the flames studied by Gibbs and Calcote [START_REF] Gibbs | Effect of molecular structure on burning velocity[END_REF] were non-adiabatic and non-equidiffusive, which are factors that have a great effect on the laminar flame burning velocity of the fuel. Therefore, their results are consistently different from those obtained with more modern techniques. Di et al. [START_REF] Di | Measurement of laminar burning velocities and Markstein lengths for diethyl ether-air mixtures at different initial pressure and temperature[END_REF] measured laminar burning velocities of diethyl ether at different temperatures, pressures and equivalence ratios using spherically propagating flames in a combustion vessel. The maximum burning velocity at atmospheric pressure and at all temperatures is recorded at φ = 1.10, which is comparable with measurements from this study and results simulated using the mechanism of Yasunaga et al. [START_REF] Yasunaga | A multiple shock tube and chemical kinetic modeling study of diethyl ether pyrolysis and oxidation[END_REF], Fig. 4. However, the results of Di et al. [START_REF] Di | Measurement of laminar burning velocities and Markstein lengths for diethyl ether-air mixtures at different initial pressure and temperature[END_REF] are quite scattered and deviate significantly from the results obtained in this study and the simulated results, especially at high initial temperatures. Under lean conditions, velocities measured by Di et al. [START_REF] Di | Measurement of laminar burning velocities and Markstein lengths for diethyl ether-air mixtures at different initial pressure and temperature[END_REF] are faster than those obtained in this study and the simulated results, whereas under rich conditions, velocities are slower than those in this study and the simulated results. Variations in these data sets may be due to the stretch correction necessary to acquire laminar flame velocities in combustion vessels with spherical flames. In the case of spherically propagating flames, two conditions occur which must be accounted for [START_REF] Saeed | Measurements of the laminar burning velocity for mixtures of methanol and air from a constant-volume vessel using a multizone model[END_REF]: the rise in pressure in the combustion vessel must be determined correctly and the equation used to measure flame velocities must be valid throughout this pressure rise.

Fig. 4. Comparison of diethyl ether data obtained in this study at 328 and 358 K with data obtained by Di et al. [START_REF] Di | Measurement of laminar burning velocities and Markstein lengths for diethyl ether-air mixtures at different initial pressure and temperature[END_REF] at 323 and 363 K. Results simulated using the mechanism of Yasunaga et al. [START_REF] Yasunaga | A multiple shock tube and chemical kinetic modeling study of diethyl ether pyrolysis and oxidation[END_REF]. Symbols: experiment; lines: simulation.

Comparison with model

Adiabatic flame velocity measurements obtained in this study were compared to simulated results obtained using the mechanism of Yasunaga et al. [START_REF] Yasunaga | A multiple shock tube and chemical kinetic modeling study of diethyl ether pyrolysis and oxidation[END_REF], updated with a new C 0 -C 3 sub-mechanism in construction at the Combustion Chemistry Centre in NUI Galway, Fig. 5 and Fig. 6. The C 0 -C 3 sub-mechanism is based on our previous work [START_REF] Conaire | A comprehensive modeling study of hydrogen oxidation[END_REF], [START_REF] Healy | Oxidation of C 1 -C 5 alkane quinternary natural gas mixtures at high pressures[END_REF] and [START_REF] Lowry | Laminar flame speed measurements and modeling of pure alkanes and alkane blends at elevated pressures[END_REF], with the most recent updates described in [START_REF] Kéromnès | 7th US national meeting of the combustion institute[END_REF] and [START_REF] Metcalfe | Proceedings of the European Combustion Meeting[END_REF]. The experimental results are well-predicted by the model and the overall trend in flame velocities is reproduced. Deviations between the simulated and the experimental results are within the limits of experimental errors. All simulations were performed using Chemkin-PRO [START_REF]CHEMKIN-PRO 15101[END_REF]. Sensitivity analyses were performed using the sensitivity function available in Chemkin-PRO [START_REF]CHEMKIN-PRO 15101[END_REF] and show the effect of a reaction on the flow rate, which corresponds to the burning velocity. These analyses were carried out at an equivalence ratio of 1.10 and an initial temperature of 298 K, Fig. 7, and at an equivalence ratio of 1.15 and an initial temperature of 398 K, Fig. 8. Reactions with negative sensitivity coefficients inhibit reactivity, while reactions with positive sensitivity coefficients promote reactivity. The chain branching reaction between hydrogen atom and molecular oxygen to produce an oxygen atom and a hydroxyl radical is omitted for clarity and has a value of +0.406 at 298 K and φ = 1.10 and a value of 0.411 at 398 K and φ = 1.15. Of the 20 most important reactions, those with smaller molecules, such as methyl, ethyl and carbon monoxide, have the largest effect on the flame burning velocities of diethyl ether. It is important to note that only one reaction involving diethyl ether has a significant effect on flame burning velocities, that between diethyl ether and hydrogen atom to form hydrogen molecule and the secondary radical, C 2 H 5 OC 2 H 4 -s. The reaction between hydrogen atom and hydroxyl radical to produce water is the most inhibiting reaction at both conditions because it is in direct competition with the chain branching reaction to produce oxygen atom and hydroxyl radical.

Conclusions

Adiabatic flame burning velocities of diethyl ether in air have been studied on a flat-flame burner at atmospheric pressure and six temperatures between 298 and 398 K. Equivalence ratios for these experiments ranged from 0.55 to 1.60. It was found that the adiabatic burning velocities increased with respect to temperature. The maximum flame velocity was reported to be at φ = 1.10 at 298, 328, 358 and 388 K and at φ = 1.15 at 378 and 398 K. It was suggested that the maximum velocities were shifting to φ = 1.15 at higher temperatures, but this was not the case at 388 K, due to uncertainties. These results were compared with experimental results from the literature [START_REF] Gibbs | Effect of molecular structure on burning velocity[END_REF] and [START_REF] Di | Measurement of laminar burning velocities and Markstein lengths for diethyl ether-air mixtures at different initial pressure and temperature[END_REF] and were found to deviate significantly. The results of Gibbs and Calcote [START_REF] Gibbs | Effect of molecular structure on burning velocity[END_REF] were non-adiabatic and non-equidiffusive and stretching of the flame was not accounted for. The results of Di et al. [START_REF] Di | Measurement of laminar burning velocities and Markstein lengths for diethyl ether-air mixtures at different initial pressure and temperature[END_REF] differed from the findings of this study, possibly due to the fact that stretch correction was necessary to obtain the laminar burning velocities. A recently constructed chemical kinetic mechanism for diethyl ether oxidation and pyrolysis [START_REF] Yasunaga | A multiple shock tube and chemical kinetic modeling study of diethyl ether pyrolysis and oxidation[END_REF] was compared to experimental results and excellent agreement was achieved. It was found that the flame burning velocities reported in this study compared well to the simulated results. Sensitivity analyses were performed and small species, such as methyl, ethyl and carbon monoxide, were reported to have the greatest effect on the laminar burning velocities of diethyl ether.

Fig. 3 .

 3 Fig. 3. Comparison of diethyl ether data obtained in this study with data obtained by Gibbs and Calcote [23] at 298 K.

Fig. 5 .

 5 Fig. 5. Comparison of diethyl ether data obtained in this study with results simulated using the mechanism of Yasunaga et al. [20] at 298, 358 and 388 K. Symbols: experiment; lines: simulation.

Fig. 6 .

 6 Fig.6. Comparison of diethyl ether data obtained in this study with results simulated using the mechanism of Yasunaga et al.[START_REF] Yasunaga | A multiple shock tube and chemical kinetic modeling study of diethyl ether pyrolysis and oxidation[END_REF] at 328, 378 and 398 K. Symbols: experiment; lines: simulation.

Fig. 7 .

 7 Fig. 7. Sensitivity coefficients showing the 20 most important reactions in diethyl ether flame burning velocities at 298 K and φ = 1.10.

Fig. 8 .

 8 Fig. 8. Sensitivity coefficients showing the 20 most important reactions in diethyl ether flame burning velocities at 398 K and φ = 1.15.

  DEE in a shock-tube are a four-centred elimination (reaction 1), C-O bond fission (reaction 2), and hydrogen abstraction to the secondary radical, C 2 H 5 OC 2 H 4 -s (reaction 3).

	C 2 H 5 OC 2 H 5 = C 2 H 5 OH+CH 4	(1)
	C 2 H 5 OC 2 H 5 = C 2 H 5 O+ C 2 H 5	(2)
	C 2 H 5 OC 2 H 5 +H= C 2 H 5 OC 2 H 4 -s+H

Acknowledgements

We would like to acknowledge the support of Science Foundation Ireland under grant number [08/IN1./I2055], the Region Lorraine, the European Commission through the "Clean ICE" Advanced Research Grant of the European Research Council and COST ActionCM0901 towards this work. Support from Judith Würmel in technical writing and Hervé Le Gall in technical assistance during experiments is also greatly appreciated.


Éric Thiébaut

Ferréol Soulez

Multi-wavelength imaging algorithm for optical interferometry

Keywords: optical interferometry, image reconstruction, inverse problem, sparsity

Optical interferometers provide multiple wavelength measurements. In order to fully exploit the spectral and spatial resolution of these instruments, new algorithms for image reconstruction have to be developed. Early attempts to deal with multi-chromatic interferometric data have consisted in recovering a gray image of the object or independent monochromatic images in some spectral bandwidths. The main challenge is now to recover the full 3-D (spatio-spectral) brightness distribution of the astronomical target given all the available data. We describe a new approach to implement multi-wavelength image reconstruction in the case where the observed scene is a collection of point-like sources. We show the gain in image quality (both spatially and spectrally) achieved by globally taking into account all the data instead of dealing with independent spectral slices. This is achieved thanks to a regularization which favors spatially sparsity and spectral grouping of the sources. Since the objective function is not differentiable, we had to develop a specialized optimization algorithm which also takes into account the non-negativity of the brightness distribution.

CONTEXT

The objective of stellar interferometric imaging is to recover an approximation the specific brightness distribution I λ (θ) of the observed astronomical object given measurements which are provide sparse samples of the spatial Fourier transform of I λ (θ). The reconstruction of a monochromatic image from optical interferometry data is a challenging task which has been the subject of fruitful research and resulted in various algorithms (e.g., MiRA, 1 BSMEM, 2, 3 WISARD, 4 etc.). When dealing with multi-spectral data, a first possibility is to process each wavelength independently and reconstruct a monochromatic image for each subset of measurements from a given spectral channel. For instance, this is what have been done by le Bouquin et al. [START_REF] Le Bouquin | Pre-maximum spectro-imaging of the Mira star T Lep with AMBER/VLTI[END_REF] for the multi-spectral images of the Mira star T Lep. Another possibility is to exploit some assumed spectral continuity of I λ (θ) and process the multi-spectral data globally to reconstruct an approximation of the 3-D distribution I λ (θ). Although demonstrated in a different context of integral field spectral spectroscopy, this latter approach has proven to be more powerful [START_REF] Soulez | Restoration of hyperspectral astronomical data from integral field spectrograph[END_REF][START_REF] Bongard | 3-D deconvolution of hyper-spectral astronomical data[END_REF] and this paper aims at exploiting the advantages of a global multi-spectral processing of optical interferometric data.

In order to simplify the problem, we restricted our study to the cases where the complex visibilities are observed and where the observed scene is a collection of point-like sources. In some sense, this latter assumption makes our algorithm a successor of the CLEAN algorithm [START_REF] Högbom | Aperture Synthesis with a Non-Regular Distribution of Interferometer Baselines[END_REF][START_REF] Schwarz | Mathematical-statistical Description of the Iterative Beam Removing Technique (Method CLEAN)[END_REF] and the building-block method [START_REF] Hofmann | Iterative image reconstruction from the bispectrum[END_REF] developed for recovering monochromatic images from radio and optical interferometric data respectively. In addition to processing multi-variate data, we however introduce the explicit minimization of a non-differentiable regularization term so as to favor spatial sparsity of the reconstructed brightness distribution in a way which is known to be more efficient [START_REF] Donoho | For most large underdetermined systems of linear equations, the minimal ell-1 norm near-solution approximates the sparsest near-solution[END_REF][START_REF] Fornasier | Recovery algorithms for vector valued data with joint sparsity constraints[END_REF] than greedy algorithms like CLEAN [START_REF] Marsh | The objective function implicit in the CLEAN algorithm[END_REF] or the building-block method . [START_REF] Hofmann | Iterative image reconstruction from the bispectrum[END_REF] Our paper is organized as follows: we first summarize the inverse approach for image reconstruction and discuss various possibilities to impose spatial sparsity, we then detail our algorithm for minimizing the objective function; finally we present some results on simulated data and discuss the advantages of our approach.

Further author information: thiebaut@obs.univ-lyon1.fr

ALGORITHMS

Principle of image reconstruction

Following inverse approach described by, e.g. Thiébaut, [START_REF] Thiébaut | Image reconstruction with optical interferometers[END_REF] we state image reconstruction as a constrained optimization problem:

x + = arg min

x∈X {f data (x|y) + µ f prior (x)} (1)
where x ∈ R Card(x) are the sought image parameters, Card(x) is the number of parameters, X ⊂ R Card(x) is the set of acceptable parameters, y ∈ R Card(y) are the measurements, f data (x|y) is a likelihood term which enforces agreement of the model with the data, f prior (x) is a regularization term which introduces a priori constraints and µ > 0 is a so-called hyper-parameter used to tune the relative weight of the prior constraints.

Due to the voids in the u -v coverage of the observations, the constraints set by the data alone do not permit to define a unique image. The prior constraints are then required to help choosing a unique solution among all the images that are compatible with the measurements.

Constraining the solution to belong to the feasible set X is a mean to impose strict constraints such as the positivity:

X = {x ∈ R n : x ≥ 0} (2)
where the inequality x ≥ 0 is taken element-wise.

Likelihood and direct model 2.2.1 Monochromatic case

In this preliminary study, we purposely omit some of the difficulties specific to optical interferometry: we do not impose any normalization constraints on the sought distribution and we assume that complex visibilities have been measured (e.g., phase-referenced observations). In this context, assuming Gaussian noise distribution, the likelihood term writes:

f data (x|y) = (H • x -y) • W • (H • x -y) (3)
where H ∈ R Card(y)×Card(x) is the linear model matrix and W ∈ R Card(y)×Card(y) is a statistical weighting matrix, in principle W is the inverse of the covariance matrix of the measurements: W = Cov(y) -1 .

To avoid dealing with complex arithmetic, we simply consider that complex values are just pairs of real. Hence the total number of observed spatial frequencies is ≈ m/2. For monochromatic interferometric data, the model matrix then writes:

H m,n =    + cos(θ n • b m /λ) if y m is the real part of a measured complex visibility -sin(θ n • b m /λ) if y m
is the imaginary part of a measured complex visibility (4) with b m the observing baseline for the mth datum, λ the wavelength and θ n the angular direction of nth picture element (pixel):

x n ≈ I λ (θ n ) .

In Eq. ( 4), θ n • b m is the usual scalar product between the baseline vector b m and the view direction θ n .

At least because of the strict constraints imposed by the feasible set X, that is whatever the type of regularization, solving the image reconstruction problem in Eq. ( 1) must be done by an iterative optimization algorithm. Thus, owing to the size of the problem, it is strongly recommended to use a fast approximation of H. The nonuniform fast Fourier transform [START_REF] Fessler | Nonuniform Fast Fourier Transforms Using Min-Max Interpolation[END_REF] (NU-FFT) is a good choice for that purpose. In this case, the operator H has the following structure:

H = R • F • S ( 5
)
where F is the discrete Fourier transform (DFT), computed by a fast Fourier transform (FFT) algorithm, R interpolates the discrete spatial frequencies resulting form the DFT at the frequencies effectively observed and S is an apodizing operator (usually S is diagonal) which pre-compensates for the convolution by the interpolation kernel used in R.

Multi-Spectral case

Although in principle any scheme is allowed for sampling a multi-variate distribution such as I λ (θ), it is usually more practical to sample the sought distribution separately along its spatial and spectral dimensions. Hence, for multi-spectral image reconstruction we consider that the image parameters are:

x n, ≈ I λ (θ n ) (6)
for λ ∈ W the list of sampled wavelengths and θ n ∈ A the list of angular directions corresponding to the pixels.

For the sake of simplicity, we use the same wavelengths in W as the ones of the observed spectral channels and we use the notation y m, for the multi-spectral data obtained with mth baseline in th spectral. This notation is intended to clarify the equations and does not impose or assume that all baselines have been observed in all spectral channel. Note that since the both real and imaginary parts are measured, each baseline appears at least twice with this convention.

In our framework, the model of the observed complex visibilities writes:

y m, ≈ n H m,n, x n, . (7)
The model is therefore separable along the spectral dimension. With the generalization of multi-processor computers or multi-core processors, this property of the operator H may be exploited to parallelize the code to apply H (or its adjoint H ) in a very straightforward way. Using a true matrix representation, H would have a block diagonal structure. Formally, the coefficients of the multi-spectral version of H are given by:

H m,n, =    + cos(θ n • b m /λ ) if y m, is the real part of a measured complex visibility -sin(θ n • b m /λ ) if y m, is the imaginary part of a measured complex visibility (8)
In order to speed up computations, the monochromatic operator H . . can be applied using the NU-FFT algorithm.

Regularization by assuming spatial sparsity

In this paper, we focus on a particular type of astronomical targets which consist in group of point like sources. These sources may have different spectra. For instance, this is suitable for multiple stars or group of stars like in globular clusters or in the center of our galaxy. For such objects, the most effective means to regularize the problem is to favor spatially sparse distributions, that is images with as few as possible sources on a dark background. In this section, we derive expressions of the regularization term f prior (x) suitable to favor spatially sparse distributions.

Separable sparsity

It is now well known that good means to impose the sparsity of the solution in the parameter space (which means that most model parameters will be equal to zero) while approximating the data, is to use the 1 norm as the regularization term:

11 f sparse (x) = x 1 def = k, |x k, | = 1 • x , (9)
where the last equality occurs because of the positivity of the parameters. In principle, any f prior (x) = x p p , with x p the p norm of the sought parameters and 0 ≤ p ≤ 1, can be used to get a sparse distribution with p = 0 corresponding the most to a true sparsity prior. However p = 1 is the smallest p which yields a convex criterion and therefore guarantees a unique solution. In practice, the 1 norm is very effective to obtain a sparse solution. [START_REF] Donoho | For most large underdetermined systems of linear equations, the minimal ell-1 norm near-solution approximates the sparsest near-solution[END_REF] The regularization term f sparse (x) in Eq. ( 9) is completely separable. Hence, in our framework where the model is spectrally separable, the criterion defined in Eq. ( 1) is also separable along the spectral dimension. Thus the image reconstruction can be solved independently for each spectral channel. We therefore do not expect as good results in this case compared to reconstructions performed while imposing some spectral continuity. We shall however use this fully separable regularization in our tests to demonstrate the benefit of imposing non-separable priors along the spectral dimension.

Note that, if the wavelength samples used in the discrete model x of I λ (θ) do not coincide with the effective wavelengths of the observed spectral channel, spectral interpolation of the model has to be used to match the observed wavelengths. In this case, the model implicitly set some spectral correlation and the 3-D image reconstruction has to be performed globally even if the regularization does not impose any kind of spectral continuity.

Non-separable spatial-only sparsity

The following regularization favors spatial sparsity but with sources at the same location whatever the wavelength:

6 f group (x) = n x 2 n,
1/2 (10) where n is the spatial index (pixel) and is the spectral channel. The fact that such a regularization favors spatial sparsity and spectral grouping is a consequence of the triangular inequality inequality. [START_REF] Fornasier | Recovery algorithms for vector valued data with joint sparsity constraints[END_REF]

Means to impose spectral continuity

Some authors [START_REF] Soulez | Restoration of hyperspectral astronomical data from integral field spectrograph[END_REF][START_REF] Bongard | 3-D deconvolution of hyper-spectral astronomical data[END_REF][START_REF] Thiébaut | Maximum a posteriori planet detection and characterization with a nulling interferometer[END_REF] have shown the efficiency of exploiting the spectral continuity of the sources by using a spectral regularization like:

f spectral (x) = n µ n (x n, +1 -x n, ) 2 (11)
with µ n > 0 suitable normalization weights.

A rather extreme means to impose spectral continuity is to assume that the spectral energy distribution is everywhere flat. Using the regularization in Eq. ( 11), this correspond to the limit µ n → ∞, ∀n, and yields:

x n, = g n , ∀ (12)
where g is a gray image of the object. To speed up the reconstruction, only the gray image has to be reconstructed, using the model:

y m, ≈ n H m,n, g n , (13)
and, to impose the spatial sparsity, the following regularization:

f prior (g) = g 1 = n |g n | = 1 • g , (14)
since g ≥ 0.

Optimization algorithm

MiRA algorithm 1 was designed for minimizing a smooth cost function. For that purpose, limited memory quasi-Newton methods such as VMLM-B 17 are relatively efficient and easy to use as they only require the cost function and its gradient to be computed. When using non smooth regularizations to impose spatial sparsity as the ones in Eq. ( 9) and Eq. ( 10), optimization algorithms based on Newton method (that is, on a quadratic approximation of the cost function) are inefficient and we had to implement a completely different optimization strategy to solve the problem in Eq. ( 1) with non-differentiable cost functions. In our algorithm, we introduce variables splitting [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] to treat the two terms of our cost function as independently as possible and exploit the alternating direction method of multipliers [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] (ADMM) to solve the resulting constrained problem.

Introducing auxiliary variables z, minimizing the two-term cost function in Eq. ( 1) can be recast in the following equivalent constrained problem:

min x∈X,z {f data (z) + µ f prior (x)} s.t. x = z . (15)
Imposing that x ∈ X (rather than z ∈ X) is not arbitrary and our motivation for that choice will be clarified in what follows.

The augmented Lagrangian 20 of this problem writes:

L ρ (x, z, u) = f data (z) + µ f prior (x) + u • (x -z) + ρ 2 x -z 2 2 , (16)
with u the Lagrange multipliers associated with the constraints that x = z, ρ > 0 the quadratic weight of the constraints, and v 2 the Euclidean ( 2 ) norm of v.

The alternating direction method of multipliers [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] (ADMM) consists in alternatively minimizing the augmented Lagrangian for x given z and u, then for z given x and u, and finally updating the multipliers u. This scheme, adapted to our specific problem in Eq. ( 15), is detailed in the following algorithm with the convention that v (t) is the value of v at iteration number t.

Algorithm 1. Resolution of Eq. ( 15) by alternating direction method of multipliers. Given ρ > 0, choose initial variables z (0) and Lagrange multipliers u (0) . Then repeat, for t = 0, 1, 2, . . . until convergence, the following steps: a/ update variables x: with:

x (t+1) = arg min x∈X L ρ x, z (t) , u (t) = arg min x∈X f prior (x) + ρ 2 µ x -x (t) 2
x (t) = z (t) -u (t) /ρ ; (17)
b/ update variables z: with:

z (t+1) = arg min z L ρ x (t+1) , z, u (t) = arg min z f data (z) + ρ 2 z -z (t) 2
z (t) = x (t+1) + u (t) /ρ ; (18)
c/ update multipliers u:

u (t+1) = u (t) + ρ x (t+1) -z (t+1) (19)
In Appendix A, we demonstrate how the updating of the variables x (step a of the algorithm) can be done efficiently and taking into account the constraint that the parameters are positive. This is our motivation for imposing x ∈ X on the variables x and not on the variables z to which the likelihood term f data (z) is applied. Another possibility could have been to introduce additional auxiliary variables for the purpose of accounting for the feasible set. The formulae to update the variables x ∈ X for f sparse (x) and f group (x) are respectively given by Eq. ( 28) and Eq. ( 34).

The updating of variables z is equivalent to solving a linear problem of which the analytical solution is:

z (t+1) = H • W • H + ρ I -1 • H • W • y + ρ x (t) + u (t) (20)
with I the identity matrix (of suitable size). Since the augmented term is diagonal, the spectral separability is preserved and this problem can be solved independently for each spectral channel:

z (t+1) = A -1 • b (t) (21)
with:

A = H • W • H + ρ I (22) b (t) = H • W • y + ρ x (t) + u (t) (23)
where the index denotes the sub-vector or the sub-operator restricted to the coefficients corresponding to the th spectral channel. In practice and owing to the number of parameters, we (approximately) solve these problems by means of the conjugate gradients algorithm 21 with left-hand side matrix A and right-hand side vector b (t) .

One of the important settings of the ADMM method is to choose properly ρ: if it is too small, the constraints x = z will converge too slowly; while, if it is too small, the cost function will decrease too slowly. The most effective schemes consist in tuning ρ at every iteration of ADMM. [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]

RESULTS

To check the proposed algorithm, we simulated a cluster of 50 stars with random positions and luminosities and with spectra randomly taken from the library provided by Jacoby et al. [START_REF] Jacoby | A library of stellar spectra[END_REF] We used 100 spectral channels from λ = 493 nm to λ = 507 nm by steps of ∆λ = 0.14 nm. The field of view is 128 × 128 pixels of size 0.5 milliarcseconds. To simulate the observations, we took 100 random interferometric baselines with a maximum baseline of 180 m. We added Gaussian white noise to the complex visibilities with a level such that the maximum signal-to-noise ratio (SNR) is equal to 100.

For the image reconstructions, we considered three different cases: the reconstruction of a multi-spectral distribution with the regularization f sparse (x) in Eq. ( 9), or the regularization f group (x) in Eq. ( 10), and the reconstruction of a gray object g with the regularization f sparse (g). In order to set the relative weight of the priors, we choose the value of the hyper-parameter µ which yields an image which has the least mean quadratic differences with the true distribution. Once chosen the values of µ and ρ, the reconstruction of a 128 × 128 × 100 distribution from ∼ 2 × 10 4 measurements takes a few minutes on a GNU/Linux workstation with a quad-core processor at 3 GHz and using a multi-treaded version of FFTW [START_REF] Frigo | The Design and Implementation of FFTW3[END_REF] to compute the discrete Fourier transforms.

Figure 1 shows the integrated flux, i.e.

x n, , for the true distribution and for the reconstructed ones. In any cases, the spatial sparsity priors effectively yield a solution with point-like structures. However, when there is no trans-spectral constraints, only a few sources are correctly found and there are many more spurious sources. When using f group (x) or assuming a gray object the estimated integrated luminosity is much more consistent with that of the true object: all existing sources are found and the spurious sources are not only less numerous but also much fainter than the true ones.

Figure 2 shows the spectra of the two stars encircled by boxes in Fig. 1, clearly the spectra recovered with f group (x) (blue curves) are of much higher quality than the spectra estimated when treating the spectral channels independently (red curves). Compared to the true spectra (green curves) there is however a small but significant bias in the spectra obtained with f group (x). This is not unexpected as the mixed norm implemented by f group (x) results in a soft compression as shown by Eq. (32). There are different possibilities to avoid this bias. The simplest one consists in using the distribution reconstructed with f group (x) to select the significant sources and then restrict the reconstruction to be non-zero only for these sources (there is no needs for other regularization) REFS. According to the good estimation of the integrated luminosity obtained by imposing a gray object, this may be used to select the location of the sources.

DISCUSSION AND PERSPECTIVES

We have shown the importance of using trans-spectral constraints to improve the quality of the restoration of the multi-spectral brightness distribution I λ (θ) of an astronomical target from optical interferometric data. These results confirm what has been observed for other type of data (like integral field spectroscopy).

For the moment, our demonstration is restricted to specific objects which are spatially sparse (e.g. point-like sources) and must be generalized to other type of spatial distributions. Being implemented by non-differentiable cost functions, spatial sparsity requires specific optimization algorithms. We have shown that variable splitting by the alternate direction method of multipliers (ADMM) is suitable to solve the optimization problem in a short amount of time. In addition to being able to deal with non-differentiable criteria, the ADMM method leads to splitting the full problem in sub-problems that are easier to solve and that may be independent. This straightforwardly gives the opportunity of speeding up the code, e.g. by means of parallelization. This possibility remains if other priors are used, e.g. to account for a smooth spatial distribution.

To simplify the problem at hand, we considered that complex visibilities have been measured. At optical wavelengths, this is only possible with phase referencing. In order to process most existing interferometric data, we will have to modify the likelihood term f data (x) and use a non-linear method (i.e. not the linear conjugate gradients) to update the auxiliary variables z. The new algorithm that we proposed, because it treats separately the two cost functions, f data (z) and f prior (x), may however be an efficient alternative to the variable metric method used in MiRA and which treats directly the sum of the cost functions.

APPENDIX A. PROXIMITY OPERATORS FOR SPATIAL SPARSITY AND STRICT POSITIVITY

Updating of the variables in the ADMM method consists in solving problems of the form:

min x∈X f (x) + α 2 x -x 2 2 . (24)
Solving problem (24) is very close to the job of the so-called proximity operator of the function f (x) which is defined as:

18 prox f ( x, α) def = arg min x∈R N f (x) + α 2 x -x 2 2 . (25)
Proximity operators for non differentiable cost functions like f sparse (x) or f group (x) have already been derived [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] and we simply modified them to account for the additional constraint that x ∈ X. Since X is the subset of non-negative vectors of R N , we denote by: prox

+ f ( x, α) def = arg min x∈R N , x≥0 f (x) + α 2 x -x 2 2 , (26)
the modified proximity operator that we use to update the variables x in our algorithm while accounting for x ∈ X.

A.1 Separable Sparsity

The proximity operator for the 1 norm is the soft thresholding operator:

prox 1 ( x, α) n, =            x n, -1/α if x n, > +1/α ; x n, + 1/α if x n, < -1/α ; 0 else . (27)
Imposing the positivity is straightforward and yields:

prox + 1 ( x, α) n, =    x n, -1/α if x n, > +1/α ; 0 else . (28)
This shows that if 1/α ≥ max n, x n, , the output of the proximity operator is zero everywhere.

A.2 Spatio-Spectral Regularization

First we note that the function f group (x) given by Eq. ( 10) is separable with respect to the pixel index n. Thus all computations can be done independently for the spectrum of each pixel. For variables x such that the function f group (x) is differentiable, the partial derivatives of f group (x) are:

∂f group (x) ∂x n, = x n, β n with β n def = x 2 n, . (29)
Assuming for the moment that β n > 0, the problem of minimizing f group (x) is convex and the solution is obtained by finding the root of the partial derivatives:

∂ ∂x n, f group (x) + α 2 x -x 2 2 =
x n, β n + α (x n, -x n, ) = 0

=⇒ x n, = x n, 1 + (α β n ) -1 (30)
Solving for β n yields:

β n = x 2 n, = x 2 n, 1 + (α β n ) -1 =⇒ β n = x 2 n, - 1 α . ( 31
)
The non-differentiable case occurs when the expression above yields a value of β n which is not strictly positive, that is when x 2 n, ≤ 1/α, in which case the minimum of the cost function is given by x n, = 0, ∀ . Finally, the proximity operator of f group (x) is:

prox fgroup ( x, α) n, =        1 -α x 2 n, -1
x n, if α x 2 n, > 1 ; 0 else.

(32)

When α

x 2 n, > 1, the factor 1 -α x 2 n, -1

is strictly positive and smaller than one, hence the name soft compressor given to this operator.

Requiring that x n, > 0 implies that β n > 0, since α > 0, the solution given by Eq. (30) can only be strictly positive if x n, is itself strictly positive; otherwise, x n, must be equal to zero. Thus:

x n, = max(0, x n, ) 1 + (α

β n ) -1 , (33)
where max(0, x) is applied element-wise. The rest of the reasonning is similar than with no positivity constraints. The proximity operator of f group (x) modified to account for positivity is finally: prox + fgroup ( x, α) = prox fgroup (max(0, x), α) .

(34)

2

 2

2

 2

Figure 1 .

 1 Figure 1. Integrated flux for the star cluster. From top to bottom, left to right: true object; reconstruction with spatial sparsity and assuming a gray object; reconstruction with fully separable sparsity prior; reconstruction with semi-separable sparsity prior. The spectra of the sources encircled by the boxes are shown in Fig. 2.

Figure 2 .

 2 Figure 2. Spectra of two selected sources. Each panel show the spectra of one of the sources encircled by the boxes in Fig. 1. In green: the true spectra; in red : the spectra in the reconstruction with fully separable sparsity prior; in blue: reconstruction with semi-separable sparsity prior.

ACKNOWLEDGMENTS

This work is supported by the French ANR (Agence Nationale de la Recherche), Éric Thiébaut works for the MiTiV project (Méthodes Inverses de Traitement en Imagerie du Vivant, ANR-09-EMER-008) and Ferréol Soulez is funded by the POLCA project (Percées astrophysiques grâce au traitement de données interférométriques polychromatiques, ANR-10-BLAN-0511).

Our algorithm has been implemented and tested with Yorick (http://yorick.sourceforge.net/) which is available for free.