Benjamin Camus
email: benjamin.camus@inria.fr

Fanny Dufossé
email: fanny.dufosse@inria.fr

Anne-Cécile Orgerie
email: anne-cecile.orgerie@irisa.fr

A stochastic approach for optimizing green energy consumption in distributed clouds

Keywords: Data centers, distributed clouds, energy efficiency, renewable energy, scheduling, on/off techniques

The energy drawn by Cloud data centers is reaching worrying levels, thus inciting providers to install on-site green energy producers, such as photovoltaic panels. Considering distributed Clouds, workload managers need to geographically allocate virtual machines according to the green production in order not to waste energy. In this paper, we propose SAGITTA: a Stochastic Approach for Green consumption In disTributed daTA centers. We show that compared to the optimal solution, SAGITTA consumes 4% more brown energy, and wastes only 3.14% of the available green energy, while a traditional round-robin solution consumes 14.4% more energy overall than optimum, and wastes 28.83% of the available green energy.

INTRODUCTION

The rapid increase of demand for Internet services leads Cloud providers to build more and more data centers for hosting these services. The data centers that constitute the Cloud infrastructures are usually geographically distributed for security reasons or to offer lower latency for their clients. This infrastructure increase comes with a dramatic growth of the power consumption globally drawn by data centers. As an example, in 2014, data centers in the U.S. consumed an estimated 70 billion kWh, representing about 1.8% of total U.S. electricity consumption [START_REF] Shehabi | United States Data Center Energy Usage Report[END_REF].

To reduce this impact, Cloud providers resort to renewable energy sources which are either on-site or off-site [START_REF] Tripathi | Optimizing Green Energy, Cost, and Availability in Distributed Data Centers[END_REF]. Such energy sources are mostly intermittent by nature (wind, sun, etc.) with high variations, and periods of time without any production (during night for instance for photovoltaic panels). Energy storage devices can help to overcome this issue. But, they constitute a costly investment and they intrinsically lose part of the energy stored [START_REF] Wang | Energy storage in datacenters: What, where, and how much?[END_REF]. Thus, without storage, renewable energy has to be consumed upon production or it is wasted. In this context, optimizing renewable energy consumption requires to know local availability for the distributed cloud infrastructure, in order to adequately allocate computing resources to incoming user requests. The goal is to geographically distribute the workload among the data centers so that, it fits at best the on-site renewable energy production, which is variable and not known in advance.

In this paper, we consider the problem of scheduling workload across multiple data centers for minimizing renewable energy loss. To solve this issue, we propose SAGITTA: a Stochastic Approach for Green consumption In disTributed daTA centers. SAGITTA uses a stochastic approach for estimating renewable energy production, and a greedy heuristic for allocating resources to the incoming user requests. Our simulation-based results show the efficiency of SAGITTA compared to classical allocation approaches. Indeed, compared to the optimal solution, SAGITTA consumes 4% more energy overall, and wastes only 3.14% of the available green energy, while a classical round-robin solution consumes 14.4% more energy overall than optimum, and wastes 28.83% of the available green energy.

The remainder of the paper is structured as follows. Related work is presented in Section 2. Section 3 details the SAGITTA approach. A simulationbased evaluation is conducted, simulation conditions are described in Section 4 and results are provided in Section 5. Future work is discussed in Section 6.

RELATED WORK

Cloud infrastructures consist in geographically distributed data centers which are linked through com-munication networks [START_REF] Wang | Scientific Cloud Computing: Early Definition and Experience[END_REF]. With the emergence of the Future Internet and the dawning of new IT models such as cloud computing, the usage of data centers, and consequently their power consumption, increases dramatically. As an example, for 2010, Google used 900,000 servers which consumed 1.9 billion kWh of electricity [START_REF] Koomey | Growth in Data Center Electricity Use 2005 to 2010[END_REF]. Other major Cloud companies present similar figures and similar issues [START_REF] Katz | Tech Titans Building Boom[END_REF].

Virtualization technology and its ability to pool resources through transparent sharing should have minimized worldwide data center consumption. But, the energy consumption of state-of-the-art servers grows inexorably as they embed more and more powerful cores and advanced features and technologies. Consequently, the global data center consumption keeps increasing rapidly [START_REF] Shehabi | United States Data Center Energy Usage Report[END_REF]. This situation raises major environmental, economic and social concerns.

The first way to save energy at a data center level consists in locating it close to where the electricity is generated, hence minimizing transmission losses. For example, Western North Carolina, USA, attracts data centers with its low electricity prices due to abundant capacity of coal and nuclear power following the departure of the region's textile and furniture manufacturing (Greenpeace, 2011). This region has three super-size data centers from Google, Apple and Facebook with respective power demands of 60 to 100 MW, 100 MW and 40 MW (Greenpeace, 2011).

Other companies opt for greener sources of energy. For example, Quincy (Washington, USA) supplies electricity to data facilities from Yahoo, Microsoft, Dell and Amazon with its low-cost hydroelectrics left behind following the shutting down of the region's aluminum industry (Greenpeace, 2011). Several renewable energy sources like wind power, solar energy, hydro-power, bio-energy, geothermal power and marine power can be considered to power up super-sized facilities.

In spite of these approaches, numerous data facilities have already been built and cannot be moved. Cloud infrastructures, on the other hand, can still take advantage of multiple locations to use green sources of energy with approaches such as follow-the-sun and follow-the-wind [START_REF] Figuerola | Converged Optical Network Infrastructures in Support of Future Internet and Grid Services Using IaaS to Reduce GHG Emissions[END_REF]. As sun and wind provide renewable sources of energy whose capacity fluctuates over time, the rationale is to place computing jobs on resources using renewable energy, and migrate jobs as renewable energy becomes available on resources in other locations. However, the migration cost, in terms of both energy and performance, may be prohibitive [START_REF] Callau-Zori | An experiment-driven energy consumption model for virtual machine management systems[END_REF].

Within the data center itself, a range of tech-nologies can be utilized to make cloud computing infrastructures more energy efficient, including better cooling technologies, temperature-aware scheduling [START_REF] Fan | Power provisioning for a warehouse-sized computer[END_REF], Dynamic Voltage and Frequency Scaling (DVFS) [START_REF] Snowdon | Power Management and Dynamic Voltage Scaling: Myths and Facts[END_REF], and resource virtualization [START_REF] Talaber | Using Virtualization to Improve Data Center Efficiency[END_REF]. The use of Virtual Machines [START_REF] Barham | Xen and the Art of Virtualization[END_REF] brings several benefits including environment and performance isolation; improved resource utilization by enabling workload consolidation; and resource provisioning on demand. Nevertheless, such technologies should be analyzed and used carefully for actually improving the energy-efficiency of computing infrastructures [START_REF] Miyoshi | Critical power slope: understanding the runtime effects of frequency scaling[END_REF].

Concerning green energy integration, Ren et al. have proposed an online scheduling algorithm which optimizes the energy cost and fairness among different data centers subject to queuing delay constraints [START_REF] Ren | Provably-Efficient Job Scheduling for Energy and Fairness in Geographically Distributed Data Centers[END_REF]. While their work is based on a distributed Cloud model similar to ours, they aim at minimizing the cost of the consumed electricity, instead of the wasted renewable energy in our case. Tripathi et al. have presented a mixed integer linear programming formulation for capacity planning while minimizing the total cost of ownership [START_REF] Tripathi | Optimizing Green Energy, Cost, and Availability in Distributed Data Centers[END_REF]. Their model schedules demand considering the availability of green energy and its price variation to lower the total cost of ownership. Finally, a literature review of renewable energy integration in data centers can be found in [START_REF] Deng | Harnessing renewable energy in cloud datacenters: opportunities and challenges[END_REF].

SAGITTA

In this section, we present our approach named SAGITTA: a Stochastic Approach for Green consumption In disTributed daTA centers. First, the Cloud model and assumptions are described in Section 3.1. Then, Section 3.2 proposes the problem formulation. The details for computing the expected green and brown consumption are provided in Section 3.3. Finally, SAGITTA's algorithms are presented in Section 3.4.

Cloud model

We consider a distributed Cloud infrastructure comprising several data centers geographically distributed and powered by the regular electrical grid on one side and on-site photovoltaic panels (PV) on the other side. The user management of the Cloud is assumed to be centralized. This Cloud model is shown on Figure 1.

Incoming users requests can arrive at any time. Each request requires to be computed by a virtual ma-

Problem formulation

We consider a system of M data centers spread over a large area. A data center DC i is characterized by its number S i of servers. Servers are considered homogeneous over the different data centers, in term of computing capabilities and energy consumption.

As for the application model, we consider identical VMs submitted at unpredictable rate. The VMs are supposed to be executable in less than one time slot. We can thus describe both computing and memory requirement of VMs by the number C of VMs that a server can complete in a single time slot. We consider that a server consumes at full capacity a power of P s .

Finally, the energy consumption of a data center DC i is proportional to its number of servers ON at current time slot t, U i (t). The total power consumed by the system is thus

M ∑ i=1 P s ×U i (t).
This power requirement is to be compared with the green power produced at each data center. We model the green power available at time slot t in data center DC i as a random variable G i (t) that follows a truncated normal distribution of mean Eg i (t) and standard deviation p i (t), with lower limit 0. Thus, the brown power consumed at time slot t in DC i is equal to max(P s ×U i (t) -G i (t), 0). Our problem consists in allocating VMs to data centers, in order to minimize the consumption of brown energy. VMs are allocated by time slots. Then, our objective is to turn ON the adequate number of servers on the better locations for this criteria. We denote N(t) the number of waiting VMs at time slot t. We thus need to have enough servers ON for all waiting VMs at time slot t:

M ∑ i=1 U i (t) ≥ N(t)/C.
All these notations are summarized in Table 1.

Expected green and brown consumption

We now evaluate the expected brown power consumption of data center DC i at time t with n S servers ON, Ec i (n S ,t). We first evaluate the density function of the random variable of the green power generation of

DC i at time t G i (t).
Let X be a random variable following a normal distribution of parameters Eg i (t) and p i (t), density function

φ(t) = 1 p i (t) √ 2π e -1 2 t-Eg i (t) p i (t)
2 and distributive function

Φ(t) = 1 2 1 + erf t -Eg i (t) p i (t) √ 2 .
Then, for x > 0,

