Centralized energy optimization at district level Simon Arberet, Yves Stauffer, Max Boegli and Emmanuel Onillon Swiss Center for Electronics and Microtechnology (CSEM) Neuchˆatel, Switzerland [email protected]

Abstract—Optimizing energy flow at district level in order to minimize exploitation costs was carried out in the framework of the European project AMBASSADOR. The proposed solution relies on model predictive control and aims at controlling the battery set-point. The present article focuses on the optimization problem formulation, predictive model development, test condition presentation and associated result analysis. The obtained savings are in the range of 10 to 20%. However, in very favorable cases, savings of up to 80% can be obtained. The provided analysis explains this huge difference. In addition, given the importance of the battery in the current district optimization, it will be shown, that the present methodology, can also be used to find the optimal battery sizing for given boundary conditions.

Building (consumption)

Photovoltaic panels (PV) Connection to grid

Index Terms—district, energy, MPC, battery, renewable. + -

Battery

I. I NTRODUCTION Given the costs of electricity and the penetration of renewable energy sources, we propose in this article a novel method for district energy flow optimization in order for instance to maximize the benefits linked to the usage of storage. First, the test site district in Greece is described. It is composed of: several buildings, renewable productions (photovoltaic and wind turbines) and electric storage (battery), as shown in Figure 1 [1]. Second, a constraint optimization problem is formulated in order to optimize the energy flow at district level. The optimization takes advantage of the usage flexibility provided by the storage elements (battery) as well as of the production/consumption predictions provided by the different district components to perform model predictive control (MPC). Every 15 minutes, the optimal battery set point is computed based on the predicted consumption and production over a horizon of 24 hours. The optimization problem is formulated as a convex optimization problem and solved using CVX [2], [3]. Third, the district components models, which predict the production and consumption for the considered elements are presented. Renewable production elements include photovoltaic and wind turbine. Prediction components of renewable productions implement support vector regression (SVR) [4] between weather conditions (e.g. solar irradiance and wind speed) and produced power over the last 4-7 days. Then the energy production over the next 24 hours is predicted based on This work has been developed under the project AMBASSADOR, that receives funding from the European Union Seventh Framework Programme Grant Agreement No. 314175. c 2016 IEEE 978-1-4673-8463-6/16/$31.00

Wind Turbines (WT)

Figure 1. Greek test site map view

weather forecast. A similar approach is taken for the prediction of energy consumption elements, which include datacenters, chemical plants, administrative buildings and cafeteria. In addition to the weather forecast (e.g. air temperature) the prediction of energy consumption takes into account working schedule and process planning. Fourth, the various simulation conditions and associated results are presented and analysed. We assess the performance of our algorithm in term of money saved in different scenarios, i.e. different electricity tariffs profiles (variable vs flat tariffs, and different ratios between the selling and buying tariffs), and if the battery is or not allowed to charge using the energy from the grid (i.e. only self-storage is allowed). Simulation results are shown in Figure 11, in very favorable cases, savings of up to 80% are obtained. However, savings of 10 to 20% are more common. II. TEST SITE DESCRIPTION In the context of the AMBASSADOR project, energy flow management will be deployed in a test site in Greece. The premises of the Lavrion Technological and Cultural Parc (LTCP) will be used as a final validation scenario. This test site covers about 250’000 m2 and is composed of several buildings as well as standalone equipment including simple consumers, producers (wind-turbines or photovoltaic) and

storage. In addition, some buildings include storage and/or production. Describing the complete test site is not relevant for the presented work, if desired refer to [5], [6]. In order to assess the functionality of the proposed algorithms, a subset of the available buildings will be employed so as to reflect more specific usage cases. In that context, two scenarios were identified. A so called, northern user case, that is based on consumption, storage and wind turbine production and a southern user case, where the production is based on photovoltaic (PV) panels. A. Southern user case The southern case will be illustrated by the administrative building of the LTCP, shown in Figure 2. It is a one storey neoclassic building with a total surface area of about 690 m2 . The electric production is provided by PV (15kWp) and the storage by a battery (14.4kWh).

