Model Predictive Control for the Energy Management in a District of Buildings Equipped with Building Integrated Photovoltaic Systems and Batteries

: This paper introduces a Model Predictive Control (MPC) strategy for the optimal energy management of a district whose buildings are equipped with vertically placed Building Integrated Photovoltaic (BIPV) systems and Battery Energy Storage Systems (BESS). The vertically placed BIPV systems are able to cover larger areas of buildings’ surfaces, as compared with conventional rooftop PV systems, and reach their peak of production during winter and spring, which renders them suit-able for energy harvesting especially in urban areas. Driven by both these relative advantages, the proposed strategy aims to maximize the district’s autonomy from the external grid, which is achieved through the cooperation of interactive buildings. Therefore, the major contribution of this study is the management and optimal cooperation of a group of buildings, each of which is equipped with its own system of vertical BIPV panels and BESS, carried out by an MPC strategy. The proposed control scheme consists of three main components, i.e., the forecaster, the optimizer and the district, which interact periodically with each other. In order to quantitatively evaluate the benefits of the proposed MPC strategy and the implementation of vertical BIPV and BESS, a hypothetical five-node distribution network located in Greece for four representative days of the year was examined, followed by a sensitivity analysis to examine the effect of the system configuration on its performance.


Introduction
Human-induced climate change is a significant problem that has economic, policy and social dimensions. One of the most important reasons for this phenomenon is the increased fossil fuel consumption over the past half century [1]. Buildings have attracted much attention in the context of decarbonizing the energy system due to the fact that they consume approximately 30% of global electricity and represent approximately 28% of total global energy-related CO2 emissions [2]. To this end, the combination of Renewable Energy Systems (RES) with Battery Energy Storage Systems (BESS) is considered to be an effective solution that not only reduces the energy consumption and carbon footprint of the buildings, but also provides potential for high controllability, dealing with demand diversity and peak shaving issues [3][4][5].
The majority of RES applications in buildings are considered to be photovoltaic (PV) systems [6], due to their high production potential and ease of implementation [7]. The implementation methods of PV panels in a building can be distinguished into two main categories: Building integrated PV (BIPV) and the building attached PV (BAPV) [8]. The additional RES, the renewable energy production is higher, but the proposed technologies are not easily applied in urban areas. In the cases with only one RES, this is of conventional rooftop PV type, and not of vertical PV type. Furthermore, the interest of this work is focused on the publications that are related to a district-level investigation. In this context, similar publications are provided below. The authors of [33] present a novel MPC-based strategy for a BESS station, aiming to improve its economic performance in a power grid with high PV penetration. The presented control strategy is designed from the point of view of the power grid operator and the goal is to decrease the equivalent cost. The forecaster uses the Weighted Mean Average (WMA), which is a simple, yet effective statistic method. The simulation results verify that the presented control strategy achieves great economic performance, effective mitigation of the rapid PV power variation and extension of the system's lifetime, since BESS can maintain a reasonable operating condition. In [34] the optimization model concerns a Peer-to-Peer (P2P) trading environment, including rooftop PV systems and BESS. The objective of the optimization is the minimization of the net energy cost of all households within the planning horizon. The results highlight the economic benefits of distributed household renewable generation in a P2P energy trading market. However, the model presented does not include either a MPC-based strategy or a forecaster or vertical PVs and the objective function is related to cost minimization and not system autonomy. The authors of [35] have developed a decentralized power management scheme based on MPC for the case of microgrids (MG), including PV and BESS. The purpose of the study is to maintain the power balance within the MG and to share load demand among distributed generation (DG) systems, which differentiates it from the objective of system autonomy. Furthermore, there is no special consideration of vertical PVs. The simulation results of the study indicate that the proposed control strategy maintains the voltage and state of charge balance among the BESS. In [36] a power management strategy based on MPC for direct current (DC) MG is introduced. The proposed controller achieves power management by regulating the DC bus voltage to a nominal value and controls the power sources to operate at the maximum power point. A study model including two PV systems, a wind generator (WG) system, one BESS and a DC load are used to demonstrate the effectiveness of the proposed controller. It is noted that the forecast is formulated with the use of Autoregressive Integrated Moving Average (ARIMA) in order to effectively predict the system disturbances. Again, the scope of this work is different to the system autonomy. Furthermore, the integration of wind turbines in the overall concept makes this work not highly representative for the case of urban areas, where the wind turbines cannot be easily installed. Furthermore, the type of PV panels used in this work does not include the vertically placed ones. The authors of [37] have developed an uncertainty-resistant stochastic MPC-based model for a MG, which consists of electrical and thermal units, including PV panels and one BESS. The key target of the controller is to minimize the MG operation cost. The uncertainty variables are represented by typical scenarios. The simulation results demonstrate the proposed strategy's robustness derived from its uncertainty-resistant approach. However, the system does not consider vertical PV panels, introduces the use of only one BESS and incorporates a cost minimization objective function. The authors of [38] present a heuristic method for the energy management of a small grid-connected system that includes vertically placed BIPV panels, wind turbines, electric vehicles (EV), two BESS and buildings. The developed heuristic batteryprotective control is compared with a traditional hierarchical control strategy to present its reliability and robustness in terms of multi-criteria performance improvements, including equivalent CO2 emissions, import cost, energy flexibility and battery capacity. However, the integration of wind turbines surely increases the supply of renewable energy, but it is a scenario that cannot be easily applied in densely populated territories, thus, making the investigation concept not representative for the case of urban areas. In [39] MPC is applied on multiple residential MGs, including PV and storage systems, taking into consideration the thermal and electrical needs of each MG. The purpose of the controller is the cost minimization without considering the autonomy of the MGs from the main grid. Although the controller can achieve its purpose, the implemented strategy does not include the development of a forecaster and limits its innovation to the development of the optimizer. Finally, the PV systems are not vertically placed, which significantly affects the production curve throughout the year, the installed capacity and the applicability of the proposed system in urban areas. Similarly, the authors of [40] propose MPC for the energy management of a distribution network that includes MG subsystems, equipped with PV panels, storage systems and wind generators. The controller aims to minimize the overall grid's cost. Again, the investigation does not focus on the autonomy of the district from the grid and the integration of wind turbines surely increases the supply of renewable energy, thus, examining a scenario that cannot be easily applied in densely populated territories (urban areas). Furthermore, in this case, no forecaster has been developed and no vertically placed panels have been considered, too. In [41] an isolated MG consisting of diesel generators, PV systems, storage systems and other distributed resources is proposed to be controlled by a MPC strategy. The objective of the optimizer is the cost minimization and the strategy only includes the development of the optimizer, without developing the respective forecaster. Therefore, the differences of this work with [41] are the same with the ones mentioned for publication [40]. On the contrary, in [42] an MPC strategy, including both an optimizer and a forecaster is presented. In fact, the developed forecaster utilizes an adaptive autoregression algorithm in order to predict the demand of a grid-connected MG, while the optimizer minimizes the cost, thus, not examining the autonomy of the system from the grid. The MG, where the MPC is applied on, contains PV systems, wind generators and one BESS. The proposed strategy has also been experimentally tested. However, the installation/utilization of wind generators, as already thoroughly described, cannot be performed in urban areas where the concept of energy communities is more feasible to be realized, as the current work considers. The authors of [43] have approached the concept of MPC on MG architectures considering the resilience of the MG as objective of the optimizer. The MG, where the MPC is applied, includes PV systems, one BESS and a variety of distributed resources. The developed forecaster predicts the voltages and currents of the MG, in order to enhance the system's resilience. Nevertheless, the existence of multiple distributed sources increases the provision of renewable energy from technological systems that cannot be easily used in densely populated areas. Furthermore, there is no examination of vertical PVs and the objective function is related to the system stability and not the system autonomy from the grid, which better resembles the concept of energy community. In [44] an existing local distribution system in Osimo, Italy is presented. The district includes PV panels, as well as other distributed resources and considers the possibility of utilizing the BESS of EV, not only as loads, but also as sources in order to maximize the district's autonomy and reduce the fuel based energy consumption. Even though the results are encouraging, the proposed methodology is not based on an optimization algorithm, but the case study is modeled by using a deterministic hour-based simulation model. However, the concept of considering the EV as possible energy sources and not solely as consumers consists in a subject that, in the future, can also be included in the scenario presented in this work as an additional scenario. The maximization of autonomy has also been studied by [45]. In fact, a district of buildings equipped with PV panels on their rooftops and one BESS each is simulated. The purpose of the study is to examine the district from the aspect of autonomy and emissions. Nevertheless, the methodology is not based on a MPC strategy, optimizes the BESS control without taking the possible future parameters into consideration and assumes only rooftop PV panels.
In order to examine the prospects of the BIPV implementation in urban areas, this paper proposes an MPC-based method that is applied on buildings of the same district, each equipped with vertically placed BIPV and a BESS, focusing on the cooperation of these buildings in order to gain as much autonomy from the external grid as possible. Therefore, the major contribution of this study is the management and optimal cooperation of a group of buildings, each of which is equipped with its own system of vertical BIPV panels and BESS, carried out by an MPC strategy. For this purpose, both an optimizer and a forecaster have been developed. The optimizer of the applied controller aims to minimize the energy exchange between the external grid and the interactive buildings with respect to the system's constraints. For the purpose of forecasting the unknown parameters, i.e., the PV energy production and the load of each building, a forecaster has been developed. The forecaster is based on a variation of the WMA method, using Trigg's tracking signal [46]. The selection of the WMA, as a basis of the forecaster, is justified by its extensive use, simplicity and accuracy [47]. Finally, the selection of vertically placed BIPV systems give the opportunity for higher installed capacity and greater feasibility and constructability as regards urban areas. The proposed MPC-based method is applied in a hypothetical distribution network in Greece, consisting of five, (5) interactive buildings including four (4) residences, and (1) one school building, as presented in Figure 1. Four (4) representative days of the year are examined in order to evaluate the cost and energy savings throughout the year and demonstrate the effectiveness of the proposed approach. Finally, two types of sensitivity analysis have been carried out, using different number of PV panels and total number of available BESS, to demonstrate the effect of selected critical design parameters on the system behavior.
The main aspects of this paper are the following:  A new MPC-based method is introduced, aiming to control a district whose buildings are interconnected in a local grid and equipped with vertically placed BIPV and BESS, thus forming a small-scale energy community.  The vertical placement of BIPV systems gives the opportunity to: (a) Consider higher installed capacity due to the provided area as compared with the case of rooftop PV installation, (b) consider a production profile with high peak and total production during the colder months of the year (instead of summer, which is the respective period for rooftop PV systems) due to the angle of the installed panels, and (c) study a system configuration that is more feasible to be implemented in urban areas compared to other renewable sources.  The developed optimizer of the proposed MPC-based method aims to achieve the maximum possible autonomy of the district, resulted by the maximum possible cooperation among the parties.  The optimizer of the controller ensures the parity of the district's five BESS.  The developed statistical forecaster utilizes tracking signals for the detection of anomalies of the unknown parameters and adjusts the forecasted values according to the smoothed error.  The controller is tested under a range of different weather conditions. The results of two types of sensitivity analysis are presented to highlight the system's response to parameter changes.