P(G i (t) < x) = P(X < x|X > 0) = P(0<X<x) P(X>0)
and the density function of G i (t) equals

φ i (t) = φ(t) P(X > 0)
.

Let B i (t) be the random variable of the brown consumption of DC i at time slot t. For the sake of simplicity, we denote P = n S × P s the power consumed by DC i at time t. Then, We now compute this integral:

Ec i (n S ,t) = E(B i (t)|G i (t) ≥ P)P(G i (t) ≥ P) +E(B i (t)|G i (t) < P)P(G i (t) < P) = E(B i (t)|G i (t) < P) × P(G i (t) < P) = P - P 0 xφ i (x)dx P(G i (t)<P) × P(G i (t) < P) = P × P(G i (t) < P) -P 0 xφ i (x)dx = P × P(0<X<P) P(X>0) - P 0 xφ(x)dx P(X>0) = P × Φ(P)-Φ(0) 1-Φ(0) - P 0 xφ(x)dx P(X>0)
P 0 xφ(x)dx = P 0 x p i (t) √ 2π e -1 2 x-Eg i (t) p i (t) 2 dx = P 0 1 p i (t) √ 2π (x -Eg i (t)) e -1 2 x-Eg i (t) p i (t) 2 dx + P 0 Eg i (t) p i (t) √ 2π e -1 2 x-Eg i (t) p i (t) 2 dx = p i (t) √ 2π   -e -1 2 x-Eg i (t) p i (t) 2 P 0   +Eg i (t)P(0 < X < P) = p i (t) 2 (φ(0) -φ(P)) + Eg i (t)(Φ(P) -Φ(0))
Finally, we obtain:

Ec i (n S ,t) = (P-Eg i (t)) Φ(P) -Φ(0) 1 -Φ(0) -p i (t) 2 φ(0) -φ(P) 1 -Φ(0) , with φ(x) = 1 p i (t) √ 2π e -1 2 x-Eg i (t) p i (t) 2 , P = n S × P s and Φ(x) = 1 2 1 + erf x-Eg i (t) p i (t) √ 2 .

Algorithms description

Our SAGITTA approach uses several algorithms to take decisions and allocate VMs to physical servers. These algorithms are designed to determine at any time slot on which data center to turn ON and OFF servers. At each time slot, our constraint is to turn ON the minimum number of servers that allows for executing all requested VMs, that is N(t)/C . The general algorithm (Algorithm 1) is designed as follows. It first determines if the number of servers available is under or over the requested number. If there is not enough servers ON, Algorithm 2 determines the location of servers to switch on. If some servers are unnecessary, Algorithm 3 determines where servers should be shut down. These decisions are done regarding the expected green energy

Algorithm 1: General algorithm if ∑ 1≤i≤M U i (t) < N(t) C then Switch on decision; (Algorithm 2) else if ∑ 1≤i≤M U i (t) > N(t) C then Switch off decision; (Algorithm 3) end if Transfer decision; (Algorithm 4) for 1 ≤ i ≤ M do
Let U i (t) servers on and fill them, switch off the rest; end for production in the different data centers. More precisely, Algorithm 2 compares the expected extra cost in brown energy consumption c i induced by an additional server ON on any datacenters, and selects the data center with minimum expected extra cost. The variable U i (t) is then incremented, but the servers are only switched on at end of Algorithm 1, when all decisions are taken on any data centers. The same way, Algorithm 3 selects one by one the servers to switch OFF to maximize the expected gain.

Finally, Algorithm 4 evaluates if the brown power consumption could be reduced by transferring the available processing power from one data center to another. More precisely, the algorithm determines some location where a fixed number of servers is turned off, and a new location where the same number of servers is turned on. One server is selected for switch OFF on the data center of maximum gain and another one to switch ON on the data center of minimum cost, if the gain on the first data center exceed the cost on the second one. However, an excessive Algorithm 2: Switch on decision

for 1 ≤ i ≤ M do if U i (t) < S i then Compute c i = Ec i (U i (t) + 1) -Ec i (U i (t)); else c i = C × P s + 1; end if end for while ∑ 1≤i≤M U i (t) < N(t) C do Find j such that c j = min 1≤i≤M c i ; U j (t) + +; Recompute c j ; end while Algorithm 3: Switch off decision for 1 ≤ i ≤ M do if U i (t) > 0 then Compute g i = Ec i (U i (t)) -Ec i (U i (t) -1); else g i = -1; end if end for while ∑ 1≤i≤M U i (t) > N(t) C do Find j such that g j = max 1≤i≤M g i ;
U j (t) --; Recompute g j ; end while turnover of workload between data centers could degrade the quality of the proposed solution. To avoid this, an additional criteria Z is added. Varying cases for this criterion are tested in Section 5.3.

After running Algorithm 4, general Algorithm 1 applies all these decisions. The selected number of servers are turned ON and OFF and all VMs are allocated to available servers.

VALIDATION FRAMEWORK

We evaluate our algorithm through a modeling and simulation (M&S) process. In the following, we first give an overview of the whole cloud implementation model (Section 4.1). We then detail our implementation of the data centers (Section 4.2), of the green power production (Section 4.3), of the cloud workload (Section 4.4), of the algorithm implementation (Section 4.5), and the different simulations performed (Section 4.6).

Algorithm 4: Transfer decision

for 1 ≤ i ≤ M do Compute g i ; Compute c i ; end for while Z and max 1≤i≤M

g i > min 1≤ j≤M c j do Find k such that g k = max 1≤i≤M g i ; Find l such that c l = min 1≤ j≤M c j ; U k (t) --; U l (t) + +; Recompute g k ; Recompute c l ; end while

Simulation overview

The whole cloud implementation model is described in Figure 2. We simulate data centers using the DCSim (Data Center Simulator) discrete-event M&S tool [START_REF] Tighe | DCSim: A data centre simulation tool for evaluating dynamic virtualized resource management[END_REF]. This simulator provides the power consumption of each data center as a function of time.

We implement our algorithm in an ad-hoc way using the Java language into a simulated cloud controller. This simulator receives as inputs the green power production for each data center as well as the cloud workload (i.e. the number of VMs to deploy on the cloud for each time slot). Based on these inputs and on SAGITTA's algorithms, the controller generates for each server the VM allocation and the instructions which are directly sent to the simulated data center manager.

Note that we do not explicitly model the brown power production as we assume it to be infinite (at the scale of the cloud). We also ignore the telecommunication network as we assume it to have negligible impact on the system functioning (we assume network to be oversized for our scenario), and an almost constant power consumption over time if no energysaving technique is applied [START_REF] Orgerie | A Survey on Techniques for Improving the Energy Efficiency of Large-Scale Distributed Systems[END_REF]. Finally, we do not take into account here the energy consumed by the data centers' cooling systems.