correspond to the amount of charge and discharge of the battery during each time step of the horizon. By convention we define that a positive value corresponds to charging while a negative value corresponds to discharging. We describe in the following the optimization problem formulation of a single battery for sake of readability, but the extension to multiple batteries was implemented and is straighforward (the cost of each battery is just added in the objective function in a vectorized way). The objective function of problem (1) contains the following terms (also depicted in Figure 3): 1)

:=leakage=(unwanted=discharg)

a)

:=min=/=max=power

b)

:=min=/=max=state=of=charge

+ 2)

:=efficiency=losses

+ 3)

:=usage=cost=(number=of=cycles)

constraints

+ 4)

:=energy=cost=(buy=or=sell)

5)

:=remaining=energy=value

Administrativepbuildingp(NTUA)

+ =

total=cost=(€) Figure 3. The different terms of the objective function and constraints which appear in the optimization problem (1). Topbuilding

Inverter

To/frompgrid

Inverter

Optimizedpset-point

Storage:pbatteries

Production:pPV

Figure 2. administrative building with PV and battery to illustrate the southern case [6]

B. Northern user case For the northern user case, the datacenter combined with wind turbines and storage is used. When compared to the southern case described in section II-A this corresponds to a consumption/production/storage increase of a factor 20. This choice was done in order to simulate a small district rather than a single building. III. OPTIMIZATION The optimization is based on a model predictive control (MPC) approach where the energy flow is optimized over a horizon of 24 hours and a sampling interval of 15 minutes, resulting in M = 96 timeslots. This optimization, which considers all the M timeslots simultaneously (i.e. the decision variable x is a vector of size M ), is rerun every 15 minutes with updated measurements. The only controlable elements of the district being the batteries, the optimization variables

1) the battery leakage cost, where l (no unit) is the battery leakage parameter, and co [e/Wh] is the vector of the electricity (selling) tariffs over the horizon. 2) the battery charging and discharging cost, where e (no unit) is the battery efficiency and where ca [e/Wh] is a cost value which was chosen1 to be the average between ci and co . ci being the electricity tariff vector (over the horizon) of buying from the grid. 3) the battery life cost, where the cost factor u [e/Wh] is computed as the battery cost Cb divided by its max number of cycles Nc , normalized by the energy range of the battery β −α, where α [Wh] and β [Wh] are respectively the minimum and the maximum battery energy state of charge. More formally: u := (Cb /Nc )/(β − α). 4) the cost of buying from and selling to the grid, co [e/Wh] and ci [e/Wh], being respectively the electricity tariffs vector (over the horizon) of selling to, and buying from the grid. d = (d(1), . . . , d(M ))| is the district power balance vector (without including the batteries) in [W]2 , where by convention d(t) > 0 means that the consumption is higher than the production at time sample t and the opposite is true when d(t) < 0. 5) the value of the stored energy at the end of the horizon, where cr [e/Wh] is a weight to promote the charging of the battery to account for what remains in the battery 1 as the tariff of buying and selling are different and vary across the time, it is difficult and somehow arbitrary to value the cost of energy loss. 2 τ := 900/3600 in (1) is a constant to convert power [W] to energy [Wh]. The objective function should be multiplied by this factor because the electricity tariffs are in [e/Wh] and not in [e/W], but as it is a constant, it doesn’t have any influence on the optimization, so we skipped it for simplicity.

at the end of the horizon. This term is here to take into account the future beyond the horizon and avoid the optimization to empty the battery at the end of the horizon. The constraints of problem (1) are: a) the battery charging power max (η [W]) and discharging power max (µ [W]). b) the battery energy state of charge min (α = 0 [Wh]) and max (β [Wh]). 2) efficiency losses