Methodology of the Proposed Controller
The proposed MPC-based method consists of three main parts, i.e., (a) the forecaster, (b) the optimizer, and (c) the district. A schematic representation of the controller, which is developed in Python, is presented in Figure 2. The three main components interact periodically with each other, considering that the planning horizon is the time interval from the current time step, t, until the end of the day. The loop described in Figure 2 is repeated for 24 times in total, once per hour. The specific tasks of each component of the control system are briefly described below:

 Forecaster
The forecaster's purpose is to provide the optimizer with forecasts of the PV production as well as forecasts of the load of each building. At the beginning of the day, the daily forecasted values are produced based on recent historical data, using the four previous days. At each hour of the day, the actual consumption and production values of the district are measured and imported to the forecaster. If the real values tend to deviate from the forecasted ones, then the forecast of the planning horizon needs to be adjusted accordingly. The measured values and the forecasts are then imported to the optimizer in order to optimize the objective of the model for the planning horizon.

 Optimizer
The optimizer consists of an objective function followed by constraints. In this case, the objective function has been chosen to minimize the daily energy exchange between the external grid and the interactive buildings. By achieving this goal, the maximum feasible autonomy of the interactive buildings is accomplished. The constraints simulate limitations of the power system such as the capacities of the BESS, their technical characteristics, etc. The optimizer calculates the optimal charge/discharge of each BESS, in order to achieve the maximum autonomy of the district for the planning horizon.

 District