In order to perform the simulations, we connect all these heterogeneous models using the MECSYCO (Multi-agent Environment for Complex-SYstem-COsimulation) M&S platform (Camus et al., 2016a;[START_REF] Camus | MECSYCO: a Multi-agent DEVS Wrapping Platform for the Cosimulation of Complex Systems[END_REF] which is based on the DEVS (Discrete-EVent System specification) formalism [START_REF] Zeigler | Theory of Modeling and Simulation: Integrating Discrete Event and Continuous Complex Dynamic Systems[END_REF]. We have defined a DEVS interface for DCSim, and implemented it in MECSYCO.

Data center simulation

Our cloud consists in five homogeneous data centers composed of five clusters. Each of these clusters contains 80 homogeneous nodes, so overall, the cloud comprises a total of 400 servers. The characteristics of each server are based on the Taurus servers of the French experimental testbed Grid'50001 . These Taurus servers are equipped with 2 Intel Xeon E5-2630 CPU with 6 cores each, 32GB memory, 598GB storage and a 10 Gigabit Ethernet interface. In order to determine the power consumption of each node, we implement the power model of [START_REF] Li | Opportunistic Scheduling in Clouds Partially Powered by Green Energy[END_REF], which is based on real measurements made on Taurus nodes. These measurements notably state that a Taurus server consumes 8W when powered OFF, 97W when idle, and 220W at 100% CPU load (i.e. P s = 220W for our algorithm).

Within this cloud, we deploy homogeneous VMs that are equivalent to the Amazon EC2 "large" flavor2 -i.e. each VM requires 4 CPU cores, 8GB memory and 80GB storage. Hence, three VMs can be simultaneously running on one node. For the sake of simplicity, we assume that, when deployed, a VM always works at full capacity. In the same way, we neglect the delays for the VM to start/stop. All the VMs are automatically deleted at the end of each time slot. A time slot lasts five minutes in our simulations.

Green power production

In order to feed the controller during the simulation, we use real recordings of green power production and real workload traces. We get the former from the Photovolta project3 of the University of Nantes. These recordings correspond to the power produced by a sin-gle Sanyo HIP-240-HDE4 photovoltaic panel updated every five minutes over one week. In order to have heterogeneous trajectories between data centers (and thus to represent solar irradiance differences between sites spread across a country), we select recordings starting at different dates, namely: 4th of September 2016, 2nd of February 2014, 8th of June 2014, 22nd of June 2015 and 21st of December 2014. We consider here that 30 photovoltaic panels (for a surface of 165.6m 2 ) are installed at each data center. Then we scale these photovoltaic signals accordingly.

Workload input

We use the normalized ClarkNet HTTP trace of [START_REF] Tighe | DCSim: A data centre simulation tool for evaluating dynamic virtualized resource management[END_REF] for our cloud workload, shown in Figure 3. This workload trace spans over one week. We scale this workload to 98% of the cloud total capacity (i.e. the maximal workload peak represents 98% of the total computing capacity of the cloud). The trace peaks are synchronized with the photovoltaic signal ones to have proper day-night cycles in our simulation.

Algorithm implementation

The controller implementing our SAGITTA approach is run at each time slot (i.e. each five minutes). It saves all the data received from the green power sources during the current day. The controller computes at each time slot the standard deviations p i (t) using this history. It computes each expected green power production Eg i (t) by averaging a reference green power production trajectory (the Photovolta project recording of the 20th of August 2013 in our case which is the day with the best yield) scaled according to the last green power production received from i. More precisely, we denote P re f (t) the green power production at corresponding hour the day of reference (see Figure 4). We obtain the following formula:

Eg i (t) = max 0, PV i (t -1) + P re f (t) -P re f (t -1) 2 .
Note that we consider with this formula that Eg i (t) is equal to the average between the green power production received at t -1 and the one estimated at t. Thus, we take into account that the green power trajectory changes during the time slot, and not only at its beginning.

For implementing the transfer decision algorithm, we use the reallocation constraint:

Z = (i < 1000) ∧ max 1≤ j≤M g j -min 1≤k≤M c k > 1 ,
with i, the number of transfers (i.e. while-loop iterations) already performed for that time slot.

In order to minimize the number of ON/OFF cycles for the servers, the controller fills in priority the hosts already ON. Therefore, from a time slot to the next one, the controller keeps trace of the employed servers.

Figure 4: Expected green power production computation for a time slot from t to t + 1.

Simulated approaches

We compare SAGITTA performance against two Round-Robin inspired algorithms:

• Round-Robin-VM distributes the VMs fairly between the data centers regardless their green power production.

• Round-Robin-DC starts filling with VMs the first data center (in an arbitrary predefined order). If this data center becomes full, the algorithm starts using the next one, and so on.

Like SAGITTA, these two algorithms employ in priority the nodes already ON.

As the performance of Round-Robin-DC strongly depends on the order of the data centers, we test two opposite configurations corresponding to the best and the worst possible contexts. To define these contexts, we sort the photovoltaic traces according to the total amount of green energy they provide. We assign then the traces to the data centers following this order. The best context corresponds to the case where the phototraces are sorted in a decreasing order. Thus, the first data center (i.e. the one filled in priority) will be supplied by best photovoltaic power trajectory. The worst context corresponds then to the case where the traces are sorted in an increasing order (i.e. the data center with the worst green power supply will always be filled first).

To properly evaluate the performance of the three algorithms, it is important to note that the green power available is not always sufficient to supply the cloud needs in our simulation. That is why we also compute the optimal cumulative brown energy consumption which corresponds to the best performance reachable regarding our cloud configuration. We determine this optimal performance based on the optimal brown power consumption of the cloud at time t, P B (t) given by the equation:

P B (t) = max P tot (t) -∑ 1≤i≤M min(PV i (t), P S × S i ) , 0 ,
with :

• P tot (t) the total power consumption of the cloud at time t, • PV i (t) the photovoltaic power production for DC i at time t.

Results

Based on the simulation framework described in the previous section, several experiments were run to validate our proposed approach. First, simulations are conducted without switching ON/OFF costs (i.e. switching does not take time nor energy) to evaluate the allocation algorithms in an ideal context (Section 5.1). An optimal theoretical lower bound is determined this way. Then, new simulations are performed with switching costs in order to fairly compare SAGITTA against state-of-the-art approaches (Section 5.2). The influence of the transfer parameter Z is analyzed (Section 5.3), as well as the influence of the green energy forecast (Section 5.4). Various green production scenarios are studied to estimate the impact of green energy location on SAGITTA's performance (Section 5.5). Finally, the scalability of SAGITTA is evaluated by increasing the number of data centers (Section 5.6).

Without switching ON/OFF costs

We simulate the cloud behavior over one week. First, the power costs of switching ON/OFF the servers are not integrated in order to have a fair comparison with the ideal unreachable case (given by P B (t)) which does not take into account these costs. Our simulation estimates that this cloud consumes a total of 4.96 MWh over the simulated week. Figure 5 shows the cumulative brown energy consumption of the cloud over time for the previously described scheduling algorithms. SAGITTA presents a consumption 4% above the ideal, and significantly better than Round-Robin-VM (28.8% above the optimal) and Round-Robin-DC (14.4% above the optimal in the best case, and 69.6% in the worst case). As shown in Table 2, SAGITTA stands out even more clearly from the other algorithms when considering the percentage of available green energy they waste. We compute these percentages based on the ratio of available green power wasted at time t, W (t) given by the following equation:

W (t) = ∑ 1≤i≤M min(P i (t) -PV i (t), 0) -P B (t) P tot (t) -P B (t)
With P i (t) the power consumption of DC i . It is worth noting that, due to its transfer decision (i.e. Algorithm 4), SAGITTA switches ON/OFF significantly more nodes than the other algorithms: 33,792 switches ON for SAGITTA against 29,606 switches ON for the other algorithms. This difference on the number of switches should have an impact on the overall cloud power consumption. This effect is not visible in this first set of simulations (shown in Figure 5) as they ignore the nodes powering OFF/ON costs.

To sum up, compared to the ideal allocation, SAGITTA consumes 4% more brown energy and wastes 3.14% of green energy (while the ideal allocation does not waste any). For both criteria, green energy waste and brown energy consumption, SAGITTA outperforms traditional approaches based on round-robin allocation.

With switching ON/OFF costs

The second set of simulation integrates the switching ON/OFF costs and estimates their impact on the algorithms' energy consumption to reflect this point. Following the data collected by [START_REF] Rais | Impact of Shutdown Techniques for Energy-Efficient Cloud Data Centers[END_REF] on the Taurus cluster, we add a static energy consumption penalty of 5.28 Wh (consumed in 150 seconds) for each switch-ON command, and 0.56 Wh (consumed in 10 seconds) for each switch-OFF command sent. As shown in Table 3, even when considering these penalties, simulations show that SAGITTA performs better than the other solutions with a difference of at least 10%. Figure 6 shows the power consumption over time of each data center in the simulated cloud using SAGITTA. This figure also shows the number of transfers made by Algorithm 4 -a negative (respectively positive) value meaning that the algorithm switches off (respectively on) hosts. This plot highlights the usefulness of the transfer algorithm. For instance, at time 173,700 s. which corresponds to early morning, DC 2 starts producing green energy slightly earlier than DC 0. SAGITTA takes then advantage of this situation by performing 24 transfers from DC 0 to DC 2. Transfers are highly correlated with discontinuities in the green power production trajectories. Thus, the transfer decision may enable adapting the VM allocation, and consequently the energy consumption, to unforeseen increases and decreases of the green power production. In the absence of transfer, the switch on and off decisions enable adapting the DC workload to their green power productioni.e. the data centers with higher power production are generally more used than the others.

For the sake of simplicity, in the following, we will consider the best case for the Round-Robin-DC algorithm (with data centers ranked by their overall green energy production). All the simulations in the next sections also include the switching ON/OFF costs.

Influence of the transfer parameter

Now, the influence of the transfer parameter is studied. When comparing with the previous simulations, for SAGITTA's case, the switching costs add 6.5% of the overall consumption. Concerning the difference between SAGITTA's power consumption and the other ones, the difference is reduced when taking into account the switching costs. This situation comes from the transfer decisions, and in particular from Z, the transfer decision criteria (used in Algorithm 4).

We redefine as follows the Z constraint in order for the transfer decision to take into account the switching energy costs:

Z = (i < 1000) ∧ max 1≤i≤M g i -min 1≤ j≤M c j × 300 3600 > α
with α, the average brown energy cost of a transfer. Figure 7 compiles the results of 74 simulations using different values of α. It shows that even when considering the switching ON/OFF penalties, SAGITTA performs better for all the α values with at least 2.77 MWh (and 2.75 MWh at best) of brown energy consumed against 3.38 MWh for Round-Robin-VM, and 3.02 MWh (respectively 4.4 MWh) for Round-Robin-DC in the best (respectively worst) context. However, one can note that, the transfer decision has a relatively small impact on SAGITTA overall performance: at best, it only saves up to 3.04 kWh of brown energy, and performs transfers only 5% of the time (in this case, it performs an average of 6 transfers per time slot). Moreover, we observe that, in the absence of a precise estimation of the green energy production, it is safer to overestimate α: then the risk is to lose the small benefit of the transfer decision. At the opposite, when α is underestimated, the transfer decision may degrade SAGITTA performance -i.e. it increases the brown power consumption (up to 19.83 kWh at worst) by inducing too many transfers. Thus α value is not inconsequential and should be properly set if transfers are considered.

Influence of the green energy forecast

One basis of the SAGITTA approach is the green energy production forecast. The value Eg i (t), namely the expected PV production in DC i at time slot t is computed regarding the electricity production at time slot t -1. This approach permits a simple computation for the value Eg i (t) to parametrize the probability law of green energy production. However, this formula estimates the electricity production regarding only the previous time slot, despite of the high volatility of solar energy. We experiment in this section an evaluation of Eg i (t) on a sliding window of PV production values. We target here the optimal size of the window, and the weight to give to the values of the different time slots of the window.

We propose several solutions to determine Eg i (t) on a sliding window of size s. For the sake of simplicity, we denote g i (t) = PV i (t) -P re f (t), with P re f (t) the daily production at same hour, the day of reference. We then make a weighted average value of values g i (t), with weight p i :

Eg i (t) max   0, PV i (t -1) + ∑ s k=1 (g i (t-k)×p s-k ) ∑ s k=1 p k + P re f (t) 2   .
The first variant CST1 uses constant weigths p k = 1 for recent and old values. In the second variant ADD1, the values of p k increase linearly:

p k = k + 1.
Finally, the values of p k are multiplied by 2 at each step in PROD1: p k = 2 k . In these variants, the computation includes values corresponding to the night, when PV i (t) and P re f (t) are both null. This impacts the estimation with useless values. Then, in variants CST2, ADD2 and PROD2, all values g i (t) corresponding to P re f (t) = 0 are removed from the computation. Results of these computations are detailed in Figure 8. Denote that in this experiments, the optimal value of α determined in Section 5.3 is applied. The first unexpected result is the very low values of the optimal size of the sliding window. Regardless of the variant, the best size of the window is always 2, with a slight reduction of the brown energy consumed. The good performance of algorithms PROD1 and PROD2 can be related to the large weight given to the earliest production values in the computation. The weight given to early values has indeed a large impact on the variants' performance.

Influence of green energy production

Cloud providers need to adequately dimension their on-site photovoltaic panel farms. This issue involves a trade-off between the financial cost of installing and operating photovoltaic panels, and the financial gains they are bringing in terms of green energy produced and thus, electricity that has not to be bought from the regular grid. We perform a set of experiments to determine influence of green energy production on SAGITTA performance. As shown in Figure 9, the number of photovoltaic panels (PV) varies per data center and the total brown power consumption is recorded over one week. We can see that, as soon as green energy is available, SAGITTA consumes clearly less brown energy than the other approaches.

Figure 9: Influence of green energy production on brown energy consumption.

Figure 9 also shows that up to about 25 photovoltaic panels, the brown energy consumption curves have a steeper slope, leading to higher gains per photovoltaic panels. For more than 25 photovoltaic panels, the energy gains are lower per added panel. When reaching 45 panels, the green energy production exceeds the total energy consumption of the data center (represented by the case with 0 panel). However, this production is concentrated during the day (as shown in Figure 6), whereas the workload, and consequently the energy consumption, spans over the day and the night. Thus, when reaching a number of photovoltaic panels whose production covers most of the Cloud energy consumption during daylight, adding panels can only save the energy consumption peaks at the beginning and the end of the day (when panels produce less energy), and their buying cost can thus exceed the monetary gains they generate.

Scalability of SAGITTA

In order to check if the SAGITTA's energy savings scale up, we simulate the power consumption of distributed clouds with a larger number of data centers. For these different clouds, we progressively increase the number of data centers, and so the number of green power sources (still taken from the Photovolta project), while maintaining the same total number of nodes (and so an unchanged input workload). The total photovoltaic energy production is also kept as steady as possible by progressively decreasing the number of photovoltaic panels per data centers. Yet we decided not to consider fractions of panels, so the number of panels slightly varies between the scenarios to keep whole numbers. The compositions of these clouds are summed up in Table 4. As shown in Figure 10, the simulation results discloses that SAGITTA scales up: it maintains its energy gains in larger clouds, and always consumes less brown energy than the other approaches. From a computing time point of view, in our simulation environment, it takes 9 seconds to execute SAGITTA over the whole workload trace (representing one week) for the case with 5 data centers, and 28 seconds for the case with 40 data centers. While this computing time is increased by a factor of 3 (when increasing the data center number by a factor of 8), it still remains inconsequential for the scalability of SAGITTA.

CONCLUSION

In this paper, we propose SAGITTA: a Stochastic Approach for Green consumption In disTributed daTA centers. It aims at allocating virtual machines in an energy-efficient way for a distributed cloud comprising several data centers that are geographically dis- tributed and that embed on-site photovoltaic panels. To reduce brown energy consumption, SAGITTA employs a stochastic approach to estimate the expected green energy consumption and to adequately allocate virtual machines on data centers depending on their green energy production. It also switches off unused servers to save energy, while taking into account the energy cost of switching on and off servers.

We conducted a simulation-based evaluation using real workload traces, wattmeter measurements on testbed servers and real production traces from photovoltaic panels. Traditional approaches do not consider the expected green energy consumption for taking virtual machine allocation decisions in distributed clouds. As a consequence, they may overestimate the green energy availability -and take nonefficient scheduling decisions -or underestimate it, and thus waste this energy. The results show that SAGITTA can allocate virtual machines in a more energy-efficient way than traditional approaches, like round-robin. In particular, it wastes only 3.14% of green energy when not considering the switching on/off energy costs. It also exhibits good results in terms of brown energy consumption with a difference of 4% with the optimal when not considering the switching on/off energy costs. When considering the switching on/off energy costs, SAGITTA consumes 10% more brown energy than the theoretical lower bound, which is the ideal allocation not taking into account the switching on/off energy costs. We study the influence of the green energy production on SAGITTA's energy gains and show that, in all cases, it outperforms traditional approaches. The results also show that SAGITTA can smoothly scale with the number of data centers belonging to the cloud.

We plan to extend this work by considering the impact of network devices on the energy consumption and integrating the ability to dynamically migrate virtual machines from one site to another.

Figure 1 :

 1 Figure 1: Considered cloud model

Figure 2 :

 2 Figure 2: Bloc diagram view of the cloud model.

Figure 3 :

 3 Figure 3: The input workload used in the experiments.

Figure 5 :

 5 Figure 5: Cumulative brown energy consumption of the cloud generated by the different allocation approaches.

Figure 6 :

 6 Figure 6: Power consumption per data center with SAGITTA and transfer decisions.

Figure 7 :

 7 Figure 7: Impact of α estimation on SAGITTA performance.

Figure 8 :

 8 Figure 8: Influence of Eg i estimation.

Figure 10 :

 10 Figure 10: Brown energy consumption of SAGITTA with increasing number of data centers.

Table 1 :

 1 Table of Notations Number of machines ON at current time slot on data center i (output) G i (t) Random variable of the green power produced at time slot t Eg i (t) Expected green power generation at data center i during time slot t (input) p i (t) Standard deviation of green power generation on data center i (input) w Workload portion (number of VM): 0 < w ≤ N(t) (input) Ec i (n S ,t) Expected brown consumption of data center i with n S servers ON at time slot t

	Notation	Definition
	Constants
	M	Number of data centers
	DC i	Data center number i
	S i	Number of servers in data center i
	C	Maximum number of VMs in parallel on a server
	P s	Maximum power consumption of a server
	Variables
	N(t)	Number of incoming VMs for time slot t (input)
	U i (t)
	Parameters
	Z	Constraint for reallocation

Table 2 :

 2 Percentage of total available green energy wasted.

		SAGITTA	Round-Robin-VM	Round-Robin-DC
	Best	3.14%	28.83%	14.15%
	Worst	3.14%	28.83%	70.27%

Table 3 :

 3 Total cumulative brown energy consumption when including switching ON/OFF energy costs.

		SAGITTA	Round-Robin-VM	Round-Robin-DC
	Best	2.77 MWh	3.38 MWh	3.02 MWh
	Worst	2.77 MWh	3.38 MWh	4.4 MWh

Table 4 :

 4 The considered cloud scenarios with increasing number of data centers.

	Number of data centers	5	10	15	20	25	30	35	40
	Total number of nodes	400	400	400	400	400	400	400	400
	Number of photovoltaic panels per data centers	30	14	9	7	6	5	4	3

https://www.grid5000.fr

https://aws.amazon.com/ec2/

http://photovolta2.univ-nantes.fr

ACKNOWLEDGEMENTS

The authors would like to thank Yunbo Li for the energy traces of real datacenter servers. The authors would also like to thank Matthieu Simonin and Nathalie Bertrand for their proofreading of the mathematical proofs. This work has been supported by the Inria exploratory research project COSMIC (Coordinated Optimization of SMart grIds and Clouds).


Dr Madleen Sallmann

Fabio Oldenburg

Dr Beatrice Braun

Marius Réglier

Ariane Jalila Simaan

Prof. Dr Christian Limberg
email: christian.limberg@chemie.hu-berlin.de

B. Sc F Oldenburg

Dr M Regliér

Dr J Simaan

A Structural and Functional Model for the 1-Aminocyclopropane-1-carboxylic Acid Oxidase

Keywords: aminocyclopropane-1-carboxylic acid, oxidase, iron, enzyme model, oxygen

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Ethylene is a gaseous hormone for plants that induces their ripening. [START_REF] Abeles | Ethylene in Plant Biology[END_REF] It is produced from the amino acid 1-aminocyclopropane-1-carboxylic acid, through oxidation with O2 as the oxidant and ascorbate as the coreductant, mediated by the enzyme 1-aminocyclopropane-1-carboxylic acid oxidase (ACCO). [START_REF] Adams | Proc. Natl. Acad. Sci[END_REF] CO2 and HCN are generated concomitantly (Scheme 1). For not well-understood reasons, the enzyme also requires CO2 (or bicarbonate ions) to develop its reactivity.

Scheme 1. The biosynthesis of ethylene.

The crystal structure of the ACCO from Petunia hybrida has been determined in 2004 by Schofield and coworkers. [START_REF] Zhang | [END_REF] It belongs to the family of non-heme iron enzymes featuring a structural motif termed the 2-His-1-carboxylate facial triad: The active center contains an iron(II) ion coordinated by two histidine moieties and one aspartate residue in a facial arrangement, leaving the other three sites of the iron ion vacant for the binding of ACC and potentially O2. The mechanism, by which the complex conversion of Scheme 1 is realized, is still discussed controversially, but it is generally assumed that the oxidation of ACC into ethylene follows a radical mechanism including two successive single electron oxidation steps. [4] Spectroscopic studies have revealed that in the first step of the reaction the substrate ACC binds to the Fe II ion in a bidentate mode most likely followed by O2 activation to generate an Fe III -superoxide intermediate. [5] For the events after O2 binding several suggestions have been made that all involve an Fe IV =O intermediate but differ in the sequence of the electron transfer steps with ascorbate and in the nature of the intermediates responsible for the successive oxidation of ACC. Rocklin et al. have proposed that binding of ACC and O2 is followed by one electron reduction to provide an iron(III)-hydroperoxide intermediate that performs the first ACC oxidation step leading to an Fe IV =O species responsible for the second oxidation step. [6] One electron from ascorbate would then complete the catalytic cycle. Alternatively, Mirica et. al. have proposed an ascorbatedependent formation of the Fe IV =O intermediate that would act as the first oxidizing species. [7] Many questions thus remain unanswered regarding the interaction mode with the substrate, the role of the different cofactors/co-substrates (ascorbic acid, dioxygen and carbon dioxide) and the catalytic mechanism.

Bioinorganic model studies can provide valuable information, both with regard to structural issues and requirements for reactivity. However, there are hardly any compounds known which may be regarded as models of ACCO, applying the minimum condition, that they contain ACC and that they display substantial ACCO-like activities. [8,9] One iron complex has been reported, which is dinuclear, though, contains iron(III) ions and two bridging ACC ligands. [8] While iron(II) complexes have long remained unknown, mononuclear 1:1 metal-ACC complexes have been accessed for copper(II). [9] An Fe II -ACC complex has very recently been described, but it can hardly be considered as a functional replicate of the ACCO: the iron centre is coordinated octahedrally by six donor atoms, and hence oxidants have not access to the iron centre. This explains that no reactivity with O2 was observed and that the yield of ethylene produced in reactions with H2O2 was only 7% higher as compared to the blank experiments. [10] To mimic mononuclear non-heme iron dioxygenases based on a 3-His structural motif we have in recent years successfully employed the tris(pyrazolyl)borato (Tp) ligand system, [11] which in the past has also been employed already to simulate representatives based on the 2-His-1-carboxylate structural motif. [12] It thus appeared a promising choice also for the development of replicates of the ACCO. One problem that has emerged in the past with regard to the modeling of this enzyme is the tendency of ACC to bridge two metal centers, so that the model character of the resulting complexes is reduced (vide supra). Hence, we chose phenyl residues at the 3-positions of the pyrazole donors within Tp to create a shielding reaction pocket approaching the situation in the enzyme. A suitable starting material was therefore the complex [Tp Me,Ph FeCl] [13] (Tp Me,Phe = 3-Phenyl-5-methylhydrido-trispyrazol-1-ylborato), which after dissolution in dichloromethane was reacted with the potassium salt of the amino acid 1-amino-cyclopropane-1-carboxylic acid (see Scheme 2). After work-up a light yellow solid was isolated, which was characterized by elemental analysis, IR and 1 H NMR spectroscopy (see Supplementary information). The IR spectrum showed characteristic bands for the NH2 stretching vibrations as well as for the carboxylate absorptions, and the n(BH) band, which is typical for Tp complexes and sensitive to the environment, appeared at 2549 cm -1 . While all analytical and spectroscopic data pointed to the envisaged complex [Tp Me,Ph FeACC], 1, all attempts to crystallize this compound led to crystals of [Tp Me,Ph 2Fe]. To clarify whether the latter complex represents an impurity that has favorable crystallization properties, or whether it forms via ligand exchange during attempts to crystallize 1, it was prepared independently and characterized by 1 H NMR spectroscopy. Comparison showed that [Tp Me,Ph 2Fe] was not present in freshly prepared samples of 1 and neither it was after 7 h, but after storing of solutions for 5 days [Tp Me,Ph 2Fe] could be detected in significant amounts. This clearly showed that 1 is stable enough in solution for reactivity studies but does not readily crystallize, and crystallization can neither be achieved by allotting additional time as this leads to ligand exchange. Mössbauer spectroscopy was used to further characterize the complex. The spectrum of 1 in acetonitrile was recorded and is characterized by quadrupoledublett with an isomeric shift δ = 1.0906 mm/s and an electronic quadrupole splitting of ∆Eq=2.7 mm/s typical of high-spin iron(II) species.

In order to obtain structural information we have thus prepared a corresponding nickel complex [Tp Me,Ph NiACC], 2, in the same way setting out from [Tp Me,Ph NiBr] (see Supplementary information). After work-up of the blue-green reaction solution a light green solid could be isolated, which showed an almost identical IR spectrum as compared to the one displayed by 1, with respect to band shape, pattern and intensities. Since this indicated that the structure of 2 is rather similar to the one of 1, too, a structural investigation of 2 was of high interest. 2 . thf could be crystallized through overlaying of a thf solution with hexane, and Figure 2 shows the molecular structure determined by a single crystal X-ray analysis. As expected the ACC ligand binds as a chelating ligand. The nickel(II) center is thus coordinated by the three N-atoms of the Tp Me,Ph ligand as well as by the amine and carboxylate functions of the O-deprotonated ACC creating a ligand sphere that is inbetween trigonal bipyramidal and square pyramidal (τ = 0.48). While in previous models the C-C bond of the cyclopropyl unit appeared shortened, [9] such an effect was not noticeable in case of 2: The C3-C4 bond length in 2 (1.496(3) Å) was almost identical with the one found for free ACC [14] (1.490 -1.497 Å).

Based on the striking resemblance of the IR spectra exhibited by 1 and 2 a corresponding structure is suggested also for 1, which thus features one free coordination site for the potential binding and activation of O2 or alternative oxidants. [START_REF] Lee | For the proposal/detection of iron(III) superoxide intermediates compare: a)[END_REF] Accordingly, reactivity studies were carried out. It has been suggested that a prerequisite for O2 reactivity is E1/2(Fe III/II ) < -0.1 V vs Fc + /Fc. [16] A CV investigation of 1 dissolved in dichloromethane indicated a reversible redox event, but the oxidation and reduction peaks were separated by 410 mV indicating reversible structural changes upon oxidation/reduction (see Supplementary information). The oxidation peak occurred at 0.03 V, which appeared not negative enough for O2 reactivity. On the other hand reactivity trends observed in previous work for TpFe-based Dke1 models could not be rationalized on the basis of redox potentials, [17] and indeed 1 was found to react with O2.