1) leakage

}| z }| { z
{ minimize co | (lLx) + ca | x+ (1 − e) + x− (1 − 1/e) + x∈RM

over the last 4-7 days. The input conditions include weather conditions (e.g. solar irradiance and wind speed) and a time reference (1-24h cycle). Then the energy production over the next 24 hours is predicted based on weather forecast and a time reference. The time reference adds prediction robustness towards uncertain weather forecast. Input/ouput signals for the prediction component models for production elements are listed in Table I. An example of PV power prediction is shown in Figure 4. The achieved root mean square of the prediction error over the next 24 hours is lower than 10% (see Figure 6a). Table I I/O OF PRODUCTION PREDICTION COMPONENT MODELS , INCL . PV S AND W IND T URBINES

ukxk1 + co | (x + d)− + ci | (x + d)+ − cr 1| x {z } | {z } | | {z }

3) usage cost

4) energy cost

I/O

5) remain.

input

subject to:  a)  b) output

solar irradiance (measured) for PV or wind speed (measured) for wind turbine time reference (daily cycle) power production (measured)

solar irradiance (forecast) for PV or wind speed (forecast) for wind turbine time reference (daily cycle) SVR kernel coefficients power production (predicted)

SVR kernel coefficients

Training

50

Prediction

50

Reference Power Predicted Power

Measured Power

40

30

30

P (kW)

40

20 10 0 -7

20 10 0

-6

-5

-4 -3 Time (day)

-2

-1

0

0

6

07-May-2013 16:00:00

1500

12 Time (h)

18

24

1500 Measured Irradiance

Forecast Irradiance

1000

1000 I (W/m 2 )

where: | • x = (x(1), . . . , x(M )) is the battery power vector in [W]. The battery is charging at sample time t when x(t) > 0 and discharging when x(t) < 0. x+ := max(x, 0) and x− := min(x, 0) are respectively the element-wise positive and the negative parts of x. • L is a low triangular matrix of size M × M whose aim is to integrate the power to obtain the state of charge vector. • 1 is a vector of one coefficients used in order to compute the state of charge at the end of the horizon by summing all the charging and discharging battery events. • s [Wh] is the current (initial) battery state of charge. We also implemented an optional constraint (not shown here) which aims to prevent charging the battery with energy coming from the grid. This rule is sometime imposed by the grid in order to prevent the district to reduce its cost by trading electricity. The optimization was implemented in CVX [2].

P (kW)

L(x+ e + x− /e)τ + s ≥ α L(x+ e + x− /e)τ + s ≤ β

(1)

Prediction (24 h)

I (W/m 2 )

x≥µ x≤η

Training (7 days)

500

0 -7

500

0 -6

-5

-4 -3 Time (day)

-2

-1

0

0

6

12 Time (h)

18

24

IV. PREDICTION MODELS The district component models predict the production and consumption of the considered district elements. The production and consumption predictions are then aggregated into the district power balance vector d ∈ RM in the optimization problem (1). The power balance is predicted over a horizon of 24 hours and a sampling interval of 15 minutes, resulting in M = 96 timeslots. Similar to the optimization case, the prediction is updated every 15 min over the same receding horizon of length M . Renewable production elements include photovoltaic and wind turbine. Prediction components of renewable productions implement support vector regression (SVR) between input conditions (e.g. weather conditions) and produced power

Figure 4. PV production prediction. The training is done over the past 7 days and the prediction is estimated for the next 24 hours.

Energy consumption elements include datacenters, chemical plants, administrative buildings and a cafeteria. Prediction components of consumption elements implement SVR between input conditions and consumed power over the last 4-7 days. The input conditions include weather conditions (e.g. air temperature), working schedule, process planning and time reference (daily and weekly cycles). Then the energy consumption over the next 24 hours is predicted based on weather forecast, working schedule, process planning and a time reference. Input/output signals for the prediction component models for consumption elements are listed in Table II.

Table II I/O OF CONSUMPTION PREDICTION COMPONENT MODELS , INCL . CAFETERIA , OFFICES , DATACENTERS , CHEMICAL PLANTS

input

output