The district consists of five buildings, which are connected to the external grid. The set points of the five BESS for each hour of the day are controlled by the optimizer.

Formulation of the Optimizer
The basis of the proposed solution constitutes a Mixed-Integer Linear Programming (MILP) problem, since the model is formulated via continuous and binary variables. In this implementation, the objective is to promote the cooperation of the buildings that are equipped with production (PV) and storage (BESS) units. It is assumed that the buildings are located close to each other so that their cooperation is feasible and valuable. The goal of the cooperation is to gain as much autonomy from the external grid as possible. At each time step, t, the optimization problem is solved, considering the planning horizon. By solving the optimization problem, the optimal charge/discharge schedule of the BESS is calculated. For each time step, t, the first control step of the schedule is executed. This means that each BESS gets charged/discharged with the power that the optimizer defines as optimal for that time, t, taking into consideration the forecast of the planning horizon. The BESS and PV system of each building are modelled separately, as presented in Constraints (2)- (12) and are utilized in order to optimally cover the district's demand. It is noted that all values of the formulas of the proposed methodology are calculated per unit, where the base power is equal to 1 kW. The objective function and the constraints of the optimization problem are presented below: The objective Function (1) minimizes the total interaction between the district and the external grid for the planning horizon. It includes the power imported from the external grid at time t, F t , as well as the power exported to the external grid at time t, T t . The interaction with the external grid is described through these two non-negative variables so that the optimizer aims at total interaction as close to zero as possible: Constraint (2) represents the power balance in the district at time t. P m,t pv is the power produced by the PV of building m, at time t, and may be utilized in order to cover the demand or charge the BESS (of building m or more of the district's buildings), or even be injected to the external grid. D m,t is the power discharged from the BESS of building m, at time t, and may be utilized in order to cover the demand (of building m or more of the district's buildings). C m,t is the power charged to the BESS of building m, at time t, and may originate from the district's PV panels (i.e., panels of building m or even panels of the other buildings). Finally, L m,t is the load of building m, at time t. Overall, Constraint (2) indicates that the total demand of the district may be covered by the external grid, the district's PV panels and the district's BESS. It takes into consideration the power charged to/discharged from the BESS, as well as possible flow towards the external grid (which is possible but against the objective of the optimizer). Since the proposed strategy aims at the optimal synergy, Constraint (2) states that the PV and BESS of each one of the buildings may be utilized for the benefit of the district: Constraint (3) represents the energy balance of the BESS of building m, at time t, where Β m,t is the battery energy of the BESS of building m, at time t and h m is its roundtrip efficiency. It is assumed that the charge and discharge efficiency are equal to each other and their product is equal to the round-trip efficiency, h m . Overall, Constraint (3) models the energy stored in the BESS of building m, at time t, considering the energy stored in this BESS in the previous time step as well as the power charged to/discharged from it: The battery energy of the BESS of building m, at time t, Β m,t , is limited by the technical limits of Constraint (4), where Β m min is the minimum battery energy for the BESS of building m, while Β m max is its maximum. This Constraint models the limitations of battery energy of each BESS, in order not to exceed the upper and lower limits: Constraint (5) initializes the state of charge of the BESS of building m, Β m initial , where t is the current time step and Β m initial is equal to the battery energy of BESS m in the beginning of the current planning horizon. Constraint (6) denotes that the BESS of building m cannot be charged and discharged at the same time. The binary variables u m,t ch and u m,t dch are set to 1 if the BESS of building m is being charged/discharged at time t respectively. This Constraint models the direction of the power flow of each BESS at each time step: Constraint (7) represents the battery energy of the BESS of building m after the end of the last hour t, f m , and Constraint (8) Constraints (9) and (10) Constraints (11) and (12) ensure the parity of the district's BESS. This means that the charge/discharge of each BESS is according to its capacity. For example, if all BESS have the same capacities, they should be charged/discharged equally at each time step. This ensures that all BESS wear out equally: It is noted that the PV production aims to cover as much of the energy demand as possible. When the production exceeds the demand, the exceeding energy charges the BESS. It should be mentioned that if the exceeding production was high enough to violate the technical limits of the BESS, the remaining production would flow towards the external grid.
The optimizer reaches the optimal decisions through the Basic Open-source Nonlinear Mixed INteger (BONMIN) solver.
In relation to the modelling of the PV panels, it is noted that the total radiation reaching the vertical surfaces is calculated with the use of the Perez model [48], taking into consideration their azimuth and angle. Afterwards, the calculated radiation, the ambient temperature and the technical specifications of the panels are used, in order to calculate the PV production of each building m, at time t, which is the input of the optimizer. The procedure described above is carried out through AixLib [49], which provides tools that calculate PV production out of meteorological data. Overall, Aixlib calculates the production of vertical BIPV systems, providing the controller with the necessary data, i.e., the PV production of each building throughout the day, in order to calculate the optimal actions.
Finally, the limitations of the proposed strategy are, (a) the fact that the objective function aims explicitly at the autonomy of the district, not taking into consideration factors such as stability, power quality, etc., which are mostly related to lower (hierarchically) control levels, (b) the time step, which is selected to be equal to 1 h, and (c) the distance between the buildings, which need to be located close to each other so that the losses can be omitted.

Formulation of the Forecaster
When it comes to optimization problems some of the parameters may be unknown. In this case, the unknown parameters are the PV energy production and the load of each building, assuming a time interval of one (1) hour. These parameters are usually forecasted. But no matter how accurate a forecasting method may be, it is possible that the forecasted parameters have significant deviation from the real values. In order to address this issue, the error needs to be detected as soon as it happens so that the forecasted values for the rest of the planning horizon can be adjusted to reality as closely as possible.
The proposed forecaster is based on WMA, utilizing Trigg's tracking signal and adjustment formulas of the forecasted values. The developed method includes sensitivity parameters and thresholds that allow the selection of a stricter or looser approach of the time series that need to be forecasted. The flowchart of the method is presented in Figure  3. The forecasted values in the beginning of the day are calculated using WMA on recent time series, giving more weight on the more recent than the older ones. The mathematical formulas of WMA forecasting method are presented in Equations (13)- (15). Equation (13) is the forecast of random parameter y. For this forecast the recent data x i are taken into account multiplied by their weights, w i , where i is the index of the days included in the history, ranging from i = 1 (the oldest day) to N i (the most recent day). It is essential that the summation of weights is equal to 1, as shown in Equation (14). Parameter x i is considered to be less recent than parameter x i+1 , so its value is less important for the prediction of y than the value of x i+1 . This means that x i should have less weight in the prediction of y than x i+1 , as expressed in Equation (15): w i <w i+1 , ∀ i<N i (15) In this case, the history that is taken into account consists of the past four days. This means that N i = 4 and the selected weights are the following: w 1 = 0.1, w = 0.2, w = 0.3 and w 4 = 0.4 [50,51]. For the detection of significant deviations from the initial forecast, at each hour, t, the real values of the load and PV production are measured at the district and then imported to the forecaster. These values are compared to the forecasted ones using Trigg's tracking signal. This index, Trigg t , detects significant deviations between forecasted and real-time series using the formula shown in Equation (16), where e j is the error of the forecasted value at time step j, j is the index of time steps since the previous adjustment of the forecaster and e 0 is the initial error of Trigg's index. For the selection of the parameters, i.e., a and b, the limitations of (17), (18) and (19) should be taken into account [52]: 0.05≤a≤1 (17) 0.05≤b≤0.5 (18) b≤a (19) In this case the selected parameters are the following: a = 0.2 and b = 0.2 [52]. If Trigg's index is greater than a certain threshold, which in this case is set equal to 0.3 (according to tests), it means that the forecast does not match the real time series. In order to resolve this issue, the forecasted values for the rest of the planning horizon need to be adjusted. The adjustment should take into consideration the random nature of the parameters that need to be forecasted. It should be noted that the production of PV is more predictable than the load of a building, especially when it comes to residential loads. This means that the adjustment formula of a building's PV production should be different than the adjustment formula of its load. In fact, since the PV production is not as random as that of the load, the adjustment of the PV forecast should be greater than the adjustment of the load forecast. The formulas of adjustment are presented in Equations (20) and (21), as follows: Equation (20) represents the adjustment formula regarding the forecast of the PV production, whereas Equation (21) represents the adjustment formula regarding the forecast of the load, where F m,t+1 old is the old forecast and F m,t+1 new is the new forecast. In the case of the PV production, the top limit is taken equal to the PV rated power, PV m max , omitting for simplicity reasons the losses in the PV installation. In the case of the load, the bottom limit is set equal to zero. It is noted that the adjustment takes into consideration the weighted error of the forecasts made by the previous adjustment. Parameter c indicates the sensitivity of the adjustment. In this case it is set to 20%. It is noted that every time that the forecasted series are adjusted, Trigg's index is reset in order to evaluate the performance of the adjusted forecasting formula. This means that and j are reset equal to 1, as they were in the beginning of the day, and e 0 is reset equal to the absolute average value of the past three forecasts.

Case Studies
The days that have been selected for the simulation of the implementation of the MPC-based control on a hypothetical district located in Greece are the 10th of March, the 26th of July, the 7th of October and the 16th of December, assuming more or less a day in each season of the year. The selection of the simulated day of each month is based on the profile of radiation. More specifically (considering vertical placement, facing the South) the measured radiation is close to the daily average values of radiation on vertical façade facing the South for the period 2005-2016, calculated by [53], presented in Figure 4. It is noted that the vertical position of the BIPV panels causes higher absorption of radiation during winter than during summer. The Typical Meteorological Year (TMY) weather data can be found in [54] and the load data can be found in [55]. The performance specifications of each implemented PV panel are presented in Table 1 [56]. The four residential buildings are considered to be equipped with 12 panels each, whereas the school is considered to be equipped with 14 panels. The selected number of PV panels is based on the average wall surface of newly built residences, according to [57], considered as representative of the EU. Also, each building is considered to be equipped with a BESS, whose performance specifications are presented in Table 2 [58]. It is noted that, in the beginning of each simulation, the state of charge of each BESS is considered to be equal to the minimum, which is 0.5 kWh, and represents the worst-case scenario for the building's autonomy. The district is part of the CIGRE benchmark system [59] and is simulated with the use of Modelica™ [60]. The results of the simulations regarding the distribution of the PV production (i.e., the part of total PV production directly feeding the load of the district and the part of total PV production charging the five BESS of the district) of each simulated day are presented in Figure 5 and the values of the total PV production are presented in Table 3. The results regarding the energy mix that covers the total load of the district, including all buildings, are presented in Figure 6 and the total values are included in Table 4. The energy mix may include energy imported from the external grid, PV production feeding directly the district's load and energy discharged from the five BESS of the district.    It is noted that the lowest production occurs in the representative day of July, i.e., 64.85 kWh. This observation has also been discussed by the authors of [61], and can be attributed to the vertical position of the PV panels and the reliance of their production mostly on diffuse radiation.
From the energy mix diagrams of Figure 6, it is observed that, in general, the grid contributes to the district during the first hours of each day. This is the period when the PV panels do not produce energy (because there is no daylight) and the five BESS of the district cannot cover the load, as they are assumed to be in the minimum state of charge. It is also noted that during the middle of the day, which is the time when the PV panels reach their peak production, the load of the five buildings is covered mostly with energy from the PV production. After the PV production is reduced to zero, the demand is initially almost exclusively covered with energy provided by the BESS, which have been charged from the exceeding PV production (during noon), as presented in Figure 5, to be followed by the grid when the BESS reach their minimum limit. Overall, the energy mixes presented in Figure 6 show how the mismatch between the overall PV production profile and the district's load profile is corrected through the utilization of the district's five BESS, controlled by the proposed MPC strategy. In this way, when the overall PV production is higher than the district's demand, the surplus energy optimally charges the five BESS and when the demand exceeds the PV production, the five BESS are optimally discharged, as long as they have sufficient stored energy to provide.
The maximum total PV production occurs in the simulated day of March, i.e., 82.18 kWh, and covers directly higher demand than the rest of the simulated days, i.e., 44.75 kWh, due to the fact that the curve of the PV production matches the load peak in the morning-noon hours of March. Also, in the simulated day of March, the highest amount of energy provided by the BESS occurs, i.e., 33.69 kWh, followed by the respective value for the simulated day of December, i.e., 29.81 kWh. This is attributed to the high demand observed in the evening hours for the simulated days of March and December, as presented in Figure 6, after the PV production has ceased. In general, the load of the simulated days of March, i.e., 96.46 kWh, and December, i.e., 94.97 kWh, appears to be higher than the load of the rest of the simulated days, mainly because of the thermal needs that are covered through electrical means, such as heaters, heat pumps etc. In fact, the load of these days is high enough to also require the contribution of the grid in the evening hours, as presented in Figure 6, calculated as equal to 1.89 kWh for the simulated day of March and 5.51 kWh for the simulated day of December.
For the simulated day of October, it is noted that in spite of the high PV production, i.e., 77.93 kWh, only a small portion of energy feeds the load directly, i.e., 33.78 kWh, while the remaining energy, i.e., 44.15 kWh, is used to charge the BESS. From all simulated days, this is the highest amount of energy headed towards the BESS and the lowest amount of energy feeding directly the load. Both are attributed to the mismatch between the profile of PV production and the load profile, which is shown in the diagrams for the simulated day of October in Figures 5 and 6. This mismatch is covered by importing more energy from the external grid during the morning hours, which decreases the daily energy savings and increases the usable energy stored in the BESS at the end of the day, i.e., 21.34 kWh.
The daily savings are presented in Table 5, utilizing Equations (22)- (24). Equation (22) denotes the total energy savings, E tot , i.e., the total energy directly utilized from PV panels, E PV , and BESS, E BESS , respectively, in order to cover the total load. Equation (23) denotes the percentage of the energy savings, which is calculated as the total energy savings, E tot , divided by the total consumption, L tot . Equation (24) denotes the daily cost saved, CS, due to the implementation of the BIPV and the BESS. It should be noted that the percentage of cost savings is the same as the percentage of energy savings, calculated by (23). The price of the energy, p, is derived from the electricity prices (including taxes) for household consumers in Greece and is set constant throughout the day and equal to 0.1551 €/kWh, according to [62]. It is estimated (taking into account the average daily load of the simulated days, equal to 81.29 kWh) that before the installation of the adaptable/dynamic building envelopes, the average daily energy cost of the entire district is equal to 12.61 €.
CS=(E PV + E BESS )*p (24)  The highest percentage of daily energy savings, i.e., 81%, is observed in the simulated days of March and July. In the case of March, which also has the highest total energy savings, i.e., 78.44 kWh, this is attributed to the high PV production and the fact that the demand profile matches the PV production profile. In the case of July, this is attributed to the low total demand, i.e., 67.12 kWh. On the other side, the lowest percentage of daily energy savings, i.e., 76%, is observed in the simulated day of December, due to the high load. In addition, the lowest daily energy savings are presented in the day of October, i.e., 52.17 kWh, due to the mismatch between the profiles of PV production and load, which causes the grid to cover the demand. As regards the economic aspects, the average daily cost savings of the district, based on the four simulated days, are calculated as equal to 9.98 €. This means that the average percentage of cost savings of the district is equal to 79%.

Sensitivity Analysis
The sensitivity analysis consists of two scenarios. The first scenario, scenario A, assumes that the residential consumers are equipped with 8 PV panels (instead of 12) and the school is equipped with 12 PV panels (instead of 14). The second scenario, scenario B, assumes that the first residential consumer is not equipped with a BESS, which means that there are only 4 BESS instead of 5 in the district. The comparison of the results of the sensitivity analysis with the reference case is included in parentheses in the following respective Tables.
The results of the simulation regarding the distribution of PV production (i.e., the part of total PV production directly feeding the load of the district and the part of total PV production charging the five BESS of the district) of scenario A for each simulated day are presented in Figure 7. The results regarding the energy mix that covers the total load of the district, including all buildings, in the case of scenario A, are presented in Figure 8 and the total values are included and compared to the reference case in Table 6. The energy mix of scenario A may include energy imported from the external grid, PV production feeding directly the load and energy supplied by the five BESS of the district. The daily savings of scenario A are presented in Table 7, utilizing Equations (22)     Overall the profile of the energy mix of each simulated day of scenario A, presented in Figure 8, follows the trend of the respective energy mix of the reference case (presented in Figure 6). This means that in both cases, the main power supply during the early hours of the simulated days is the external grid, due to the lack of PV production and adequate battery energy. Furthermore, in both cases the demand in noon is mostly covered by PV production. Finally, as in the reference case, scenario A, the five BESS (which are charged when the PV production exceeds the demand, as presented in Figure 7) are optimally discharged in the evening, when the PV production is insufficient, as presented in Figure 8.
As presented in Figure 7 and Table 6, the total and the maximum PV production generated by the district's buildings in the case of scenario A is lower than in the reference case for all simulated days, i.e., 29%, due to the fact that there are fewer PV panels. Consequently, the BESS are charged with up to 61% less energy than in the reference case, depending on the simulated day (maximum percentage reduction observed for the simulated day of July), and the PV production is mostly used in order to cover the demand directly, as shown in Figure 7.
In all cases, except for the simulated day of October, the total discharge of the five BESS is lower than in the reference case, ranging from 43% to 60%, as the BESS have received less energy from the PV panels. However, for October, the total discharge of the BESS towards the load, i.e., 19.52 kWh, is slightly higher than in the reference case, i.e., 18.39 kWh. This happens because the PV production in both cases charges the BESS with enough power to cover the load for the rest of the day, but in scenario A, the PV production is lower than in the reference. This causes the BESS to cover the difference, thus, increasing the total energy provided by them.
As presented in Figure 8, only in the case of October, the BESS energy is enough to cover the load for the rest of the day, after the PV energy production is reduced down to zero. In all simulated cases, due to lower PV production, more energy needs to be imported from the grid, i.e., 14% to 120% increase, in order to cover the demand, especially in the cases of March, i.e., 120%, and December, i.e., 87%, due to the high total load. The higher contribution of the grid can be observed in the diagrams of Figure 8, especially after the PV panels stop producing power, when the energy provided by the BESS is in most cases insufficient.
Also, in the case of October, it should be highlighted that the energy stored in the BESS for the next day is significantly lower than in the reference case, 87%, due to the reduction of PV production and the higher cover of load from BESS. For the rest of the simulated days, i.e., the simulated days of March, July and December, the BESS have no energy stored for the next day.
As expected, the savings of scenario A are lower than in the reference case because of the lower production of the PV. The energy savings of the simulated day of March remain the highest, i.e., 56.86 kWh, due to the size of the PV production, i.e., 58.38 kWh. However, the highest percentage of energy savings, i.e., 75%, is observed for the simulated day of October. This is attributed to the high contribution of the BESS, i.e., 19.52 kWh, and to the fact that the load is lower than the rest of the simulated days, i.e., 66.59 kWh. For scenario A, the average daily cost savings are calculated as equal to 7.91€. This means that the average percentage of cost savings of the district is equal to 63%.
When it comes to scenario B, the results of the simulation regarding the PV production are identical to the reference case ones due to the fact that the number of panels remains the same. However, the energy mix that covers the total load is different because there are four BESS instead of five. The results regarding the energy mix, which covers the total load of the district in the case of scenario B, are presented in Figure 9. The values of total energy provided to the district by each BESS are concluded in Table 8.  As expected, scenario B presents the greatest amount of energy provided by each BESS compared to the reference case, as the PV production and the load remain the same and the number of BESS decreases from five to four. This causes the four BESS to manage more energy. Out of all simulated days, the highest amount of energy provided by each BESS occurs in the simulated day of March and is calculated as equal to 8.42 kWh. It is noted that the parity of the BESS, ensured by Constraints (11) and (12), results in the equal increase of the discharge of each remaing BESS by 25%. It should be mentioned that significant usage of BESS wears them down sooner and decreases their lifetime.
The total contribution of the grid remains the same as in the reference scenario for each simulated day, because the PV production remains the same, the same amount of PV production is provided to the four BESS and the same total amount of energy is provided by the four BESS to the district, without challenging their technical limits. It is noted that since the PV production, the total energy provided by BESS and the total energy imported from the grid remain the same, the savings remain the same as in the reference case. Nevertheless, if the BESS were even fewer or if the demand was higher, it is possible that the technical limits of the BESS would be challenged, resulting in lower daily savings.

Conclusions
This paper presents a control methodology for the energy management of a district that consists of buildings that are equipped with vertical BIPV and BESS. The proposed MPC-based method aims to achieve the maximum possible autonomy of the district through the cooperation of the interactive buildings. At each time-step, the method seeks the operation schedule that fulfills the objective of the controller, which is the minimum possible energy exchange between the district and the external grid, taking into consideration the limitations of the model and the forecast of the uncertain parameters for the planning horizon, i.e., the PV production and the load of each building. The operation schedule is applied to the district, which consists of a five-node distribution network, based on the CIGRE distribution benchmark and is renewed at each time-step, according to the controller. Τhe derived results for four representative days of the year show that the cooperation between the buildings has quite satisfying impact, in terms of autonomy from the external grid, with daily energy savings ranging from 76% up to 81%, which verifies the effectiveness of the proposed technology and method. Additionally, two scenarios of sensitivity analysis have been conducted, the first one regarding the number of PV panels and the second one regarding the number of BESS. The results of the sensitivity analysis indicate the effect of changes of the selected parameters on the system behavior during the controller's implementation where the calculated daily energy savings range from 55% up to 81%.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available because they constitute bearing proprietary information.

Conflicts of Interest:
The authors declare no conflict of interest.