Complex 1 was dissolved in DMF and added to O2saturated DMF placed in hermetically sealed vials. Analysis of the headspace gas by GC revealed the presence of ethylene, and the conversion yield reached a maximum of ca. 17 % after a few minutes at 70°C. Importantly, no ethylene was detected when ACC alone (in form of its NBu4 + salt) was placed in aerated DMF for one hour indicating that the observed activity is not due to free ACC (potentially released from the complex) in solution. Furthermore, ethylene production was found to be 3 times lower (conversion yields of ca. 6%) when Fe(ClO4)2•xH2O in combination with (NBu4)ACC was employed instead of complex 1 under the same conditions. Notably, the nickel complex 2 does not react with O2, which suggests an initiation of the reaction by binding of O2 at the iron center in the first step.

Under single turn-over conditions and in absence of ascorbate, Rocklin et al. found that only 0.35 mole of ethylene per mole of ACCO are formed, and hence they proposed that the electrons needed for catalysis are provided by a fraction of the initial ACCO enzyme. [6] Also in case of our experiments it is conceivable that the required electrons are provided by a fraction of complex 1 therefore limiting the reaction yield as observed for the enzymatic system. Attempts to provide electrons by adding different reductants (ascorbate or benzoine) remained unsuccessful, though. To obtained information about the fate of the residual skeleton of 1 we have performed GC/MS studies with the gasphase after the reaction, which revealed CO2 as a further product, as one should expect. [19] Monitoring the UV/Vis spectrum of a DMF solution of 1 with time at r.t. and at 70°C showed that after O2 addition the main changes (general increase of absorption in the range of 300-500 nm) occur within 0-5 min. (see SI). After this time the spectrum remains nearly constant displaying only discrete changes, while ethylene yield increases within ca. 2 h. Hence, there is probably an initial fast step (perhaps only partially as part of an equilibrium), followed by a much slower reaction producing ethylene. At r.t. ethylene yield reaches a plateau at ca. 7 %. This suggests that a decomposition reaction competes with the productive path, which, however, can be accelerated to a significant extent by warming to 70°C (à 17%).