Training (7 days)

Prediction (24 h)

air temperature (measured) working schedule process planning time reference (daily / weekly cycles) power consumption (measured)

air temperature (forecast) working schedule process planning time reference (daily / weekly cycles)

8

8 P (kW)

P (kW)

10

4

2

2

Reference Power Predicted Power

0 -5

-4 -3 Time (day)

-2

-1

0

0

6

12 Time (h)

18

24

07-May-2013 16:00:00 1.5

1.5 Daily cycle Weekly cycle

Daily cycle Weekly cycle

Timing

1

Timing

1

0.5

0 -7

0.5

0 -6

-5

-4 -3 Time (day)

-2

-1

0

145

0

6

12 Time (h)

18

130

135

140

145

Figure 6. Normalized root mean square error (NRMSE) for the prediction of a) PV production and b) Cafeteria consumption. The NRMSE is averaged over a 24 hours prediction horizon. The error evolution is shown for 20 days.

6

4

-6

140

Time (day)

Prediction

Measured Power

10

0 -7

135

Time (day) b) Cafeteria Consumption Prediction Error

5

0 125

power consumption (predicted)

12

6

10

130

SVR kernel coefficients

SVR kernel coefficients

Training

12

a) PV Production Prediction Error

5

0 125

NRMSE (%)

I/O

10

NRMSE (%)

An example of the cafeteria consumption prediction is shown in Figure 5. The achieved root mean square of the prediction error over the next 24 hours is lower than 10% (see Figure 6b).

24

Figure 5. Cafeteria consumption prediction. The training is done over the past 7 days and the prediction is estimated for the next 24 hours.

The SVR prediction is implemented with the StatLSSVM toolbox [4] for nonparametric regression using Gaussian kernels. Two parameters need to be tuned in the SVR algorithm: 1) the kernel bandwidth, which defines the regression smoothness or granularity, and 2) a regularization term, which defines the regression robustness or a tradeoff between fitting accuracy and stability. Parameter tuning can be done either manually or by an automatic procedure. V. SIMULATION , RESULTS AND ANALYSIS

The consumption over the 30 days is 13’760 kWh, which corresponds to an average of 459 kWh/day, and the battery size is 500 kWh with a maximum power charging of 200 kW. Thus the battery size corresponds approximately to an average day of consumption, which corresponds approximately to an average day of production. We evaluate the performance of our MPC method over the 30 days of the simulation, with different electricity tariffs profiles (flat tariffs and variable tariffs). In each case, the reference is a simulation of the district without using the battery. 1) Variable tariffs: Figure 7 shows the different tariffs profiles for the variable tariffs case. In every cases, the tariff from grid, i.e. the tariff of buying electricity from the grid is 0.1 euro per kWh during the night and 0.14 euro per kWh during the day. As depicted on Figure 7, we tested three different tariffs to grid (i.e. tariffs of selling electricity to the grid), a null tariff, a tariff equal to the tariff from grid, and a tariff equals to 0.8 times the tariff from grid. For each tariff to grid, we tested two scenarios, one which directly uses the optimization problem (1), and a second which implements a battery constraint to prevent charging the battery with electricity from the grid. It results in six different scenarios which are summed up in Table III. Table III SCENARIOS AND PERFORMANCE IMPROVEMENTS FOR VARIABLE TARIFFS ON THE NORTHERN DISTRICT

A. Northern district with wind turbines (WT) The north district introduced in section II-B is composed of a datacenter building, some wind turbines, and a battery. We simulate this district for a duration of 30 days. The district production (by the wind turbines) over the 30 days is 24’414 kWh, which corresponds to an average of 814 kWh/day.

1

2

3

4

5

6

tariff to grid factor battery constraint

scenario #

0 yes

0 no

0.8 yes

0.8 no

1 yes

1 no

improvement [%]

53

53

5

8

3

22

ation can then exploit the variation of the tariffs in time in order to realize capital gains.

0.15 0.14