For many O2 activating non-heme iron enzymes and models thereof the initial O2 binding step with formation of iron(III) superoxide intermediates has been found to be endergonic in theoretical investigations. The subsequent steps are usually exothermic but in case of low-molecular weight analogues often characterized by substantial barriers, which -together with the endergonic O2 binding step -decrease the reaction rates. If enzymes are considered, which -as the ACCO -in parallel to O2 consume electrons to reach substrate oxidation and feature peroxide intermediates, this problem can be circumvented by employing reduced forms of dioxygen, i.e. for instance H2O2 or Oatom transfer reagents. Significant ethylene production was observed when complex 1 was placed in the presence of PhIO or mCPBA. However, The latter reagents significantly oxidize ACC already in the absence of iron of complex (ca. 75% and 35 % yield within 2-3 h respectively), so that their employment does not provide any information (see SI). Hence, different concentrations of hydrogen peroxide were tested, and interestingly, for a 0.5 mM solution of complex 1 the ACC conversion into ethylene reaches 65% yield employing 10-20 mM of hydrogen peroxide after 20 min. Under the same conditions, ACC alone was hardly oxidized to produce ethylene (less than 2% conversion yield) confirming that the observed activity is not due to released ACC in solution. Also the reactivity of 2 towards H2O2 was found to be rather limited (less than 2% conversion yield), emphasizing the importance of the nature of the metal center. Finally, ethylene production from a 0.5 mM solution of Fe(ClO4)2•xH2O/(NBu4)ACC in DMF in the presence of 10-20 mM H2O2 was found to be 6 times lower (conversion yields of ca. 10-12%) than that observed for complex 1, which demonstrates the importance also of the co-ligand to control the reactivity.

Altogether, these results imply that the observed ethylene production from complex 1 after reaction with O2 or hydrogen peroxide crucially depends on the suitability of the Tp Me,Ph Fe moiety to mimic the 2-His-1-carboxylate iron core within the ACCO. It has been rather well established that H2O2 can react with Fe II or Fe III complexes to provide the corresponding, reactive Fe III -OOH intermediates. [START_REF] Aj | Structure and Bonding[END_REF] Reaction of H2O2 with complex 1 can thus be expected to first result in the oxidation of the Fe II to Fe III and then in the formation of the Fe III -OOH intermediate. In the presence of dioxygen, the formation of the Fe III -OO ° intermediate followed by one electron reduction by a further equivalent of 1 and protonation (for instance through residual water in DMF) could lead to a Fe III -OOH intermediate as well. Our results thus suggest that a Fe III -OOH species is indeed involved in the catalytic cycle of the ACCO and that its formation probably precedes the oxidation of the bound substrate, in agreement with the two proposed mechanisms for the enzymatic system.

In conclusion we have described here the first Fe II -based low-molecular weight analogue which faithfully mimics the structure of the ACCO active site and at the same time also simulates the function. Indeed, this Fe II complex in contact with O2 shows an oxidase activity (17%), which is remarkable as the enzymatic one is only ca. 35% under single turn-over conditions. Future studies, including various alternative oxygenation reagents, electron and proton donors as well as theoretical calculations, will focus on the further elucidation of the mechanism of this biomimetic oxidation.

Entry for the Table of Contents

Scheme 2 .

 2 Scheme 2. Synthesis of [Tp Me,Ph FeACC], 1.

Figure 1 .

 1 Figure 1. Comparison of the IR spectra of 1 (blue) and 2 (red) as recorded for KBr discs.

Figure 2 .

 2 Figure 2. Molecular structure of 2 . thf (all hydrogen atoms have been omitted for clarity). Selected bond lengths (Å) and angles (°): Ni-O1 1.9585(11), Ni-N1 2.0859(14), Ni-N2 2.0993(13), Ni-N5 2.0698(13), Ni-N7 2.0510(13), O1-Ni-N1 82.19(5), N2-Ni-N5 87.29(5), N2-Ni-N7 90.03(5), N5-Ni N7 95.17(5), O1-Ni-N2 95.59(5), O1-Ni-N5 146.53(5), O1-Ni-N7 118.12(5), N1-Ni-N2 175.66(5), N1-Ni-N5 92.51(5), N1-Ni-N7 94.31(5).

  The depicted complex mimics the structure and the function of the ACCO. It reacts with both O2 and H2O2, respectively, to yield ethylene, which allows for inferences with respect to the enzymatic mechanism.Madleen Sallmann, Fabio Oldenburg, Beatrice Braun, Marius Réglier, A. Jalila Simaan*, Christian Limberg* TOC Graphic here))

Acknowledgements

We are grateful to the Humboldt-Universität zu Berlin for financial support. We also thank the Cluster of Excellence "Unifying Concepts in Catalysis" funded by the DFG for helpful discussions. The EU-COST Network (CM1305) on spin state reactivity is acknowledged for support.


Manuel Marin
email: mmarin@ulg.ac.be

Federico Milano

David Defour

Midpoint-Radius Interval-based Method to Deal with Uncertainty in Power Flow Analysis

Keywords: Uncertainty modelling, Power ow analysis, Interval arithmetic, Midpoint-radius

This paper presents a novel method based on midpoint-radius interval arithmetic to deal with uncertainties in the power ow problem. The proposed technique aims at nding a balance between accuracy and computational eciency. It relies on an original decoupling of the interval power ow equations into midpoint and radius parts. This representation allows avoiding the factorisation of an interval Jacobian matrix. Moreover, the proposed formulation is combined with an optimisation problem in order to prevent overestimation of the interval solution while preserving uncertainty. The proposed technique proves to be more ecient than existing approaches based on interval and ane arithmetic and as accurate as the conventional Monte Carlo method.

Introduction

Motivations

Uncertainty is inherent to any physical systems. This is particularly true for power systems, where uncertainty can have several causes, e.g., imprecise demand forecast, price variability, renewable energy generation, economic growth, 5 industry placement, and line aging [1,[START_REF] Cortés-Carmona | Fuzzy arithmetic for the dc load ow[END_REF]. Failing to properly account for uncertainties can, in some cases, lead to erroneous estimations or insecure operating conditions. Therefore, a reliable tool to handle several possible scenarios and combinations of scenarios is crucial to provide a clear understanding of the expected behaviour of the grid. This paper focuses on how to properly account for uncertainties in power ow (PF) analysis.

State of Art

In the literature, uncertainty in PF analysis has been handled mainly by two types of methods: probabilistic and interval-based.

The probabilistic approach relies on solving multiple instances of the PF problem for several (typically randomly generated) possible scenarios, and then aggregate results. The Monte Carlo method is the most common probabilistic approach. The Monte Carlo method is adequate for o-line analysis and is assumed to yield the correct results, provided that a suciently large amount of samples are considered [START_REF] Chen | Probabilistic load ow: A review[END_REF]. However, the computational burden of the Monte Carlo method can be unsuitable for practical purposes, real-time analysis and preventive and/or corrective control actions [START_REF] Su | Probabilistic load-ow computation using point estimate method[END_REF]. For an extensive survey of probabilistic PF methods, the interested reader can refer to [START_REF] Chen | Probabilistic load ow: A review[END_REF] and [START_REF] Billinton | Bibliography on the application of probability methods in power system reliability evaluation 1996-1999[END_REF].

Interval-based methods rely on using intervals to model the system, according to a possibility distribution obtained from experience and historical data.

Sentences such as load between 0.5 and 1 pu and generation around 0.9 pu can be easily translated into intervals. In [START_REF] Wang | Interval arithmetic in power ow analysis[END_REF], the interval Newton method is directly applied to a case of PF analysis with 5 buses, assuming small uncertainty in the nodal injected powers. In [1], the authors use ane arithmetic to keep track of correlation between inputs, a feature that is absent from traditional interval arithmetic. While interval arithmetic has a low computational burden with respect to probabilistic methods, the major drawback is its tendency to overestimate the intervals of the solution, especially if input parameters are characterised by wide intervals. Wide intervals can make the solution either of little practical interest or useless.

Contributions

The technique proposed in this paper deals with uncertainty in PF analysis and balances computational burden and accuracy of results. With this aim, we utilise a midpoint-radius representation of intervals and separate the solution of the PF problem into a standard PF problem for the midpoint and an interval problem for the radius. The latter is solved through a carefully designed optimisation problem, where the constraints ensure that the solution reects all the uncertainty of the input, while the objective function helps to prevent overestimation. The concept of linking the numerical results with uncertainty in the input data is borrowed from backward error analysis [START_REF] Wilkinson | Modern error analysis[END_REF].

The proposed technique enhances the one presented in [START_REF] Marin | Power ow analysis under uncertainty using symmetric fuzzy arithmetic[END_REF] where the midpoint and radius problems were solved together, thus leading to a higher computational burden and lower accuracy than the technique proposed in this paper.

The proposed method is tested on the IEEE 57 bus test case system, proving to yield results as accurate as the Monte Carlo method. A study of the computational burden of analogous methods is performed in order to show that the proposed method is competitive with state-of-the-art interval-based techniques and much more ecient than the Monte Carlo method.

Organisation

The remainder of the article is organised as follows. Section 2 reviews probabilistic and interval-based approaches to deal with uncertainty in PF analysis. Section 3 describes the proposed interval method for PF analysis. Section 4 presents a study of the computational complexity of various methods for PF analysis with uncertainty. Section 5 presents a case-study including a comparison in terms of accuracy with the Monte Carlo method. Finally, Section 6 duly draws conclusions and outlines future work.

Uncertainty in PF analysis

This section reviews two main approaches used for dealing with uncertainty in PF analysis. These are the probabilistic approach and the interval-based approach.

Probabilistic approach

The probabilistic approach models uncertainty as random variables with a certain probability distribution. The above relies on statistical data to obtain the probability distribution of the inputs. A probabilistic PF model is dened by extending the PF equations to random variables. The equations are solved to obtain the distribution of the unknown variables. The solution method can be numerical, such as Monte Carlo method, or analytic. In the following, we focus exclusively on the Monte Carlo method as it is the most commonly used and is considered to be the most accurate approach [START_REF] Chen | Probabilistic load ow: A review[END_REF].

Monte Carlo method. This is a numerical method to approximate the distribution of an unknown random variable. The method relies on the law of large numbers and sampling and consists of the following steps:

1. Create a number of scenarios by taking samples of known random variables.

2. For every scenario, compute a sample of the unknown variables using a deterministic model, e.g., simulation.

3. Aggregate the results into some relevant parameters.

Algorithm 1 illustrates the Monte Carlo method for probabilistic PF analysis.

The algorithm computes the sample mean vector of bus voltages v, given the distribution functions of bus power injections F S (•). The function rng(•) returns a random number in the interval [0, 1], and is used for sampling purposes.

Typically, a high number of samples is needed to achieve accurate results.

For this reason, the method can become cumbersome if applied to real-world problems involved in the operation of power systems. Common applications are power system planning and reliability analysis [START_REF] Billinton | Reliability assessment of large electric power systems[END_REF], or as a validation tool to test other techniques. Several examples can be found of the latter, in which the results from Monte Carlo method are considered the correct ones, e.g., [START_REF] Vaccaro | A range arithmetic-based optimization model for power ow analysis under interval uncertainty[END_REF] and [START_REF] Saric | Ellipsoidal approximation to uncertainty propagation in boundary power ow[END_REF].

Algorithm 1 Monte Carlo method for probabilistic PF analysis.

Input: Number of samples, n s . Distribution functions of bus power injections, F S (•). Deterministic PF equations, s = f (v).

Output: Sample mean of bus voltages, v.

1: for h in 1, . . . , n s do 2: r = rng()

3: s (h) = F -1 S (r) {Sampling} 4: v (h) = f -1 s (h) {Deterministic PF} 5: end for 6: v = 1 n h v (h) {Sample mean}

Interval-based approach

The interval-based approach models uncertainty as intervals without a prob-95 ability distribution. This allows using expert knowledge in the denition of input intervals, in case statistical data is lacking. An interval PF model is dened by extending the PF equations to interval variables. The equations are solved in order to compute interval bus voltages. These methods are self-validated, as interval operations respect the property of isotonicity [START_REF] Alvarado | Uncertainty in power system modeling and computation[END_REF]. 100

Interval-based Newton method. This is a method for bounding the zeros of a dierentiable function f (•) : R n → R n . Given an initial interval guess [x] (0) , the method computes a series [x]

(k) , k = 1, 2, . . ., such that, x ∈ [x] (0) , f (x) = 0 =⇒ x ∈ [x] (k) ⊂ [x] (k-1) ⊂ . . . ⊂ [x] (0)
.

(1)

The method relies on the mean value theorem, which is applied with vectors in the current interval. This leads to a new interval that include all the zeros.

However, the new interval might be overlapping the current one. Therefore, both are intersected in order to compute the next interval in the series. Given

the current interval [x] (k)
, the mean value theorem states that,

f (x) ∈ f (y) + J f x∈[x] (k) (x -y), ∀x, y ∈ [x] (k) , (2)
where J f is the Jacobian matrix of f (•). Note that the Jacobian matrix is evaluated on the whole interval [x]

(k)

. Thus, the result is an interval matrix.

Enforcing f (x) = 0, and choosing y as the midpoint of [x]

(k)

(denoted by x(k) ),
yields the following expression.

x

∈ x(k) -J f x∈[x] (k) -1 f (x (k) ), (3)
which leads to the interval Newton iteration,

[x] (k+1) = [x] (k) x(k) -J f x∈[x] (k) -1 f (x (k) ) . (4)
The interval Newton method is applied to interval PF analysis in the same way as the Newton-Raphson (NR) method is applied to deterministic PF analysis. This idea was introduced in [START_REF] Wang | Interval arithmetic in power ow analysis[END_REF].

Note that the interval Newton iteration requires inverting an interval matrix (or, alternatively, solving a system of interval linear equations). Typical fac-105 torisation techniques include interval Gaussian elimination and the Krawczyk's method [START_REF] Moore | Methods and Applications of Interval Analysis[END_REF]. However, depending on the structure of the matrix, these methods can either not converge, take too long or simply deliver an impractical result [START_REF] Barboza | Towards interval analysis of the load uncertainty in power electric systems[END_REF].

Ane arithmetic method. A novel method for interval PF analysis has been proposed in [1]. The method relies on ane arithmetic, a renement of interval arithmetic which translates intervals into ane forms as follows:

[x] = x + x 1 1 + x 2 2 + . . . + x n n , (5)
where x is the midpoint of [x], the x i 's are known partial deviations, and the i 's are noise symbols representing sources of uncertainty. By denition, x i ∈ R 110 and i ∈ [-1, 1]. In ane arithmetic, the same noise symbol may appear in dierent expressions thus referring to the same source of uncertainty. Therefore, ane arithmetic is able to circumvent the so-called dependency problem, unlike traditional interval arithmetic (see Appendix 7).

In [1], the interval PF equations are re-written using ane forms for bus voltage magnitudes and phases, so that the problem is reduced to the interval computation of the noise symbols through solving a linear system:

A [x] = [y] -[b] , (6)
where A is a real matrix of known partial deviations, [x] is an interval vector of unknown noise symbols, [y] is an interval vector of specied nodal power injections, and [b] is an interval vector of known noise symbols introduced by the ane arithmetic computational process. Note that knowing the noise symbols allows one to compute intervals for all the dependent variables by using the ane forms. The solution of ( 6) is obtained as the solution of two linear optimisation problems, one for the lower bounds and one for the upper bounds.

The method is tested on the IEEE-57 bus test case system, and the resulting intervals are found to be fairly good approximations of the ones computed with the Monte Carlo method, with a little overestimation introduced by the ane arithmetic formulation.

Range arithmetic method. Yet another approach has been recently presented in [START_REF] Vaccaro | A range arithmetic-based optimization model for power ow analysis under interval uncertainty[END_REF], by the same authors of [1], this time using range arithmetic. Once again, the method is tested on the IEEE-57 bus test case system and results are a little more accurate than with ane arithmetic.

Proposed method for interval PF analysis

The proposed method for interval PF analysis relies on the midpoint-radius representation, which allows the PF equations to be solved separately for the midpoint and for the radius. For an introduction to midpoint-radius interval arithmetic, see Appendix 7.

In order to introduce our method, let us rst recall the crisp PF equations as follows:

s i = v i k∈N (Y ik v k ) * , ∀ i ∈ N , (7)
where s i is the complex power injection at bus i, v i is the complex voltage at bus i, Y ik is the element of the admittance matrix linking buses i and k, and N is the set of system buses. Let us also recall the denitions of complex power and voltage, as follows:

s = p + jq, (8) v = v∠θ, ( 9
)
where p is the real power injection, q is the reactive power injection, v is the voltage magnitude, θ is the voltage phase angle, and j is the imaginary unit. Equation ( 7) is translated in interval form as follows:

[s i ] = [v i ] k∈N (Y ik [v k ]) * , ∀ i ∈ N , (10)
where [s i ] and [v i ] are now complex intervals. The proposed method assumes that there is no uncertainty associated with line parameters, so that the element of the admittance matrix Y ik is crisp. Note that this is a common assumption in similar approaches [START_REF] Vaccaro | A range arithmetic-based optimization model for power ow analysis under interval uncertainty[END_REF].

Using the midpoint-radius representation, [s] and [v] are expressed in the following form:

[s] = š, ρs , [START_REF] Saric | Ellipsoidal approximation to uncertainty propagation in boundary power ow[END_REF] [v] = v, ρv , [START_REF] Alvarado | Uncertainty in power system modeling and computation[END_REF] where š = p + j q, [START_REF] Moore | Methods and Applications of Interval Analysis[END_REF] v = v∠ θ, [START_REF] Barboza | Towards interval analysis of the load uncertainty in power electric systems[END_REF] and ρs, ρv ∈ R.

Note that the above denitions imply that both [s] and [v] dene discs in the complex plane. In this way, the uncertainty associated with both real and reactive power is jointly dened by the power radius ρs. Similarly, the uncertainty associated with both voltage magnitude and phase angle is jointly dened by the voltage radius ρv (see [START_REF] Rump | Fast and parallel interval arithmetic[END_REF] for a complete description of complex intervals in the midpoint-radius format).

At this point we adapt traditional notions of crisp PF analysis to interval PF, resulting in the following assumptions:

• the real power injection and voltage magnitude midpoints, p and v, are specied on all generator buses;

• the real and reactive power injection midpoints, p and q, are specied on all load buses;

• the voltage magnitude and phase angle midpoints, v and θ, are specied on the slack bus.

In addition, we introduce the following novel assumptions specically regarding interval PF:

• the power injection radius, ρs, is specied on all system buses;

• the voltage radius, ρv, is specied on all generator buses.

The above set of assumptions is summarised in Table 1. Generator p, v, ρs, ρv θ, q Load p, q, ρs v, θ, ρv

By applying the rules of midpoint-radius interval arithmetic (see Appendix 7), equation ( 10) is split in midpoint and radius parts. This yields the following system of equations:

ši = vi k∈N (Y ik vk ) * , ∀ i ∈ N , ( 15a
)
ρs i = (|v i | + ρv i ) k∈N |Y ik |ρv k + ρv i k∈N Y ik vk , ∀ i ∈ N . (15b)
Note that equation (15a) represents a crisp PF problem on the midpoints of all variables. Solving this problem, e.g., by using the NR method, yields vi , the midpoint of the interval solution. In order to compute the radius ρv i on load buses, we propose the following multi-objective optimisation problem:

min ρv i , ∀ i ∈ N (16a) s.t. (|v i | + ρv i ) k∈N |Y ik |ρv k + ρv i k∈N Y ik vk ≥ ρs SP i , ∀ i ∈ N , ( 16b
)
ρv i ≥ ρv SP i , ∀ i ∈ N P V (16c) ρv i ≥ 0, ∀ i ∈ N P Q (16d)
where N P V is the subset of generator buses, N P Q is the subset of load buses, ρs SP i is the power injection radius specied at bus i ∈ N , and ρv SP i is the voltage radius specied at bus i ∈ N P V .

The problem presented above is the core of our method for interval PF analysis, as follows. Starting from the bottom, constraint (16d) states that the voltage radius on load buses must be non-negative, following a basic assumption in interval algebra. Constraint (16c) states the voltage radius on generator buses must be greater than or equal to the specied, so that the interval PF solution eectively takes into account all the uncertainty on the bus voltages.

Similarly, constraint (16b) states that the computed power injection radius on all system buses must be greater than or equal to the specied, so that the solution eectively takes into account all the uncertainty on the power injections.

Note that the left-hand side of (16b) is directly taken from (15b). Finally, the objective function (16a) intends to prevent overestimation of the solution intervals through minimizing the radius.

We can nd similarities and dierences between the proposed method [START_REF] Gall | Powers of tensors and fast matrix multiplication[END_REF] and the one from [1]. For example, the non-negativity of the voltage radius in our method seems to have the same eect as the bounds on the noise symbols in the method from [1]. A clear dierence is that in the proposed method, the intervals are expressed in the midpoint-radius format, whereas in [1] they are expressed in the lower-upper format. This allows us to solve one optimisation problem rather than two. Another big dierence lies in the representation of complex intervals, as disks in the proposed method, and rectangles in the method from [1].

Complexity analysis

This section studies the time complexity of the proposed method as well as the Monte Carlo, classic interval arithmetic and ane arithmetic methods.

In the remainder of the section, n b denotes the number of buses in the power system.

In the Monte Carlo method, the running time is given by the number of samples multiplied by the running time of each instance of PF analysis. For illustration, let us assume that individual instances are solved through the NR method. At each iteration of the NR method, the LU factorisation of the Jacobian matrix dominates all other calculations with a running time proportional to the time of multiplying two matrices, i.e., O(n 2.376 ), as shown in recent studies [START_REF] Gall | Powers of tensors and fast matrix multiplication[END_REF]. Under the very weak assumption that the number of iterations needed to converge is much lower than the problem size, we obtain a time complexity of O(n s • n 2.376 ), where n s is the number of samples. Generally speaking, n s should be chosen in accordance with the standard deviation of results in order to achieve a target level of condence. However, n s is often determined empirically by looking at the standard deviation of results for increasing number of samples. When the standard deviation ceases to change, it means that enough samples have been taken. For instance, in the case study presented in Section 5, the ideal number of samples for a network of 57 buses is determined to be about 5,000 (i.e., two orders of magnitude larger than n b ). In the Monte Carlo method, n is the size of the PF problem, which is about 2n b . For example, if we assume n s ≈ 100 • n b , the running time of the Monte Carlo method can be estimted as 5.19

• n s • n 2.376 b ≈ 519 • n 3.376 b .
Regarding the classic interval-based Newton method, calculation time is dominated by solving an interval linear system. This can be done with the algorithm described in [START_REF] Rump | On the solution of interval linear systems[END_REF], which has an asymptotic complexity of O(n 3 ) and a running time of 5n 3 . In this case, n is the number of unknowns of the PF analysis, i.e., 2n b . Hence, the running time of the interval-based Newton method is

40 • n 3 b .
The ane arithmetic PF method proposed in [1] requires to solve two multiobjective optimisation problems. For illustration, let us assume that the interior point algorithm is used to solve these problems [START_REF] Hertog | Interior point approach to linear, quadratic and convex programming: algorithms and complexity[END_REF]. This algorithm has a time complexity of O(n 3.5 ), where n is the number of decision variables. In the ane arithmetic PF method described in [1], the number of decision variables is equal to the size of the PF problem, i.e., 2n b . Thus, the running time of the two multiobjective optimization problems can be stimated as 2

• (2n b ) 3.5 ≈ 22.6 • n 3.5 b .
Finally, for the proposed method based on interval analysis, we need to solve It is worth noticing that both our method and the ane arithmetic PF method described in [1] have the same asymptotic time complexity. However, we expect our method to be faster than the ane arithmetic approach, for the following reasons:

• Our method only needs to solve one standard PF problem and one optimisation problem, whereas the method in [1] requires to solve two optimisation problems.

• The number of decision variables in [1] is strictly higher than the number of decision variables in our problem (twice as large, in the worst case).

Table 2 summarises the above discussion. The Monte Carlo method is the slowest one as, typically, one has that n s n b . The interval-based Newton method is the fastest one but its convergence cannot be always guaranteed.

This fact consistently limit applications of the interval-based approach to the 12 PF problem. Finally, the ane arithmetic method and the one proposed in this paper are comparable, with our method being slightly faster.

Case study

This section studies the accuracy of the proposed method in comparison with Monte Carlo and (whenever possible) other approaches from literature.

Experiments are carried on the IEEE-57 and IEEE-118 bus test systems [START_REF]Power systems test case archive[END_REF],

where uncertainty is introduced by assuming that generator and load bus power injections, and generator bus voltages, oscillate within an interval. The midpoint of these intervals is specied as the value from the original case denition, whereas the radius is specied so as to reect an arbitrary level of uncertainty, as discussed below.

The proposed method is implemented in the Julia programming language [START_REF] Bezanson | Julia: A fresh approach to numerical computing[END_REF], using the package for mathematical programming JuMP [START_REF] Lubin | Computing in operations research using julia[END_REF] and the interior-point solver Ipopt [START_REF] Wächter | On the implementation of an interior-point lter line-search algorithm for large-scale nonlinear programming[END_REF]. The multiple objective in (16a) is implemented using a min-max approach, where the maximum voltage radius is minimized.

The Monte Carlo method is implemented using Dome [START_REF] Milano | A python-based software tool for power system analysis[END_REF] and assuming that all inputs follow a uniform distribution over the complex disks dened in Section 3.

IEEE-57 bus test system

Figure 1 shows the one-line diagram of the IEEE-57 bus system, which has 7 generators, 42 loads and 80 branches including transmission lines and transformers. On this system, three uncertainty scenarios are considered as presented in Table 3. Note that scenario (a) corresponds to the case analyzed in [1] and [START_REF] Vaccaro | A range arithmetic-based optimization model for power ow analysis under interval uncertainty[END_REF], where the uncertain variables are assumed to vary over rectangles in the complex plane. For this scenario, results of the Monte Carlo, ane arithmetic and range arithmetic methods are obtained from [START_REF] Vaccaro | A range arithmetic-based optimization model for power ow analysis under interval uncertainty[END_REF]. Furthermore, the midpoint of the interval PF solution in the proposed method is computed as the midpoint of the Monte Carlo solution from [START_REF] Vaccaro | A range arithmetic-based optimization model for power ow analysis under interval uncertainty[END_REF]. For scenarios (b) and (c), the midpoint of the interval solution and the Monte Carlo solution are computed using Dome [START_REF] Milano | A python-based software tool for power system analysis[END_REF]. Figure 2(a) shows the bus voltage magnitude intervals obtained with the proposed method, the ane arithmetic method, the range arithmetic method and the Monte Carlo method for the medium uncertainty scenario (a). We observe that the intervals computed by the proposed method are very close to the ones by Monte Carlo at the vast majority of buses, whereas with the ane and range arithmetic methods there seems to be a consistent overestimation.

Nonetheless, the proposed method is unable to match the 0% oscillation specied for the voltage magnitude at generator buses (namely, buses 2, 3, 6, 8, 9 and 12), whereas this feature is guaranteed by default in the ane and range arithmetic methods. All in all, the performance of the proposed method in this scenario is 14 Voltage magnitud (pu) (a)

Figure 2: Bus voltage magnitude intervals obtained with dierent methods on the IEEE-57 bus test system. (a): medium uncertainty scenario (see Table 3).

Bus number 3). Voltage magnitud (pu) (a) Figure 3: Bus voltage magnitude intervals obtained with dierent methods on the IEEE-118 bus test system for a scenario of medium uncertainty, i.e., 20% oscillation on the generator and load bus power injections, and 1% oscillation on the generator bus voltages. However, in scenario (c) there is a clear overestimation of the interval results at buses of low identier. By looking at Fig. 1, we observe that buses with lower identier are typically closer to generators, which in practice may help keeping the voltage magnitude within controlled bounds. Unfortunately, the proposed method is unable to capture this peculiarity which in turn becomes clear by 290 running the heavy Monte Carlo method.