2) Flat tariffs: Figure 9 shows the different flat tariffs profiles used for this experiment. In every cases, the tariff from grid is 0.12 euro per kWh during the whole day/night. As depicted in Figure 9, we tested four different tariffs to grid, a null tariff, a tariff equal to the tariff from grid, and tariffs equals to 0.2 and 0.8 times the tariff from grid. For each tariff to grid, we tested the two scenarios with and without battery constraint. It results in eight different scenarios which are summed up in Table IV.

price [euros/kWh]

0.112 0.1 0.08 from grid to grid (ttg=1*tfg) to grid (ttg=0.8*tfg) to grid (ttg=0*tfg)

0.05

0

0

5

10

15

20

0.15

times [hours] 0.12

reference

price [euros/kWh]

Figure 7. variable tariffs Ambassador MPC

1500 1000

Cost (EUR)

500

0.1 from grid to grid (ttg=1*tfg) to grid (ttg=0.8*tfg) to grid (ttg=0.2*tfg) to grid (ttg=0*tfg)

0.05

0.024

0 1

2

3

4

5

0

6

-500

0

5

10

15

20

times [hours]

-1000

Figure 9. flat tariffs

-1500 -2000 -2500

Simulation scenario

Table IV SCENARIOS AND PERFORMANCE IMPROVEMENTS FOR FLATS TARIFFS ON THE NORTHERN DISTRICT

Figure 8. variable tariffs results on the northern district scenario #

• •

when the difference betwen the tariff to grid and the tariff from grid is important as in scenarios #1 and #2, when the tariff to grid is equal to the tariff from grid and there is not the battery constraint. The battery optimiz-

2

3

4

5

6

7

8

tariff to grid factor battery constraint

0 yes

0 no

0.2 yes

0.2 no

0.8 yes

0.8 no

1 yes

1 no

improvement [%]

46

47

60

60

0

0

0

0

reference

Ambassador MPC

1000

500

Cost (EUR)

Variable tariffs analysis: Results depicted in Figure 8 show that in the two first scenarios, where the district cannot sell energy to the grid (tariffs to grid = 0), using the battery can decrease the electricity cost by two. In the remaining scenarios, where the tariff to grid is higher, the district is able to make profits by selling its production to the grid. The usage of the battery optimized by our MPC algorithm increases the profit by 5% (with battery constraint) and 8% (without battery constraint) in scenarios #3 and #4 where the tariff to grid is 0.8 times the tariff from grid, while when the tariff to grid is equal to the tariff from grid (scenarios #5 and #6), the battery usage increases the profit by 3% (with battery constraint) and 22% (without battery constraint). These results highlight the cases where the battery optimization has a large impact:

1

0 1

2

3

4

5

6

7

-500

-1000

-1500

Simulation scenario

Figure 10. flat tariffs results on the northern district

8

B. Southern building with PV panels The southern building introduced in section II-A is composed of an administrative building with PV panels, and a battery. We simulate this district for a duration of 30 days. The district production (by the PV panels) over the 30 days is 1291 kWh, which corresponds to an average of 43 kWh/day. The consumption over the 30 days is 742 kWh, which corresponds to an average of 25 kWh/day, and the battery size is 25 kWh with a maximum charging of 30 kW. Thus the battery size corresponds approximately to an average day of consumption. We evaluate the performance of our MPC method over the 30 days of the simulation, with variable and flat tariff profiles. In each case, the reference is a simulation of the district without using the battery. In both cases, the tariff to grid is equal to zero (as it is the case on the real test site of Lavrion) and the variable and flats tariffs from grid are the ones respectively depicted in Figure 7 and Figure 9. reference

Ambassador MPC

25

20

Cost (EUR)

Experiments with different battery sizes: In this last experiment, we study the influence of the battery size on the result gain, for the case we have a variable tariff from grid and a zero tariff to grid. The optimization gain obtained with different sizes of battery is shown in Figure 12. The gain obtained with the nominal battery size of 25 kWh is 83% and corresponds to scenario #2 of Figure 11. It can be observed that, with a battery of size 12.5 kWh, i.e. twice smaller than 25 kWh, the gain is only reduced by about 10%. Gain versus battery size