IEEE-118 bus test system

The IEEE-118 bus test system has 54 generators, 91 loads and 186 branches including transmission lines and transformers. On this system we consider only one scenario of medium uncertainty, i.e., 20% oscillation on the generator and load bus power injections, and 1% oscillation on the generator bus voltages.

Figure 3 shows the bus voltage magnitude intervals obtained with the proposed method and the Monte Carlo method for this scenario. Note that the dierence between the two is very small, which demonstrates the accuracy and scalability of the proposed approach.

Conclusion and future works

The paper proposes a novel method for PF analysis allowing to capture uncertainties through an interval-based approach. The proposed solution proves to be as ecient as but less demanding than methods previously proposed in the literature, such as the interval-based Newton and ane arithmetic approaches.

The proposed method can be as accurate as the Monte Carlo method in situations of moderate uncertainty, with a consistently lower computational burden. Indeed, a case study based on the IEEE-57 and IEEE-118 bus systems indicates that, if the uncertainty is below a given threshold (namely, 20%), the results obtained with the proposed method are very close to those provided by Monte Carlo.

Future works will focus on enhancing the proposed optimisation problem [START_REF] Gall | Powers of tensors and fast matrix multiplication[END_REF] in order to: a) properly deal with high uncertainty levels; b) incorporate stability considerations and reactive power limits.

In addition, tests on larger networks, both academic and non-academic, will also be performed so as to keep investigating the accuracy and scalability of the proposed method.

Interval operations and rounding

Interval operations are dened so as to respect the basic property of inclusion isotonocity, dened below. Definition 2 (Inclusion isotonicity) Let • be a basic arithmetic operation, i.e.,

• ∈ {+, -, •, /}, [x] and [y] intervals. If x • y ⊆ [x] • [y] , ∀ x ∈ [x] , ∀ y ∈ [y] ,
then • is said to be inclusion isotone. Inclusion isotonicity ensures that no possible values are left behind when performing interval operations. If oating-point numbers are used, the above means handling rounding errors inherent to oating-point arithmetic.

The IEEE-754 Standard for oating-point computation allows for self-validated interval arithmetic implementations by dening the following three rounding attributes: rounding upwards (towards innity), rounding downwards (towards minus innity), and rounding to nearest [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF]. The IEEE-754 Standard also denes the relative rounding error and the smallest representable (unnormalised) oating-point number, which are particularly useful when the midpoint-radius representation is used [START_REF] Rump | Fast and parallel interval arithmetic[END_REF][START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF].

Midpoint-radius interval arithmetic requires rounding to nearest when computing the midpoint, and rounding upwards when computing the radius. In addition, the error in computing the midpoint is to the radius in order to satisfy inclusion isotonicity. Midpoint-radius interval operations are dened as follows. It is important to note that both equations (17b) and (17c) introduce a bounded overestimation into the result. In the case of the multiplication, the overestimation is bounded by a factor of 1.5 [START_REF] Rump | Fast and parallel interval arithmetic[END_REF].

The impact of rounding on performance depends on the computing architecture. For example, for most CPU architectures, the rounding attribute is a processor mode, and thus changes in rounding mode cause the entire instruction pipeline to be ushed.

The dependency problem

The so-called dependency problem causes overestimation in interval computations where the same variable occurs more than once. Since all the instances of the same variable are taken independently by the rules of interval arithmetic, the radius of the result might be expanded to unrealistic values. For example, consider the function, f (x) = x 2 -1, where x ∈ [-1, 1]. Then, f (x) ∈ [-1, 0].

However, using interval arithmetic, f

([-1, 1]) = [-1, 1] • [-1, 1] -1 = [-2, 0].
This can become a major issue to the application of interval arithmetic to real-world problems, especially to non-linear ones. Accordingly, specic measures have to be taken to preserve the relevance of the results.

5 b

 5 one multi-objective optimisation problem with n b decision variables. Assuming once again that the interior point algorithm is used, our method has an asymptotic complexity of O(n3.5 ) and an estimated running time of n 3.. To this time, one has to add the time required to solve the crisp PF problem for midpoint values of the intervals, which is n 2.376 b .

  problem size, n b : number of buses, ns: number of Monte Carlo samples.

Figure 1 :

 1 Figure 1: One-line diagram of the IEEE 57 bus test case system.

Figure 2 :

 2 Figure 2: Bus voltage magnitude intervals obtained with dierent methods on the IEEE-57 bus test system. (b): low uncertainty scenario, (c): high uncertainty scenario (see Table3).

Figures 2

 2 Figures 2(b) and 2(c) show the bus voltage magnitude intervals obtained with the proposed method and the Monte Carlo method for the low uncertainty scenario (b) and high uncertainty scenario (c), respectively. Note that the intervals computed by the proposed method for scenario (b) are very accurate. 285

  Midpoint-radius interval arithmetic Intervals are convex sets of real numbers. They can be represented in several formats, e.g., lower bound-upper bound, midpoint-radius, and even lower bound-diameter. This article considers the midpoint-radius representation. Definition 1 (Midpoint-radius representation) Let [x] be an interval dened by [x] = {x ∈ R : |x -x| ≤ ρx}, where x, ρx ∈ R, ρx ≥ 0. Then [x] is noted x, ρx .

Definition 3 (

 3 Midpoint-radius interval operations) Let [x] = x, ρx and [y] = y, ρy . Let (•) and (•) be the rounding attributes to nearest and towards plus innity, respectively. Let be the relative rounding error, = 1 2 and η be the smallest representable (unnormalised) oating-point positive number. Then, [x] ± [y] = z, ( |z| + ρx + ρy) , z = (x ± y), (17a) [x] • [y] = z, (η + |z| + (|x| + ρx)ρy + |y|ρx) , z = (xy),

[ 1 ]

 1 A. Vaccaro, C. Cañizares, D. Villacci, An ane arithmetic-based methodology for reliable power ow analysis in the presence of data uncertainty, IEEE Transactions on Power Systems 25 (2) (2010) 624632.

Table 1 :

 1 Known and unknown variables for each type of bus in the proposed interval PF analysis method.

	Bus type	Knowns	Unknowns
	Slack	v, θ, ρs, ρv	p,	q

Table 2 :

 2 Asymptotic time complexity of dierent methods for PF analysis with inclusion of uncertainty.

	Method	Asymptotic complexity	Running time
	Monte Carlo	O(n s • n 2.376 )	5.19 • n s • n 2.376

Table 3 :

 3 Uncertainty scenarios and specied radius (as percentage of the midpoint) for the IEEE-57 bus system case study.

	Scenario	Uncertainty level	Power injection radius	Voltage radius
	(a)	Medium	20%	0%
	(b)	Low	10%	1%
	(c)	High	50%	2%