100 80

Gain (%)

Flat tariffs analysis: Results depicted in Figure 10 show that in the four first scenarios, where the tariffs to grid are low (0 and 0.2), using the battery can decrease the electricity cost by 46% in scenarios #1 and #2 and, 60% in scenarios #3 and #4. However, in the remaining scenarios, where the tariffs to grid are high (0.8 and 1), the district is making profit and there is no additional profit due from the battery usage. This fact confirms, what was already observed for the variable tariffs, that the impact of the battery optimization is important when there is an important difference between the tariff to grid and the tariff from grid. Note also that in all of these eight scenarios, the absence of the battery constraint did not help, which can be understood by the fact that the battery cannot “trade” the electricity if the tariffs do not vary in time.

83.05 %

60 40 20 0 0

5

10

15

20

25

30

35

Battery size (kWh)

Figure 12. Optimization gain versus battery size. Scenario with variable tariff from grid, zero tariff to grid and without the battery constraint. The gain obtained with the nominal battery of 25 kWh is 83%.

VI. CONCLUSION AND OUTLOOK Reducing exploitation costs of districts by using model predictive based storage element (i.e. battery) control was presented in this article. In particular, the optimization problem was illustrated and the associated test cases highlighted. Emphasis was put on the boundary conditions, which include energy tariffs and battery sizing with respect to daily production and consumption. The simulations and associated analysis clearly show what boundary conditions are favorable for battery deployment and associated smart control. Indeed, in favorable cases, savings up to 80% can be achieved. The present work will be deployed on the Greek test site within the end of 2015.

15

83% reduction

83% reduction

84% reduction

84% reduct.

10

5

0

#1, tfg 2 values, ttg=0.0, with bat. const.

#2, tfg 2 values, ttg=0.0, without bat. const.

#3, tfg flat, ttg=0.0, with bat. const.

#4, tfg flat, ttg=0.0, without bat. const.

Simulation scenario

Figure 11. Costs on the administrative building for variable tariffs (left) and flat tariffs (right). The results are the same with and without the battery constraint.

Results: Results are depicted in Figure 11. The gain obtained with our MPC algorithm is 83% for variable tariffs and 84% for flat tariffs. The influence of the battery constraint is negligible in this district as can be observed by comparing scenarios #2 and #4 with scenarios #1 and #3 respectively.

R EFERENCES [1] “AMBASSADOR Demonstration site,” http://ambassadorfp7.eu/demonstration-sites/demonstratio-site-no-2-lavrion/, accessed: July 2015. [2] M. Grant, S. Boyd, and Y. Ye, “CVX: Matlab software for disciplined convex programming,” 2008. [3] Y. Stauffer, E. Onillion, S. Arberet, E. Olivero, L. von Allmen, and D. Lindelof, “MPC based fan coil unit control algorithm,” in CISBAT, 2015. [4] K. De Brabanter and J.A.K Suykens and B. De Moor, “Nonparametric regression via StatLSSVM,” Journal of Statistical Software, vol. 55, no. 2, pp. 1–21, 2013. [Online]. Available: http://www.esat.kuleuven.be/sista/statlssvm/ [5] M. Taxiarchou, A. Peppas, K. Tsoukalas, L. Karalis, P. Bernaud, J.-L. Bergerand, and R. Khodabuccus, “D1.2 test scenarios for validation sites,” AMBASSADOR FP7 project, Public Deliverable, 2015. [Online]. Available: http://ambassador-fp7.eu/mission-objectives/deliverables/ [6] A. Peppas and K. Tsoukalas, “D3.3 final ambassador communication infrastructure,” AMBASSADOR FP7 project, Public Deliverable, 2015. [Online]. Available: http://ambassador-fp7.eu/missionobjectives/deliverables/