**Granger Causality Network Methods for Analyzing Cross-Border Electricity Trading between Greece, Italy, and Bulgaria**

**George P. Papaioannou 1,2,**†**, Christos Dikaiakos 1,3,\*,**†**, Christos Kaskouras 1,**†**, George Evangelidis 4,**† **and Fotios Georgakis 5,**†


Received: 15 January 2020; Accepted: 11 February 2020; Published: 18 February 2020

**Abstract:** Italy, Greece, and, to a lesser degree, Bulgaria have experienced fast growth in their renewable generation capacity (RESc) over the last several years. The consequences of this fact include a decrease in spot wholesale prices in electricity markets and a significant effect on cross border trading (CBT) among neighboring interconnected countries. In this work, we empirically analyzed historical data on fundamental market variables (i.e., spot prices, load, RES generation) as well as CBT data (imports, exports, commercial schedules, net transfer capacities, etc.) on the Greek, Italian, and Bulgarian electricity markets by applying the Granger causality connectivity analysis (GCCA) approach. The aim of this analysis was to detect all possible interactions among the abovementioned variables, focusing in particular on the effects of growing shares of RES generation on the commercial electricity trading among the abovementioned countries for the period 2015–2018. The key findings of this paper are summarized as the following: The RES generation in Italy, for the period examined, drives the spot prices in Greece via commercial schedules. In addition, on average, spot price fluctuations do not affect the commercial schedules of energy trading between Greece and Bulgaria.

**Keywords:** cross border trading; Granger causality; electricity trading; spot prices

#### **1. Introduction**

Over the several past years, European electricity markets have gone through a process of significant developmental changes which contributed to the further liberalization of the energy sector. In addition to this, European countries are committed to reach specific targets for 2030 and 2050 regarding the percentage of renewable energy sources (RES) participation in the energy generation mix, alongside the targets set by the European Green Deal [1]. Variable energy sources (VRES), especially solar and wind, cannot produce constant power since they are highly correlated with weather conditions.

European electricity markets are becoming increasingly integrated under the target of a single European electricity market. Greece has not yet market-coupled with any of the interconnected countries, but it is expected to do so under the "European Target Model" (ETM) mechanism by June 2020. The traded energy volume appears to have a severe degree of seasonality with human activities

slowing down towards summer and then increasing again, based on European electricity market reports [2]. In general, energy trading follows demand. The implemented policies led the traded energy volume within European countries to reach the 33% of the total consumption in the first quarter of 2019 which indicates a 4% increase compared to the corresponding period of the previous year [3]. Nevertheless, in the second quarter of 2019, a reduction of 14% compared to the same period of the previous year in trading volume appeared [4]. The aforementioned policy of increasing the number of market-coupled countries within Europe and hence enforcing the pan-European energy market was initially designed in order to achieve energy price convergence and thus ensuring the highest level of safety and security of supply [5].

On the other hand, European countries following the 2030 climate and energy framework have pledged to reduce the use of fossil fuels and increase RES to improve their share in the energy generation mix [6]. This has led to a remarkable increase of RES, particularly solar and wind power. The European Union has committed to obtaining 20% of the consumed energy until 2020 from RES with this percentage increasing to the 32% in 2030 [7]. These targets will be revaluated for an upward revision in 2023 [6]. By 2017, the respective percentage was 17.5% with some countries having already achieved their individual targets [8].

The increased share of RES has a direct effect on spot prices since, as investigated, it is highly correlated with solar and wind intermittent availability. In general, deployment of renewable (mainly wind and solar) generation capacity affects the dynamics of electricity spot prices in a negative way, while the opposite happens with its volatility [9–11]. This should be anticipated since RES are more efficient than conventional generators in terms of marginal cost, but on the other hand, they cannot guarantee a constant and secure supply, since they are highly dependent on exogenous parameters.

This combination of uncertainty and lower spot prices discourage stakeholders from investing in new capacities, either from renewables or conventional sources. However, market designs can, to some extent, smoothen the volatility of electricity prices and enhance investments [12].

In this paper, we studied the abovementioned issues by considering the interconnection of Greece with Bulgaria and Italy and the effect of one country's intermittent generation to the other's market clearing price and, more importantly, on their cross border trading (CBT). In addition to this, we tried to leverage the available data and detect any hidden (causal as we will mention in the results section) connections among CBT fundamentals. All the examined countries invest in VRES deployment while taking advantage of the interconnection. The level of the system marginal price (SMP) in Greece, which is consistently one of the highest in Europe, discourages new investments in energy capacity. However, the interest is high, since the Greek energy market has room for improvement, especially under the implementation of the ETM. Therefore, we believe that our work will highlight the effect of RES growth in these countries and how this contributes to the changes in Greek SMP and commercial schedules.

Currently, Greece has active electricity interconnections with Turkey, Albania, North Macedonia, Bulgaria, and Italy with the latter being the most active regarding total exchanges. In 2017, 22.72% of energy imports in Greece were from Italy, while the respective percentage of energy exports was 26.82% of the total. The total imports in Greece for 2018 were 11,223.913 MWh, and exports were 4983.061 MWh [11].

However, it should be mentioned that the above countries are connected to each other via physical interconnection and not by market coupling, hence comparison with other studies that investigate countries that are connected via market coupling might lack common ground. However, misfunctions of simple interconnection compared to market coupling can be detected.

The rest of this paper is organized as follows: A thorough literature review is presented in Section 2. Previous works on CBT using multiple models are illustrated in this section. In addition to this, detailed information on the methodology is provided. Section 3 briefly presents the market structure of the examined countries and the interconnection between policy and procedure. Sections 4 and 5 illustrate the data (time series) we used for the completion of the current work on the preprocessing procedure, their summary statistics, and the correlations among them. The Granger causality connectivity analysis

is briefly presented in Section 6. The methodology and the results exploitation as well as the data preprocessing and the model's validation are shown in this section. Finally, the outcomes of the simulations along with a discussion are provided in Sections 7 and 8, respectively.

#### **2. Literature Review**

Various studies have been conducted examining how the increasing penetration of renewables combined with interconnections among neighboring countries affect the dynamics of electricity spot prices domestically and cross border as well as their inter-trading commercial schedules. Since Germany is the leading country in Europe implementing renewable generation in its energy mix and one of the countries with the most cross border interconnections (CBI), it has been the center of research by the scientific community [10,13,14].

The first layer of the investigation consists of the consequences the increasing RES capacity has on the domestic SMP. Ketterer (2014) [13] used daily data on a GARCH model to examine the effect of domestic wind electricity generation on the volatility of electricity price. Specifically, the study showed that increasing wind generation reduces the price levels while at the same time increases its volatility.

In addition to the previous study, Wozabal et al. (2014) [10] using hourly data and ordinary least squares regression (OLS) showed that Intermittent Energy Sources (IES)affects SMP in a complex way. For a small to moderate quantity of IES, the price variance tends to decrease, while large quantities of IES have the opposite effect. Furthermore, in their study they highlighted some policy measures that might support variance absorbing technologies. Grid interconnections were one of the suggested measures, in the sense that they might operate as a means of ensuring the development of sufficient capacities over time.

Extensive investigation has been conducted regarding the Germany–France interconnection. German intermittent electricity generation has been proven to affect the dynamics of French spot prices. Specifically, increasing German renewable generation leads to a decrease in the level of French spot prices but has a positive effect on its volatility [12]. A step forward has been made by Haximusa (2018) [14], who showed that wind and solar generation in Germany have an ambiguous effect on the volatility of French spot prices. He considered three levels of French demand (i.e., low, medium, and high) and applied a GARCH analysis to show the effects on the spot price regarding the different demand levels. He proved that during medium and high French demands, importing energy from Germany decreases the spot price volatility, while the opposite happens when the demand is low.

Another example of strong interconnection coupled with increasing competition is the Nordic electricity market. Bask et al. (2008) [15] applied a stochastic model to generate electricity prices at Nord Pool and then used a GARCH analysis to show that prices became less sensitive to external shocks while the Nordic power market was being expanded and the degree of competition was being increased.

Denny et al. (2010) [16] investigated how the increasing penetration of wind generation affects the spot price dynamics of countries coupled to the UK and Ireland. The results indicate that increasing interconnections will reduce the level as well as the volatility of spot prices in both countries.

A paper that is very related to our work here, as far as its main topic is concerned, is the work of Zugno et al. 2013 [17], dealing with the influence of wind power generation on European cross border power flows. Wind power generation and spot price have been found to have a non-linear effect on cross border power exchange across Europe. Using Principal Component (PC) analysis (to reduce the problem's dimensionality), cross border power exchange hourly data were used as dependent variables in local polynomial regression using the PC (extracted from the matrix of cross border flow variables) as exogenous factors. The main findings in this work were that an increase in forecasted wind power generation causes a fall in the German import of power (or rise in the export), while rising spot prices show the opposite direction. Another very significant finding was that, from a global perspective, variations in wind power generation in Germany had significant effects on power flows in Europe. More specifically, import and export patterns were largely altered and loop flows

were originated. The abovementioned paper was a data-driven research study made possible due to the availability of data on wind power generation, consumption, and power flows provided by the European Network of Transmission System Operators for Electricity (ENTSO-E). In general, the higher the wind power penetration in Germany, the more this country exports to its direct neighbors to the South. Also, it appears that there is a loop flow in the power transit from Germany to Switzerland via Austria, since the flow from Austria to Switzerland is positively correlated with German wind power generation. In our work, presented in this paper, we tried to detect similar interactions using, however, Granger causality analysis, on the CBT between Italy–Greece–Bulgaria.

In the present work, we have applied Granger causality (GC) on hourly data extracted by ENTSO-E database for four years (2015–2018). To the best of our knowledge, the cross-border electricity exchange among the abovementioned countries has not yet been investigated and, hence, our work appears to provide a valuable contribution to the scientific community. Another aspect of our work is that we examined how a country's spot price (e.g., Greece) and commercial (trading) schedule was affected by VRES generation in a neighboring country (e.g., Italy) without the existence of a market coupling policy.

Granger causality, an established time series technique, was used to analyze the interactions ("cause and effects") of a representative sample of 13 European electricity spot prices for the period 2007–2012 [18]. The study applied GC via network theory and provided inferences regarding the European electricity network's state and dynamic evolution over time, therefore assessing spot price convergence and "quality" of market coupling. The causal interactions among European electricity spot prices were modeled as a connectivity network on which the spot prices constituted the nodes of the network, while the links corresponded to the significant influences among relative pair-wise price changes.

One of the early applications of causal flow modeling in electricity markets is the work by Park et al. (2006) [19], in which the authors used advanced techniques in causal analysis (vector autoregressive (VAR), vector error correction model (VECM), and directed acyclic graphs (DAGs)) to find the dynamic relationships between electricity spot prices and the prices of major electricity-generation fuel sources (i.e., oil, gas, and coal prices) in US electricity spot markets.

A multivariate version of GC was used in a paper by Narayan et al. (2009) [20] to examine the causal relationships between electricity consumption, exports and gross domestic product for a panel of Middle Eastern countries.

In their causal modeling and inference for electricity markets, Ferkingstad et al. (2011) [21] used a combination of VAR, VECM, and a linear non-Gaussian acyclic model (LiNGAM) which they call time-lagged causal flow, a concept very close to the GC. They applied their hybrid model to weekly Nordic and German electricity prices, using oil, gas, and coal prices with German wind power and Nordic water reservoir levels as exogenous. They showed that, in contemporaneous time, Nordic and German spot prices were interlinked through gas prices.

A combination of out-of-sample Granger causality tests and DAGs were also used by Yang and Zhao (2014) [22] to investigate the temporal linkages among economic growth, energy consumption, and carbon emissions in India.

The methodology applied in this work, referred onwards as Granger causality connectivity analysis (GCCA), is based on modern network theory, an efficient approach to characterizing connecting systems [23]. Our purpose was to study the dynamic evolution of the network and the nodes which included the following criteria of three countries: the spot prices; imports–exports, the demand (load), the VRES generation; and the commercial schedules of CBT. This study will inform how the spot prices of the above countries are influenced by their own interactions as well as by the interactions of the other variables. Based on the results of this study, we will be able to draw conclusions regarding the cross-border trading development of this "peripheral" sub-system of the European electricity system.

The identification of directed operation connectivity among the components of a system is a challenging work. The GCCA is a powerful tool in detecting such connectivities in complex systems such as an electricity market. Granger Causality is a way to investigate "causality" between two variables, a probabilistic account of causality, and a concept closely related to the idea of cause and effect but not the same. GC allows one to know which particular variable comes before another in a time series but does not describe a causal link in the true sense. In Econometrics, "cause" is actually realized as "Granger cause".

The most recent approach to complex systems is the contemporary network theory, originating from the small world theory of Watt and Strogatz (1988) [24], as well as its scale-free "version" of Barabasi and Albert (1999) [25]. Recently, a large number of papers has focused on various applications of complex networks [26]. We applied this fascinating approach to the CBET among three neighboring countries to extract valuable information regarding their interactions. The results are expected to be useful in policy design among other areas. All calculations in this paper were based on the theoretical approach to GCCA presented in the work of Seth (2008, 2010 2014) [27–29] and implemented by using his toolbox run on MATLAB (ver.2019a, The MathWorks Inc., Natick, MA, USA).

A review paper presenting the application of GC in energy economics research is given by Narayan and Sangth (2014) [20]. The examination of causal relationships between the electricity consumption and economic growth using linear and non-linear GC tests in the case of Turkey is presented in the paper by Nazlioglou et al. (2014) [30]. The tests reveal a bi-directional GC in both the short and long run between the electricity consumption and economic growth in Turkey. The relationship between economic output and energy use is also the focus of a work by Brun et al. (2014) [31]. Woo et al. (2006) [32] applied a Granger instantaneous causality test in order to examine the potential causal relationships between wholesale electricity and natural gas prices in California, revealing bi-directional relationships between these two markets [32]. An in-depth analysis of the causal relationships in oil markets globally is given in a thesis work by Antoniadou (2015) [33]. In this work, the GC and the Toda–Yamamoto approaches were applied, revealing that the causality between crude oil prices and natural gas prices is not bi-directional, since only crude oil prices Granger causes natural gas prices. The relationship between investor attention and crude oil prices is examined by the paper of Li et al. (2019) [34]. Using the Google search volume index (GSVI) that "captures" investor attention, the author used linear and non-linear GC tests. The results show that a bi-directional GC exists only between WTI future crude oil returns and investor attention.

The latent volatility GC for four renewables energy exchanged traded funds (ETFs) and crude oil ETF (USO) were examined in a work by Chang et al. [35]. The empirical results showed that there are significant positive latent volatility GC relationships between solar, wind, nuclear, and crude oil ETFs as well as significant volatility spillovers shocks for renewable energy ETFs. Using a GC in quantiles analysis (evaluating causal relations in each quantile of the distribution), they succeeded in discriminating between causality affecting the median and the tails of the conditional distribution and provided evidence for the existence of a bi-directional causality between changes in RES consumption and economic growth using RES consumption, oil prices, and economic activity data in the US (July 1989 to July 2016).

The GC test was also used in analyzing spillover effect of oil and natural gas prices between emerging and developed countries in a paper by Zhong et al. (2019) [36]. The main findings were that oil and natural gas markets have significant GC and that emerging markets have a strong impact on many developed markets regarding returns and volatility spillover systems.

A combination of GC test and the (recently developed) cross-quantilogram was used on crude oil, natural gas, heating oil, electricity, and gasoline market data to evaluate their directional predictability in a work by Scarcioffolo et al. (2019) [37]. They found no strong evidence to support the decoupling between crude oil and natural gas markets. Natural gas and heating oil were found to be strongly linked across all quantiles.

#### **3. Wholesale Electricity Markets: The Case of Greek, Italian, Iberian, French, and Bulgarian Markets**

#### *3.1. The Greek Market*

Greece's liberalized electricity market was established according to European Directive 96/92/EC. The Greek Wholesale Electricity Market (GEM) currently operates as a day ahead mandatory pool which is based on an optimization algorithm for both energy and ancillary services [38] taking into account the following:


Greek Wholesale Energy Market GEM is in a transitional period towards its final design, namely, ETM where both an intra-day market and a balancing market are expected to operate in order to enhance GEM's liquidity and efficiency. More specifically, the new regulatory framework for the ETM makes provisions for a day-ahead and intra-day market managed by Hellenic Energy Exchange S.A. (Athens, Greece) and a balancing market managed by Greece's Transmission System Operator (IPTO) (Athens, Greece). The new day-ahead market will be a semi-compulsory market, where orders should cover the availability and should be compatible with Price Coupling of Regions algorithm (PCR EUPHEMIA) standards, and the biddings will be on a physical asset basis except RES which could be on a portfolio basis. There are also provisions for exchange-based futures and over-the-counter (OTC) contract limits on the volumes.

#### *3.2. The Italian Market*

In Italy, Gestore del Mercato Elettrico (GME) is responsible for the operation of the power market as well as gas and environmental ones. Gestore del Mercato Elettrico runs under the framework of the Italian Regulatory Authority (ARERA) which sets the rules and activities of the abovementioned sectors. The Italian Wholesale Electricity Market essentially runs on the Italian Power Exchange (IPEX), where producers and suppliers trade blocks of energy. More specifically, GME operates a day-ahead market (MGP) based on auctions, an intra-day auction market (MI), a forward market (MTE), and a market (MPEG) for continuous trading of energy daily products. Moreover, it operates on behalf of TERNA which is the Italian Transmission System Operator the Ancillary Services Market (MSD) while at the same time issues the OTC transactions (PCE).

Italy is divided into six zonal markets based on geographical criteria. These are Northern Italy (NORD), Central Northern Italy (CNOR), Central Southern Italy (CSUD), Southern Italy (SUD), the Sicily (SICI), and Sardinia (SARD). Greece is connected to SUD, which is the South part of the Italian electricity system. The zonal price, in the case of no congestion, is unique and is the interception of the aggregated demand and supply curves which are calculated after an algorithmic procedure. In the case of congestion, the market is split into the four aforementioned regions, and the system calculates different price equilibriums. The national unique price (PUN), which is the consumer price, is the weighted average of the individual zonal prices [38].

#### *3.3. The Iberian and French Electricity Markets*

Italy has already established a market coupling relationship with France and hence with the Iberian wholesale electricity market. The Iberian wholesale electricity market MIBEL, is a joint wholesale electricity market which comprises Spain and Portugal [39]. OMI-Polo Español S.A. (OMIE) belongs to the Iberian Market Operator business group, and it is subject to the rules and regulatory framework governing Spain's electricity sector [40]. The Iberian Wholesale Electricity Market is

composed of an intra-day market and a day-ahead market producing a common spot price for both Spain and Portugal, unless there is a violation in the interconnection capacity between them; in such a case, a market splitting mechanism is activated which results in different prices for both countries. The Spanish Transmission System Operator (REE) and the Portuguese one (REN) are responsible for the technical implementation of the daily schedules. Following the gate closure, six intra-day market sessions are held four hours ahead of the physical delivery [41].

The French wholesale electricity market plays a key role in the French power system by allowing the balance between supply and demand. Most of the electricity (95%) injected into the French power system comes from nuclear and other sources (hydro, gas, coal, renewable energy sources) with a small percentage (5%) covered from imports [42]. Electricity products of the French Wholesale Electricity Market are traded on power exchanges or over the counter (OTC) via brokers or through bilateral agreements between the involved parties. The spot products are daily (day-ahead) or weekend products delivered at base load hours or during peak hours (from 8 a.m. to 8 p.m. for the working days) [43]. The abovementioned products are available on an hour or half-hour basis or as complex blocks covering several hours. The spot price of the French Wholesale Electricity Market is the price of the day-ahead market which is run on the European Power Exchange (EPEX SPOT) [44].

France is connected to the Central European System through its six interconnection lines to England, Belgium, Germany, Switzerland, Italy, and Spain. Currently the interconnection capacity between Spain and France is approximately 5% of its installed capacity which is far below the EU targets for interconnections; this limit has been set at 15% for 2030 at the EU level [45]. At the moment, the planned interconnections between France and Spain through the Gulf of Bizkaia is under consultation; this project has been characterized as a project of common interest with an interconnection capacity of 5000 MW. The latter is expected to be in operation between 2024–2025 [46,47].

#### *3.4. The Bulgarian Market*

The Bulgarian electricity market operates under the command of the Independent Bulgarian Energy Exchange (IBEX) which was established in 2014 as a 100% subsidiary of Bulgarian Energy Holding and is among the last countries introducing such an exchange market. A major change that was introduced by IBEX was the establishment of an organized day-ahead market, in 2016, to replace the pre-existing model.

Currently, IBEX operates three trading platforms: a day-ahead market, and intra-day market, and a centralized market for bilateral contracts; it is mandatory for generators with an installed capacity of 1 MW or more to sell their electricity through IBEX.

The electricity market consists of two segments.

• Regulated Market

Prices are set by the regulator (Energy and Water Regulatory Commission) and consumers are supplied based on territory. This segment includes households and small businesses connected to the low voltage distribution network.

• Free (Liberalized) Market

Electricity is freely negotiated and bought by suppliers and consumers directly from the electricity generators or via IBEX. According to the last amendment (May 2019), all RES and co-generation power plants with an installed capacity of 1 MW or more have the obligation to sell their electricity only via the power exchange.

#### *3.5. Markets Comparison*

One of the outmost criteria to measure a market's concentration and hence define its competitiveness is the Herfindahl–Hirschman index (HHI). Even though in various cases it fails to consider the complexities of various markets, it remains a reliable index. Specifically, regarding the EU energy markets, the European Commission has set the limits for the HHI as shown in Table 1 [46].


**Table 1.** Herfindahl–Hirschman index (HHI) in the power generation ''market" of the analyzed electricity markets for 2015-2018 as depicted in the national reports of regulators sent to the European Commission (EC).

\* Legend (according to the European Commission 2014 [46]): vhc = very high concentration (HHI > 5000), hc = high concentration (1800 < HHI < 5000), mc= moderate concentration (1000< HHI<1800), d=deconcentration (HHI < 1000).

Figure 1 illustrates the energy generation mix in 2018 for the examined countries.

**Figure 1.** Energy generation mix for Greece, Italy, and Bulgaria (left to right) for 2018 [47–49] (outcome of IPTO's elaboration).

From the energy generation mix above, we can observe that Bulgaria was primarily dependent on nuclear power (36%) as a source of energy, which does not appear in either the Greek generation mix nor in the Italian, and on lignite (40%). On the other hand, Italy devoted 45% of its generation mix to natural gas, and the corresponding amounts for Greece and Bulgaria were 28% and 4%, respectively.

The common element that is being illustrated is the share of hydro in the examined countries' energy mix. Specifically, hydro in Greece accounts for 18%, in Italy for 17%, and Bulgaria for 13% of the total energy generation. Finally, Greece appears to have the highest percentage of RES generation in the energy mix (31%), followed by Italy (23%) and Bulgaria (6%).

#### *3.6. Cross Border Trade in Electricity*

The future "architecture" of the European CBT will be based on capacity allocation and congestion management (CACM). According to CACM, all TSOs are required to develop deliverables towards implementing the markets coupling, the so-called single intra-day (SIDC) and single day-ahead coupling (SDAC) [50]. The regulation also outlines the methods for calculating how much capacity the participants in the market can use on cross border lines without putting the system at risk. Also, the document harmonizes cross border market operations in Europe in order to enhance the competition and the integration of renewables.

The SIDC creates a single EU cross-zonal intra-day electricity market, in which buyers and sellers of energy cooperate to trade electricity continuously on the day the energy is needed. The SIDC enhances the efficiency of intra-day trading across Europe by competition promotion, liquidity increase (facilitating the mechanism of buying and selling without affecting the price), making the sharing of energy generation resources easier [51].

The SIDC is implemented in three phases or waves. In the first wave, 14 countries went live in June 2018, in the second wave 7 countries, including Bulgaria, went live in November 2019. Italy and Greece will be included in the third wave, foreseen for Quarter 4 of 2020 [52].

The SDAC is the pan-European single day-ahead coupling serving 27 countries (at the time of the writing of this paper). Under this agreement, 33 TSOs and 16 nominated market electricity operators (NEMOs) will cooperate. The SDAC uses (PCR EUPHEMIA), that calculates, across Europe, electricity spot prices and does not allocate implicitly the auction-based cross border capacity [53]. According to the Agency for the Cooperation of Energy Regulators (ACER), SDAC has improved the level of efficiency regarding the usage of the interconnections from 60% in 2010 to 87% in 2018 [54].

On the 1st of October 2018 (first delivery date of the projects related to SDAC), the Germany–Austria bidding zone split was successfully implemented. In South East Europe (SEE), the target times for the bidding zone borders adhering to CACM are as follows [52]:


The electricity trade can be unidirectional or bidirectional over long time horizons. Electricity flows in either direction due to the demand changes in both countries, for example, when the latter is not perfectly correlated due to the different seasonal or diurnal patterns.

The correlation coefficient between Greece's and Italy's forecasted load was 0.524 and statistically significant (so the loads were not perfectly correlated, indicating the existence of strong bidirectionality of CBT).

According to Antweiler [55], higher correlations diminish trade, a very essential finding, since if demand is strongly correlated between two countries, they will both have high and low demand simultaneously, allowing a limited space for additional trade. Another very significant finding in Reference [55] is that the intensity of cross-border electricity trading increases along with the coefficient of variation (cost/total demand) and that the difference or dissimilarity in the size of the countries (e.g., Italy versus Greece) encourages more trade.

A top priority for the European Commission is the harmonization and integration of national electricity markets to a single pan-European market [5]. Energy policy, however, remains strongly a tool of national sovereignty which means that a greater level of integration corresponds to the fact that unilateral national policies can impact interconnected markets.

In this work, we investigated the impact of interconnections on the expansion of renewables promoted by fixed-in tariffs and (unlimited) priority on DA-prices between Italy, Greece, and Bulgaria, as they are revealed through CBT. The contribution of the specific hub-selection lies on the lack of market coupling among the examined countries, and hence the CBT operates as a simple interconnection. It would be of high importance to identify any differences between market-coupled countries' CBT and the specific work.

#### *3.7. Greece's Imports and Exports*

Greece relies heavily on interconnections to meet the demand. This is due to the fact of high wholesale electricity prices which remain high despite the penetration of RES into the system. Moreover, due to the economic recession investment in RES (which, in general, lowers spot prices) had been discouraged.

The average relative margin available for cross zonal trade (MACZT) between 2016 and 2018, as calculated by ACER [54], for Greece and Italy, relative to maximum admissible active power flow (Fmax), set to at least 70% by the clean energy package (CEP), is almost 100%, and the percentage of hours when MACZT is at least 70% for the GR–IT border is also almost 100%.

According to the same report [54], for the border between Bulgaria and Greece, the average DA price differential (€/MWh) and average of absolute DA price differential (€/MWh) during the period 2016–2018 is as follows (Table 2).


**Table 2.** Average and average absolute DA price differential for the Bulgaria–Greece border.

In 2018, the annual average DA prices in Greece (60.4 €/MWh) and Italy (62.04 €/MWh) were among the highest in European bidding zones, whereas Bulgaria's DA prices were among the lowest (39.89 €/MWh). This difference in DA prices justifies the fact that during this specific period, Greece was a major exporter, as shown in Figure 2, in all neighboring countries except Italy.

**Figure 2.** Greece's energy imports and exports for 2018. The net position is marked with red (outcome of IPTO's elaboration).

Figure 3 illustrates the monthly amount of energy trading between Greece and all the interconnected countries for the examined period (2015–2018).

From Figure 3 we observed a reduction in the difference between imports and exports (interconnection balance). It was continuously positive, which means that Greece exported less energy than it imported (it is a net importer) during the whole period between 2015 and 2018.

As far as the rate of implementation of the ETM for DA markets is concerned, significant progress has been made in recent years. The ETM foresees a single DA coupling, enabling the efficient use of cross-zonal capacity in the "right economic flow (direction)" (from low to high price), when there is a price differential across a bidding zone border. The level of efficient use of interconnectors in the DA market timeframe and the estimated social welfare gains still to be obtained, from extending DA market coupling per border, are the two indicators that illustrate the progress toward ETM.

The estimated social welfare gains, still to be obtained, from further extending DA market coupling in the borders of Greece–Italy and Greece–Bulgaria, during the period between 2017 and 2018, as calculated by ACER [56] are illustrated in Table 3.

**Figure 3.** Energy balance between Greece and interconnected countries.

**Table 3.** The estimated social welfare gains from further extending Day Ahead (DA) market coupling in Greece and Italian/Bulgarian borders.


Historically, Greek Wholesale Electricity Prices (SMPs) were driven by the low prices of fossil fuels and more specifically from the lignite ones. The Green Deal alongside climate change awareness made the EU impose stricter rules for CO2 emissions. The 2030 and 2050 emission targets (90% reduction compared to 1990 EU levels) led to higher prices for CO2 emission rights which, in turn, had a dramatic impact on lignite prices. The latter seriously affected the GEM, since gas prices are now compared to lignite prices which, in turn, incurs higher wholesale electricity prices according to the merit order mechanism. Also, even though Greece has already achieved the EU targets for RES penetration levels by 2020, this is not mirrored in market prices due mostly to the financial crisis of 2011 but also to the low volatility and liquidity of GEM.

#### 3.7.1. Joint Allocation Office (JAO)

The Joint Allocation Office (JAO) is a joint service company that consists of 22 TSOs in 19 countries and facilitates the electricity market by organizing auctions for cross border transmission capacity. It was created in June 2015 and was the result of the merge between the Central Allocation Office (CAO) and ASC.EU S.A. located in Germany and Luxembourg respectively, the two former regional allocation offices for cross border transmission capacities. On the 1st of October 2018, JAO became the Single Allocation Platform (SAP) of forward capacity for all European TSOs, operating in accordance with EU legal rules. Currently, JAO performs long- and short-term auctions of transmission capacity and is the duty of the national regulatory authorities to decide which auction should be performed.

#### 3.7.2. Greece–Italy (GR–IT) Interconnection

Currently, Greece and Italy operate as two separate electricity markets, while no market coupling has been applied.

Before the annual and monthly capacity allocation process, IPTO (ADMIE) calculates the GR–IT Net Transmission Capacity (NTC) and agrees on it with TERNA. The final NTC values are stored in the system and then communicated via email to JAO. After receiving the NTC values, JAO is responsible for the annual and monthly transmission rights allocation on GR–IT border. Joint Allocation Office (JAO) performs the auction and informs IPTO and TERNA of the capacity auction outcome.

Two days before the scheduling day (until 12:30 CET), IPTO receives from JAO the annual and monthly final capacity right documents which are aggregated into the long-term capacity rights.

The values of daily ATC and relative data are submitted by ADMIE and TERNA to JAO every day (until 09:30 CET). Once the daily auctions are calculated, they are sent by JAO to ADMIE and TERNA.

After having received the capacity right document, IPTO interacts with the capacity holders or their counterparties who inform IPTO of the capacity they are going to use, based on their capacity rights documents (long term or short term).

The day-ahead scheduling process, performed by IPTO, includes two phases of nominations, submission and matching.

The long-term (LT) phase: until 08:30 CET, capacity holders submit their nominations based on long-term rights. Consequently, IPTO and TERNA exchange and match the LT nominations until 09:00 CET. In the case of a mismatch, the export value prevail rule is applied.

The remaining capacity (non-nominated LT capacity rights), until the NTC of the GR–IT interconnection, comprises the offered capacity for daily auctions (performed by JAO). The IPTO sends both LT and ST capacity rights to ENEX which is the NEMO for the Greek DAM. The SMP of the Greek DAM is the hourly solution of the algorithm of the daily co-optimization of the energy production bids (taking into account the techno-economic constraints of the plants) as well as the energy bids of importers and exporters and the Load declaration of load representatives. The total commercial schedules are determined up to 99% by the Greek DAM clearing, and they have a causal relationship with the price difference between Greek DAM and South Italy marginal prices.

The short-term phase: until 14:30 CET, capacity holders submit their nominations based on both LT and ST rights. Both IPTO and TERNA exchange and match the LT and ST nominations until 15:35 CET. In the case of a mismatch, the minimum value prevail rule is applied.

#### 3.7.3. Greece–Bulgaria (GR–BG) Interconnection

The interconnection between Greece and Bulgaria operates more or less similarly to the GR–IT one, with no market coupling and through the JAO platform. The responsible party for the Bulgarian side is ESO (the Bulgarian system operator).

The auctions between Greece and Bulgaria are bidirectional and are separated to:


The NTCs, which are the auctioning basis, are agreed to between IPTO and ESO. Afterwards, the NTCs and ATCs for each interconnection and direction are published on the auction websites. If there is a change to the ATCs regarding the monthly auctions, users can be informed timely regarding the relevant auction website.

#### **4. Data and Preprocessing**

We focused on the Italian, Greek, and Bulgarian electricity markets, their actual and forecasted solar and wind power productions, demands (loads), DA-prices (spot or wholesale prices), and finally their commercial schedules and transfer capacities of cross-border trading with their interconnected countries. Below, we may refer to them as cross-border trading fundamentals (CBTFs). Table 4 below contains the description of the data which consists of hourly observations, converted to daily values (average of 24 hourly values), covering the period 2015 to 2018 (1460 data points).


#### **Table 4.** Data Overview.

<sup>1</sup> Total commercial schedule is the summation of all the agreed transactions for the delivery and receipt of power and energy among the traders in the interconnection. That schedules are also agreed to by the TSO of the two control areas. <sup>2</sup> Transfer capacity is the allocated through auctions transfer ability to traders that allows them to exchange electricity among geographic areas for each market time unit and for a given direction.

Regarding Bulgaria, there were no available data for the DA price (D21) before 2017 which meant that in the examined cases (study A and C) of this work, which included the specific variable, we limited the data range to the years 2017 to 2018 (725 data points).

In general, we used data for the South Italy region (SUD) which connects to Greece. However, when it came to Load Forecasted, we decided that it was more accurate to use the Total Load Forecasted and not the South one. This decision was based on the procedure IPEX follows to determine the PUN, in which the Load is considered as aggregated for the whole Italy and is not split into zonal regions. Moreover, the total load forecasted in Italy is highly correlated to the load forecasted in South Italy (0.70), so we can safely use the total demand.

Hence, in the forecasted load time series for both Greece and Italy (D7 and D9), we considered in all our calculations the total forecasted load in the Italian and Greek markets.

#### *Testing Examined Time Series for Stationarity*

Since the primary precondition for Granger causality analysis is that the variables must be covariance stationary (CS), also known as weak or wide-sense stationarity, we tested our data for stationarity. Augmented Dickey–Fuller (ADF) and Kwiatkowski–Phillips–Schmidt–Shin (KPSS) tests were performed for testing the null hypothesis that the examined series had unit roots which indicates stationarity or trend stationarity.

Since we used two sets of data (2015 to 2018 and 2017 to 2018 to incorporate the Bulgarian data), the tests were performed for each one of the time series and for both data sets. Table 5 illustrates the tests results for the acceptance or not of the null hypothesis. Rejection of the null hypothesis indicates data stationarity, while the inability to reject it indicates non-stationarity.

The ADF unit root test did not reject the null hypothesis for D24 for the whole series and for D9, D15, D16, and D24 for the reduced data series scenario. On the other hand, KPSS test indicates that D11 and D12 were not trend stationary when the whole series scenario was examined and D21 when the reduced series scenario was examined.

Therefore, for the abovementioned time series, we took their first difference to make them stationary. Then, we performed the tests again and all series appeared to be stationary and trend stationary.


**Table 5.** Augmented Dickey-Fuller (ADF) and Kwiatkowski–Phillips–Schmidt–Shin (KPSS) stationarity test results.

To continue further with our research, we defined three studies (models) in order to find the best models to test the Granger causality of cross-border trading among the examined countries. Our aim was to conduct an analysis for all the possible combinations, hence we separated the model structure as follows:


Table 6 shows all the variables (or nodes in the network) and the corresponding studies (models) they were included as an input.


**Table 6.** The case studies and the variables that were included 1.

<sup>1</sup> The numbers below each time series indicate the corresponding node of the specific variable to the Granger causality connectivity analysis model.

#### **5. Summary Statistics**

#### *5.1. Summary Statistics*

Since Granger causality Analysis is a multivariate (MVAR) method and since all MVAR methods rely heavily on the assumption of normality or near-normality which is often hard to achieve in practice, we proceed in presenting summary statistics of all variables, focusing on normality aspects. The results are illustrated in Table 7.


**Table 7.** Summary statistics.

<sup>1</sup> We converted Bulgarian lev to euros taking an exchange rate of 1 euro = 0.511 lev.

As illustrated in Table 7, none of the examined time series follow normal distribution. Some variables that worth highlighting regarding their skewness values are D2, D12, D14, D21, and D23 whose value appears to be higher than 1, which indicates strong positive skewness and hence data concentration to the right tail (maximum) of the distribution. This was anticipated for the wind generation (D2), since the specific variable was highly stochastic. However, regarding DA prices in Greece and Bulgaria (D12 and D21, respectively), the high value of skewness revealed a difficulty in the forecasting model construction. The same happens with the commercial schedules from Greece to Italy and to Bulgaria (D14 and D23). The contribution of this work is to provide a way of enhancing the forecasting by revealing and indicating any causal connections among CBTF.

Furthermore, from Table 7 we observe that all the variables appeared to have values of kurtosis close to three which was our benchmark of normality, except for D11, D12, D21, and D23. This indicates that the degree of uncertainty in the specific time series (DA prices) was high which contributes more to the difficulty of prediction.

The turnover ("liquidity") of the South Italy electricity market was much smaller than the turnover of the Greek Market. Thus, prices of South Italy zone are more volatile than those at Greece. As illustrated in Table 7 the maximum price at South Italy was 127.30 €/MWh, and the standard deviation was 11.5 €/MWh (mean value 49.72 €/MWh). In Greece, the maximum price was111.4 €/MWh, and standard deviation was 10.63 €/MWh (mean value 52.41 €/MWh).

Finally, observing the Jarque–Bera (JB) test results, we see that none of the examined series were normally distributed, which means that we need to further analyze the data and apply decomposing methods. It is worth mentioning the corresponding value of the JB test for the Bulgarian DA price and the commercial schedule from Greece to Bulgaria (D21 and D23, respectively) which were extremely high, revealing a deviation from the normal distribution.

#### *5.2. Correlation Coe*ffi*cients*

In order to have a "rough insight" into how the variables were interrelated or interacted with each other, we observed their correlation coefficients matrix. However, GC is actually a much more powerful tool to reveal any hidden causality, both in magnitude and direction.

The following tables (Table 8, Table 9) depict the correlation coefficient among the examined time series. We separated the results according to the models we tested to highlight each interconnection and hence make clearer their connection with the results.


**Table 8.** Correlation Matrix: Greece–Italy interconnection (2015–2018).

**Table 9.** Correlation matrix: Greece–Bulgaria interconnection (2017–2018).


Table 8 shows the correlation coefficients of the variables regarding the Greece and Italy interconnection for the period 2015–2018. As illustrated from the coefficients among D1, D2, D3, and D4, there was a relatively high correlation (either positive or negative) that reveals a similarity (similar diurnal effects) in weather conditions which was anticipated between Greece and South Italy. Specifically, there was a high positive correlation between solar (D1 and D3) and wind (D2 and D4) generation in two countries. On the other hand, correlation among solar and wind generation appeared to be negative which reveals that high solar was not followed by high wind generation and vice versa. We did not expect to observe this relation to the results of the current work, since solar and wind generation are exogenous and highly stochastic to be interpreted in the results.

The Pearson's correlation coefficient between PV generation (D1) and spot price (D11) in South Italy was not significant (0.029), while for Greece it was even smaller (0.013). This result is in compliance with results in a similar study conducted by Mayer and Luther [57] on the correlation between PV generation and spot prices in the European Power Exchange (EEX) and Amsterdam Power Exchange (APX) electricity markets. However, these weak correlational linkages between PV generation and spot prices do not "constitute" Granger Causal interactions as shown below.

Another insight from Table 8 is the strong positive correlation (0.526) between total load forecasted of the two countries (D7 and D9). This should be anticipated since the weather conditions, which are the main electricity load driver, were similar. The DA prices (D11 and D12) were also positively correlated with a coefficient of 0.302, something that can also be explained by the weather conditions similarities, since when RES generation is high, the wholesale price is or shows a tendency to be reduced.

Furthermore, a relatively moderate positive correlation (0.206) appears between the commercial schedule from Italy to Greece (D13) and wind generation in Italy (D2). However, we did not observe a similar correlation between commercial schedule from Greece to Italy (D14) and wind generation in Greece (D4). We anticipated an outcome that would partially explain this difference in the specific coefficients.

Table 9 depicts the correlation coefficients of the examined variables for the interconnection of Greece and Bulgaria. Now the examined period was 2017–2018, since as mentioned, we could not obtain previous reliable data for D21.

In contrast to the case of Greece and Italy, Table 9 did not reveal any correlation on weather conditions between Bulgaria and Greece. However, the wind and solar generation in the same country appeared to have a relatively strong negative correlation, −0.280 and −0.396 for Greece and Bulgaria, respectively.

Despite the fact of non-similar weather conditions, the total load forecasted for the two countries had a strong positive correlation of 0.626. A high positive correlation also appeared between D9 and D21, hence it would be interesting to see if a relation between Greek load and Bulgarian DA price does really exist and if the one Granger causes the other. When the demand in Greece was high, the Greek spot price also became high with a consequence to examine import opportunities of cheaper energy from Bulgaria.

A similarity with the case of Greece–Italy interconnection is shown in Table 9 with respect to RES generation and commercial schedules. Specifically, solar generation in Bulgaria (D17) was strongly correlated with commercial schedules from Bulgaria to Greece (D22) (correlation coefficient 0.426). Again, we did not observe a bidirectional relation which means that high RES generation in Greece does not imply high commercial schedules to the opposite direction.

Finally, Table 9 shows a strong relation between the Bulgarian load forecasted (D24) and the domestic DA price (D21), as expected, with a correlation coefficient of 0.304.

#### **6. Granger Causality Connectivity Analysis (GCCA)**

#### *6.1. Granger Causality Test*

The purpose of this paper was to investigate the relationships among the variables related to CBT between two countries, in the same multivariate framework, in which the well-known GC test (based on the MVAR approach) is the most widely used technique.

However, the traditional (in sample) GC test was not adequate to explore the important contemporaneous causal pattern among the variables which is used to conduct a data-determined structural decomposition of the VAR model. As it well known, the conventional VAR analysis relies heavily on the Cholesky decomposition in order to achieve a just-identified system in contemporaneous time. This decomposition is severely criticized for imposing a somehow not so realistic assumption of a recursive contemporaneous causal structure [58]. In order, therefore, to explore the important contemporaneous causal pattern, it is necessary to apply also the DAG tool [59,60].

The advantage of GCCA adopted in this paper, however, is that it combines the above two separate tools into one, "producing" a "new tool", providing similar as well as enhanced results. For information about the typical Granger causal connectivity computational technique as well as for a more theoretical approach to causal networks, the reader is referred to References [27–29].

In the current study, we applied the GC test to examine whether there is an overflow relation among the examined variables. According to GC, if a time series *x* "Granger-causes" time series *y*, then past values of *x* should contain information that helps predict *y* above and beyond the information contained in past values of *y* alone [61]. In one binary *p*-order VAR model.

$$
\begin{pmatrix} y\_{1} \\ x\_{0} \end{pmatrix} - \begin{pmatrix} \varrho\_{10} \\ \varrho\_{20} \end{pmatrix} + \begin{pmatrix} \varrho\_{11}^{(1)} & \varrho\_{12}^{(1)} \\ \varrho\_{21}^{(1)} & \varrho\_{22}^{(1)} \end{pmatrix} \begin{pmatrix} y\_{t-1} \\ x\_{t-1} \end{pmatrix} + \begin{pmatrix} \varrho\_{11}^{(2)} & \varrho\_{12}^{(2)} \\ \varrho\_{21}^{(2)} & \varrho\_{22}^{(2)} \end{pmatrix} \begin{pmatrix} y\_{t-2} \\ x\_{t-2} \end{pmatrix} + \dots \\ + \begin{pmatrix} \varrho\_{11}^{(p)} & \varrho\_{12}^{(p)} \\ \varrho\_{21}^{(p)} & \varrho\_{22}^{(p)} \end{pmatrix} \begin{pmatrix} y\_{t-p} \\ x\_{t-p} \end{pmatrix} + \begin{pmatrix} \varepsilon\_{11} \\ \varepsilon\_{21} \end{pmatrix} \tag{1}
$$

Considering the equation above, if all the coefficients <sup>ϕ</sup>(*q*) <sup>12</sup> (*q* = 1, 2, 3, ... , *p*) are zero, variable *x* is not Granger cause of *y*, which means that *x* cannot change *y*. GC test looks at whether lagged variable of one variable can be brought into alternative variable equations. In case that variable *x* can help explain variable *y*, then variable *x* is the Granger cause of variable *y*. A tool that helps judging the Granger cause is the F-test:

$$H\_0: \ \wp\_{12}^{(q)} = 0 \tag{2}$$

$$q = 1, 2, \dots, p \tag{3}$$

**Hypothesis 1.** *Existence of at least one q such that*ϕ(*q*) <sup>12</sup> -0:

$$S\_1 = \frac{(RSS\_0 - RSS\_1)/p}{RSS\_1/(T - 2p - 1)} \sim F(p, \ T - 2p - 1)\tag{4}$$

*where: RSS1 is the residual sum of squares in Equation (4) and RSS0 is the residual sum of squares in the y equation when* <sup>ϕ</sup>(*q*) <sup>12</sup> = 0.

$$RSS\_1 = \sum\_{t=1}^{T} \varepsilon\_{1t}^2 \tag{5}$$

The above statists fit in the F distribution.


The selection of *p* in VAR model is crucial while applying the GC test, since the results are strongly related to the number of lag orders. P is expected to be sufficiently large to reflect dynamic features but at the same time not excessively large since it brings several parameters for estimation and thus reducing the degrees of freedom (DOFs) of the model.

Akaike Information Criterion (*AIC*) and Schwarz Criterion (*SC*) are commonly used methods of selecting the number of lag orders. This selection should be conducted considering the quantity of lag items and DOFs. The calculations of those two methods are as follows:

$$AIC = \frac{-2l}{T} + \frac{2n}{T} \tag{6}$$

$$\text{SC} = \frac{-2l}{T} + n\frac{\ln T}{T} \tag{7}$$

where *n* is the total number of estimated parameters, *k* is the number of endogenous variables, *T* is the sample length, *d* is the number of lag orders, *l* is the logarithmic likelihood calculated as follows, by hypothesizing that the multivariate normal distribution is obeyed:

$$I = -\frac{Tk}{2}(1 + \ln 2\pi) - \frac{T}{2}\ln\left(\det\left(\frac{1}{T - m}\sum\_{t} \varepsilon\_{t} \varepsilon'\_{t}\right)\right) \tag{8}$$

Ideally, values of *AIC* and *SC* should be as low as possible.

#### *6.2. Granger Causality Network and Visualization*

Visualization of the connectivity among different time series is especially powerful when it comes to understanding the relation between different examined variables and how the one "Granger causes" the other.

In the network depiction of the corresponding G-causalities, nodes represent the different variable or system elements, in our case different time series and directed edges (arrows) represent causal interactions. Green lines depict unidirectional connections, with the direction of the arrow stating the relation between the nodes, meaning which one G-causes the other, while red lines depict bidirectional connections between two nodes. The thicker the line, the stronger the connection between the two corresponding nodes.

The same depiction, but in a matrix format, is also useful for the representation of causal interactions. In this setup, instead of lines there are colored boxes which correspond to two different nodes. The direction of the lines are replaced by the *x*- and *y*-axis with nodes of on *x*-axis being the nodes of the G-cause of the nodes on the *y*-axis. The darker the color of the box, the higher the G-causality effect.

#### 6.2.1. Causal Density

The causal density (*cd*) reflects the portion of interactions among nodes that are causally significant and provides a useful measure of system dynamical complexity. Causal density *cd* is defined as:

$$cd = \frac{\mathcal{gc}}{2N(N-1)}\tag{9}$$

where *gc* is the number of the total number of non-zero interactions observed, and *N* is the total number of nodes. The higher the value of *cd* the more the coordination of the system elements (nodes).

However, the term unit causal density (*cdu)* is often being used as a derivative of the *cd*. This is defined as the total number of significant interactions involving a specific node *i*. Nodes with relatively high values of *cdu* are considered to be hubs within the system.

#### 6.2.2. Causal Flow

Causal flow (CF) is a measure that helps defining if a node is a causal source or causal sink. It is defined as the difference between a node's in-degree and out-degree. It is another mean of network's representation which considers all the lines and their thickness in order to define if a node exerts causal influence or not. Nodes with positive CF are causal sources while those with negative CF are causal sinks.

#### *6.3. Data Preconditioning and Preprocessing in GCCA*

In order to perform a GCCA, a primary precondition must hold, that the variables must be CS (also known as wide-sense stationary). Covariance Stationarity is equivalent to saying that the first and second statistical moments (mean and variance) of each variable are constant, not varying with time. Deviations from CS can be tested by detecting "unit roots" in the data. We used the ADF test to assess the presence of a unit root in the variables (Section 4). A variable is that CS exhibits a tendency to return to a constant mean (or to a deterministic trend line). According to this mean-reversion tendency, large values will tend to be followed by smaller values, and on the opposite, small by larger ones. The ADF test detects the absence of this behavior. For non-CS variables, the following steps are: (a) linear (deterministic) trends are removed and (b) unit roots are removed by differencing [28].

#### *6.4. Model Validation in GCCA*

The inferences in GCCA are valid only in the case that a MVAR model captures well the "hidden" correlation in the data.

The amount of the variance "explained" by the model (in terms of the adjusted sum-square error, or adjusted *R*2), is the simplest indicator of the model's adequacy.

Typically, a value *R*<sup>2</sup> < 0.30 indicates that the model cannot explain an adequate amount of the variance in the data.

Another, similar to the abovementioned indicator is the model consistency, defined in Reference [28]. Within the frame of GCCA, this indicator is estimated as:

$$C = 1 - \frac{|R\_S - R\_r|}{|R\_r|} \times 100\tag{10}$$

where *Rr* is the correlation vector of the real data and *Rs* the correlation vector of the data generated by the MVAR model. Generally, if *C* < 80%, a cause for concern regarding model's ability arises.

If the model captures effectively the data, the residuals of the MVAR model must resemble a white noise process, i.e., they must be serially uncorrelated.

The Durbin–Watson statistic *d* [62], tests this behavior, and is given by:

$$d = \frac{\sum\_{t=2}^{T} (\varepsilon\_t - \varepsilon\_{t-1}^2)}{\sum\_{t=2}^{T} \varepsilon\_t^2} \tag{11}$$

If *d* < 1.0 there may be, also, cause for concern.

#### *6.5. Research Assumptions*

Correlation coefficients (Tables 8 and 9) reveal some connections among CBTF. We tried to apply a more powerful tool (GC) to detect if there was any hidden relation and to identify the direction of such a relation.

#### **7. Results**

We tried a number of MVAR models (or case studies, see Section 4, Table 6) to study the interaction of the CBT variables between Greece–Italy–Bulgaria, although there was no physical interconnection between Italy and Bulgaria. The first complete model, A, consisted of all three countries in order to detect all possible interactions, direct and indirect (transit flows). The other two models, B and C, consisted of Greece–Italy and Greece–Bulgaria, respectively, in order to emphasize the direct interactions among interconnected countries.

From Table 10, we observe that the MVAR models can explain the largest amount of the variance of the dependent variable, based on the shown validity measures, for cases A and B, i.e., the interconnections that involve the variables as shown in Table 6(as described Table 4). Variables D9, D15, D16, D21, and D24 for model A and D11 and D12 for model B were entered in the respective MVAR models as first differenced variables to guarantee their covariance stationarity (Table 5). We focused on models A and B having the largest consistencies.


**Table 10.** Models' consistency.

#### *7.1. Study A: Cross-border Trading between Greece, Italy, and Bulgaria*

In case study A (model A), nodes 1, 5, 7, 11, 12, and 16 were identified as causal sources (positive bars in the causal flow bar-chart), while nodes 3, 6, 8, 9, 17, and 18 as causal sinks (negative bars) as can be seen in Figure 4.

In general, we observed that node 5 (Italy's total load forecasted) was the most active variable (strongest driver), while node 8 (DA-Greece) was the most passive (biggest sink).

A connection was observed between nodes 1 and 3 (Italy's solar forecasted generation and Greece's solar forecasted generation, respectively) which is fairly explained due to the similar weather conditions of the two countries (diurnal effects).

Italy's load forecasting (node 5) seems to be strongly correlated with the first differences (changes) in Greece's and Bulgaria's load forecasting (nodes 6 and 18, respectively) which is explained mostly due to the similar weather conditions as well.

**Figure 4.** Study A: Demonstration of the GCCA. Matrix and network representations of the corresponding G-causalities are shown in the top panels. Causal flow and unit causal density are shown in the bottom panels.

Furthermore, we observed that DA Italy (node 7) affected the commercial schedule from Italy to Greece (node 9), since the Italian market was highly volatile and liquid, so it faced large excursions (spikes) in its wholesale electricity prices which, in turn, affected the power exchanged between the two countries.

Regarding DA Greece (node 8), which as mentioned was the most passive variable in the specific model, we observed that it was moderately correlated and affected by transfer capacity from Greece to Italy (node 11) and from Italy to Greece (node 12) which as shown in Table 8 were significantly highly correlated (0.998). A rational explanation seems that when the Greece–Italy interconnection was operational, NTC (500MW) comprised a significant value relative to the total demand in Greece. Consequently, DA Greece was highly affected whether Greece–Italy interconnection was operational or not.

A strong bi-directional relation was observed between solar generation in Bulgaria (node 13) and commercial schedules from Bulgaria to Greece which is also explained by the fact that Greece is a major importer from Bulgaria, since the latter has one of the lowest market clearing prices in EU driven by the low marginal costs of its nuclear plants and the low demand needs. So, it is rational to assume that a change in the solar generation of Bulgaria represents an excessive energy production which normally is exported by Bulgarian companies, taking into account that domestic demand for energy is covered by the conventional plants.

Finally, no relation was illustrated between the commercial schedule of Bulgaria and Greece (nodes 16 and 17) and DA prices in Bulgaria (node 15) which firstly seems rather strange. A possible explanation could be the wide margin between DA prices in Greece (higher prices) and Bulgaria. This fact leads electricity traders to commit to long-term contracts exporting energy from Bulgaria to Greece, as there is no commercial risk in that.

#### *7.2. Study B: Cross-Border Trading between Greece and Italy*

In case study B (model B), nodes 1, 2, 5, 7, and 11 were identified as casual sources (positive bars in the casual flow bar-chart), while nodes 3, 8, 9, and 10 as casual sinks (negative bars) as can be seen in Figure 5.

**Figure 5.** Study B: Demonstration of the GCCA. Matrix and network representations of the corresponding G-causalities are shown in the top panels. Causal flow and unit causal density are shown in the bottom panels.

In general, we observed that node 2 (wind generation in South Italy) was the strongest driver of the specific network with nodes 8 and 11 (DA Greece and commercial schedules from Greece to Italy, respectively) being the most remarkable sinks. This reflection may hide a relation of active and passive roles in the specific interconnection.

The solar forecasted generation, in South Italy (node 1), was strongly correlated to the Greek forecasted solar generation (node 3). A moderate correlation was also observed between South Italy's and Greece's wind generation (nodes 2 and 4). This seems reasonable, taking in consideration the geographical locations and the meteorological conditions of the two regions, South Italy and Greece, having a very similar climate, and it is something we observed from the correlation coefficient in Table 8. However, as mentioned earlier, in the case of wind and solar generation there was no point to discuss the Granger causality, since they are highly stochastic time series. The forecasted South Italy wind generation (node 2) also Granger causes strongly the commercial schedule from Greece to Italy (node 10).

We observed that node 7 (the change in DA price in South Italy), Granger causes nodes 9 (commercial schedule IT-GR) and 10 (commercial schedule GR–IT). But node 7 was influenced or Granger caused by node 5 (Italy's total load forecast). Also, node 5 had a bi-directional connection with node 6 (Greek total load forecast). This was something expected, since when the wholesale price in Italy was high, traders tried to buy cheap electricity for neighboring interconnected countries, hence the commercial schedules were directly affected. So, load forecasts in both countries Granger cause the DA price in South Italy (which means that the error in forecasting Italy's load were minimized when considering Greek load forecasts) which, in turn, Granger causes the total commercial schedules between IT–GR and vice versa. This is a very interesting result, indicating that the spot price in South Italy drives the commercial programs in the cross-border trading of both countries.

Node 4 (Greek forecast wind generation) Granger causes node 8 the changes in DA Greek price, as expected (wind generation reduces, in general, Spot prices). The changes in the DA Greek price were driven (Granger caused) by node 11 (the transfer capacity from Greece to Italy) which had a bi-directional connection to the node 12 (transfer capacity from Italy to Greece).

Regarding the casual density of variables under consideration, nodes 5 (Italian forecast load), 7 (changes in DA-price of South Italy) and 11 (transfer capacity from Greece to Italy) showcase the largest values, reflecting their total amount of casual interactivity, which is also a form of dynamical complexity in the network. These three variables were globally coordinated in their activity, within the system produced by model B. This means that they are useful in predicting each other's activity, as we have seen, but also useful in predicting other variables with different power.

The South Italian forecasted wind generation is the largest casual source, i.e., the variable (node) that exerts the strongest casual influence on the system as a whole, since Granger causes variables 10 and is strongly correlated to 4 (total commercial schedule GR–IT and Greek forecasted wind generation).

#### **8. Conclusions**

In the current study, we tried to identify and explain any causalities among energy fundamentals between Greece, Italy, and Bulgaria, using Granger causality theory. Causal connectivity implies a connection that lies on the enhancement of forecasting between two time series. That being said, and according to our knowledge, the lack of similar research on non-market coupled countries makes the current work unique and the outcomes that arise should be taken under consideration especially in the era of the implementation of the ETM.

We designed two models of which the one was reliable (highly consistent) enough to be analyzed. The model that was presented is the one that examines the interconnection between Greece and Italy. The authors could not design a reliable model for the Greece–Bulgaria interconnection. That might be due to the lack of data for the Bulgarian DA price (D21) which forced us to use half of the available data for the rest of the time series as well (in order to have the same number of observations) when examining the specific interconnection.

The results show that there were some causal connections among some of the electricity trading fundamentals which indicate that in order to facilitate the forecasting and minimize errors, we need to include them in the same model. There were also some weaknesses revealed which bring to the surface the need for market coupling among those countries.

It seems that the Italian market has a more active role in the specific interconnection followed by Greece. This could be explained by the fact that Italy is already part of the Central Western European (CWE) region which includes 19 market coupled countries, while Greece is not coupled yet with any of the interconnected countries and thus not so volatile and liquid like the Italian one

The lack of designing a reliable model that explains the causalities between Greece and Bulgaria offers opportunities for further research. Also, it would be interesting to see how the imminent market coupling between Greece and Italy will affect the relations described above.

**Author Contributions:** G.P.P. with C.D. and C.K. contributed significantly to the conceptualization, development, methodology and formal analysis of the paper. Moreover the former authors have contributed to software development, supervision and project administration as well. G.E. and F.G. contributed to Writing-Review & Editing, Data Curation and supervision. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding

**Acknowledgments:** Special thanks to ADMIE's Market Management Department for sharing their experience and knowledge with regards to the cross border trading.

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

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **An Integrated Energy Simulation Model for Buildings**

#### **Nikolaos Kampelis 1,\*, Georgios I. Papayiannis 2,3, Dionysia Kolokotsa 1, Georgios N. Galanis 2, Daniela Isidori 4, Cristina Cristalli <sup>4</sup> and Athanasios N. Yannacopoulos <sup>3</sup>**


Received: 14 January 2020; Accepted:25 February 2020; Published: 4 March 2020

**Abstract:** The operation of buildings is linked to approximately 36% of the global energy consumption, 40% of greenhouse gas emissions, and climate change. Assessing the energy consumption and efficiency of buildings is a complex task addressed by a variety of methods. Building energy modeling is among the dominant methodologies in evaluating the energy efficiency of buildings commonly applied for evaluating design and renovation energy efficiency measures. Although building energy modeling is a valuable tool, it is rarely the case that simulation results are assessed against the building's actual energy performance. In this context, the simulation results of the HVAC energy consumption in the case of a smart industrial near-zero energy building are used to explore areas of uncertainty and deviation of the building energy model against measured data. Initial model results are improved based on a trial and error approach to minimize deviation based on key identified parameters. In addition, a novel approach based on functional shape modeling and Kalman filtering is developed and applied to further minimize systematic discrepancies. Results indicate a significant initial performance gap between the initial model and the actual energy consumption. The efficiency and the effectiveness of the developed integrated model is highlighted.

**Keywords:** deformable models; electric energy demand; functional statistics; Kalman filtering; shape-invariant model

#### **1. Introduction**

Energy consumption in the building sector (combined with buildings construction) is associated with 36% of global final energy consumption and approximately 40% of total direct and indirect CO2 emissions. Associated figures rise every year mainly due to (a) improved energy access levels in developing countries, (b) increased ownership levels of energy consuming devices, and (c) the growth of global buildings floor area [1]. Measures to reduce energy consumption at building level include passive and active energy efficiency measures, storage, energy management, and building integrated renewable energy systems (e.g., solar, geothermal, and wind). Design and modeling of integrated energy systems in net zero energy buildings is discussed by Athienitis and O'Brien [2].

Assessing energy efficiency in buildings is a complex task which varies according to the aim of the analysis and the specificity of each case. In any case, expert knowledge and a set of technical and non-technical information is required. Technical information usually concerns the geometry and thermal characteristics of the building envelope, the design and operation of HVAC system, as well as data from measurements regarding indoor/outdoor climate conditions and energy end usage. Non-technical data in some cases are associated with occupancy levels, clothing levels, users behavior, perceptions, personal economics, and preferences.

Several decades now, significant research efforts have been directed towards evaluating energy efficiency in buildings. Research results have been applied in advancing state-of-the-art, knowledge and understanding, fostering new policies, better regulation, and innovations. In this context, several methods have been developed and evolved along with the creation of custom software applications built to match the needs of certain research and development fields. Various new techniques have been developed and applied to include the physical or "white box" approach, the statistical or machine learning "black box" and the hybrid "gray box" approach. A comparison of the models and applications for predicting building energy consumption in terms of complexity, easiness to use, running speed, input requirements, and accuracy is conducted by Zhao and Magoulès [3].

Physical models are developed to evaluate energy consumption in buildings by using equations of heat and mass transfer between the building and the surrounding environment as well as energy conservation among the building spaces and system components. Such physical models are classified based on the deployed approach and whether it is defined as single (well-mixed) zone [4], multi-zone, or zonal as provided in the state-of-the-art building modeling and energy prediction review by Foucquier et al. [5]. The zonal method is a simplification of CFD in the sense that the thermal zone is divided into several cells and in some cases representation in 2D is feasible. The advantage of the zonal approach is linked to the capability of dealing with large volumes in moderate computation time. SPARK [6] is such an example of a zonal approach software application. According to the multi-zone approach, each one zone or building element (e.g., wall and window) or system (e.g., HVAC system) or specific load (e.g., due to occupancy) is considered as one node for which the heat transfer equations are calculated. Each thermal zone is considered as homogeneous and represented by uniform state variables such as temperature, pressure, relative humidity, etc. The multi-zone or nodal approach is particularly effective when evaluating the energy performance of a building with many thermal zones since computational time for a year round simulation is relatively small. Therefore, the multi-zone approach is suitable for testing the impact of alternative energy efficiency measures provided that a reliable validated model has been created. The main disadvantage of this approach is related to the difficulty in creating a valid building model especially in the absence of holistic information of the building as built, systems installed, operational aspects, and data from measurements. EnergyPlus, ESP-r, and TrnSys are among the most robust and frequently used multi-zone software programs. Benchmarking between building energy simulation software programs reported in the literature is available by Harish et al. [7].

Data-driven methods on the other hand require no physical information, i.e., thermal or geometrical parameters, as they do not deploy heat transfer equations. Regression, Artificial Neural Network (ANN), Genetic Algorithm (GA), and Support Vector Machine (SVM) are some of the techniques used in building energy forecasting based solely on measurements of parameters such as temperature, relative humidity, solar radiation, wind velocity, and energy consumption/production. Machine learning techniques for estimating building energy consumption are exploited by Naji et al. [8] and by Robinson et al. [9]. The main drawback of data-driven methods concerns the interpretation of results in physical terms. Data-driven methods for building energy prediction and classification studies are reviewed by Amasyali and El-Gohary [10] and by Wei et al [11]. Hybrid or "gray box" methods combine physical modeling with data-driven techniques to counterbalance the weaknesses of the "white" and "black" box approaches. Machine learning techniques can be used for parameter estimation, e.g., by coupling a multi-zone model with GA. Another hybrid approach is to use statistical models to improve the performance of physical models with respect to end-uses, which are often unknown and hard to be modeled in a deterministic way.

In contrast with data-driven methods, detailed simulation models do not require conditions monitoring to predict the building performance. However, various sources of uncertainty can lead to significant discrepancies between model predictions and metered energy use. Such sources of uncertainty can be distinguished to specification-related, modeling-related, and scenario-related as discussed by De Wit and Augenbroe [12]. Furthermore, embedding realistic occupant behavior on building modeling and its impact on energy predictions is investigated and addressed as a significant and complex problem by Ryan and Sanquist [13]. Overall, performance gaps can be attributed to a variety of discrepancies in building modeling. A methodology to identify building energy performance problems related to operational, measurements, or simulation aspects is proposed by Maile et al. [14].

Essentially, validation is a critical step to assess the simulation model's plausibility and reliability. This is especially important for the investigation of energy efficiency measures and the assessment of the actual (not the relative) potential impact in terms of energy and cost savings, environmental performance, thermal comfort, etc. Methods for model validation (otherwise referred to as calibration) are reviewed by Coakley et al. [15], and indicators used in relevant standards (ASHRAE Guideline 14, IPMVP, FEMP) to address acceptance thresholds are discussed [16] by Royapoor and Roskilly [17]. Calibration of a building energy model based on actual measurements is extensively investigated by O'Neil et al. [18]. The building under study hosts offices and conferences rooms, and the authors perform an extensive sensitivity analysis including more than 2000 parameters automatically refined using analytic meta-model-based optimization. A similar approach is followed by O'Neil et al. [19] incorporating EnergyPlus and TRNSYS in their investigation and over 1000 parameters for model calibration.

This paper investigates (i) the validation of a Near-Zero Energy Building (NZEB) simulation model through trial and error approach of important model parameters and (ii) a novel postprocessing approach is proposed which improves the simulation model's accuracy, combining techniques of functional statistics through appropriate energy shape modeling (employing the concept of deformable models) and Kalman filtering to reduce the remaining biases. Initially, a physical model of the Leaf Lab industrial NZEB in Italy is used to conduct a simulation of the building's energy consumption for two consecutive years using weather data recorded from onsite meteorological stations. HVAC electric consumption obtained from the simulation is compared to actual measured values to establish the initial (baseline) performance gap. Subsequently, through a trial and error approach, important model parameters are identified and fine-tuned until an improved acceptable correlation between simulated and recorded HVAC electric energy consumption is reached. At a second stage, inconsistencies in the shape of the optimized energy prediction are corrected through an appropriate functional shape/reshape modeling task by comparing recent measured energy demand outputs to predictions and appropriately estimating expected deviations which are used to improve the shape model's output. At the last step, Kalman filtering is incorporated to remove remaining systematic biases where is needed. Therefore, the original physical model's results are passed through this integrated model where meaningful interventions are performed to better predict the true situation without neglecting the physical interpretation of the model output. The benefits from the proposed approach are illustrated in the later section (Section 3) with substantial reduction in the error magnitude.

#### **2. Methodology and Modeling Approaches**

In this section, we describe the modeling approaches that are used throughout the paper for the prediction of the energy demand. In particular we present (a) the initial energy simulation model and its optimized version by appropriate parameter tuning, (b) the shape model approach in which the optimized output of the energy simulation model is reshaped to further reduce systematic inefficiencies related to discrepancies from the actual shape of the energy demand, and (c) postprocessing using Kalman filtering which reduces remaining systematic errors not captured by the shape model approach. A flowchart of the methodology is presented in Figure 1.

**Figure 1.** Flowchart of the high level methodology followed.

#### *2.1. The Energy Simulation Model: General Presentation of the Model, Areas of Application, Advantages, and Shortcomings*

EnergyPlus is the simulation engine software used to conduct an integrated simulation of the building, system, and plant whereby supply and demand are matched based on successive iteration substitution following Gauss–Seidel updating [20]. Open Studio is used as the API software for developing and parameterizing the model following the principles outlined by Brackney et al. [21]. Ambient temperature, relative humidity, solar radiation, and wind speed data for 2017 and 2018 was obtained from local meteorological equipment and converted to two yearly weather files. Data of total HVAC energy consumption, for the years 2017 and 2018, are exploited for providing the baseline against which model based results are evaluated.

The simulation model contains, on the one hand, the geometry, construction components and materials of the building under study. For opaque material thickness (m), thermal conductivity (W/mK), density (kg/m3), and thermal absorptance (*dimensionless*) properties are edited. For transparent materials, such as glass in windows and sky windows thickness (m), thermal conductivity (W/mK) and optical properties, such as solar, visible, and infrared transmittance, are inserted. On the other hand, a model of the HVAC system is designed based on the installed technologies and adjusted accordingly to the actual key performance heat pump technical parameters such as Coefficients of Performance (COP), fan maximum flow power (m3/s), pressure rise, and efficiency. Other parameters such as rated total heating/cooling capacity, and rated and maximum air flow rated are automatically sized based on the software's calculations. Furthermore, with respect to the operation of the major installed systems, the simulation model takes into account the temperature set points of the HVAC system, ventilation, and infiltration rates (ACH−1) and a number of schedules to determine artificial lighting, electric equipment, and occupancy. Moreover, humidity control is exercised using appropriate schedules controlling the operation of the HVAC to ensure that the relative humidity during working hours in the various thermal zones varies between 40% and 60%. Subsequently, an intensive search of the parameters that affected the daily, monthly, and annual power distribution profiles is followed to improve the initial results of the model based by minimizing deviation from HVAC power consumption data. Through the trial and error various combinations and fine-tuning of the all of the above parameters is carried out to reach the optimum results when assessing intra-day, monthly, and annual deviation levels.

EnergyPlus simulation is based on heat balance calculations solved simultaneously with the aid of on an integration solution manager, which includes surface heat balance, air heat balance, and building systems simulation blocks. The heat balance of outside surfaces is calculated based on the equation

$$q\_{asol}^{\prime\prime} + q\_{L\mathcal{W}R}^{\prime\prime} + q\_{conv}^{\prime\prime} - q\_{ko}^{\prime\prime} = 0 \tag{1}$$

where

*q <sup>α</sup>sol* is the absorbed direct and diffuse solar (short wavelength) radiation and heat flux *q LWR* is the net long wavelength (thermal) radiation flux exchange with the air and surroundings *q conv* is the convective flux exchange with the outside air *q ko* is the conduction heat flux (q/A) into the wall

Clearly, *q <sup>α</sup>sol* is influenced by parameters such as location, surface angle and tilt, surface material, and weather conditions. *q LWR* is determined by radiation exchange between the surface and the ground, sky and air. It is dependent on the absorptivity and emissivity of the surface; the temperature of the surface, sky, ground, and air; and corresponding view factors. Assumptions such that each surface is at uniform temperature and energy flux leaving a surface is evenly distributed are considered reasonable for building energy simulation. Using the Stefan–Boltzmann Law in the above equation yields

$$q\_{\rm LNR}^{\prime\prime} = \varepsilon \sigma F\_{\rm gnd} (T\_{\rm gnd}^4 - T\_{\rm surf}^4) + \varepsilon \sigma F\_{\rm sky} (T\_{\rm sky}^4 - T\_{\rm surf}^4) + \varepsilon \sigma F\_{\rm air} (T\_{\rm air}^4 - T\_{\rm surf}^4) \tag{2}$$

where

 is the long-wave emittance of the surface *σ* is the Stefan–Boltzmann constant *Fgnd* is the view factor of wall surface to ground surface temperature *Fsky* is the view factor of wall surface to sky temperature *Fair* is the view factor of wall surface to air temperature *Tsur f* is the outside surface temperature *Tgnd* is the ground surface temperature *Tsky* is the sky temperature *Tair* is the air temperature

The above equation is converted by introducing linear radiative heat transfer coefficients such that

$$q\_{\rm LNR}^{\prime\prime} = h\_{r,\rm gnd}(T\_{\rm gnd} - T\_{\rm surf}) + h\_{r,\rm sky}(T\_{\rm sky} - T\_{\rm surf}) + h\_{r,\rm air}(T\_{\rm air} - T\_{\rm surf}) \tag{3}$$

where

$$h\_{r, \text{gnd}} = \epsilon \sigma F\_{\text{gnd}} (T\_{\text{surf}}^4 - T\_{\text{gnd}}^4) / (T\_{\text{surf}} - T\_{\text{gnd}}) \tag{4}$$

$$h\_{r, \text{gnd}} = \varepsilon \sigma F\_{\text{gnd}} (T\_{surf}^4 - T\_{\text{sky}}^4) / (T\_{surf} - T\_{\text{sky}}) \tag{5}$$

$$h\_{r, \text{gud}} = \epsilon \sigma F\_{\text{gud}} (T\_{surf}^4 - T\_{air}^4) / (T\_{surf} - T\_{air}) \tag{6}$$

Exterior convection is modeled using equation

$$q\_{conv}^{\prime\prime} = h\_{\varepsilon, \text{ext}} A (T\_{surf} - T\_{air}) \tag{7}$$

where

*q conv* is the rate of exterior convective heat transfer *hc*,*ext* is the exterior convection coefficient *A* is the surface area *Tsur f* is the surface temperature *Tair* is the outdoor air temperature

Conduction heat fluxes are modeled based on equation

$$q\_{ko}^{''}(t) = \sum\_{j=0}^{\infty} X\_j T\_{o, t - j\delta} - \sum\_{j=0}^{\infty} Y\_j T\_{i, t - j\delta} \tag{8}$$

where

*q ko*(*t*) is the conductive heat flux for the current time step *T* is temperature *i* indicates the internal element of the building *o* indicates the external element of the building *X*,*Y* are the response factors

In more detail, Conduction Transfer Functions (CTFs) as shown in (9) and (10) below are used to estimate the heat fluxes on either side of the building elements based on previous temperature values of interior and exterior surfaces as well as previous interior flux values.

$$q\_{\rm ki}^{''}(t) = -Z\_{o}T\_{i,t} - \sum\_{j=1}^{nz} Z\_{j}T\_{i,t-j\delta} + \chi\_{o}T\_{o,t} + \sum\_{j=1}^{nz} \chi\_{j}T\_{o,t-j\delta} + \sum\_{j=1}^{nq} \Phi\_{j}q\_{ko,t-j\delta}^{''} \tag{9}$$

$$q\_{k\nu}''(t) = -Y\_o T\_{i,t} - \sum\_{j=1}^{nz} Y\_j T\_{i,t-j\delta} + X\_o T\_{o,t} + \sum\_{j=1}^{nz} X\_j T\_{o,t-j\delta} + \sum\_{j=1}^{nq} \Phi\_j q\_{k\nu, t-j\delta}^{''} \tag{10}$$

where

*Xj* is the outside CTF coefficient, j = 0,1,...nz *Yj* is the cross CTF coefficient, j = 0,1,...nz *Zj* is the inside CTF coefficient, j = 0,1,...nz *φ<sup>j</sup>* is the flux CTF coefficient, j = 0,1,...nq *Ti* is the inside surface temperature *To* is the outside surface temperature *q ko* is the conduction heat flux on the outside face *q ki* is the conduction heat flux on the inside face

In addition, for each thermal zone EnergyPlus simulation is based on an integration of energy and moisture balance as shown in the Equation (11) below.

$$\mathbf{C}\_{z}\frac{dT\_{z}}{dt} = \sum\_{i=1}^{N\_{ol}}\mathbf{Q}\_{i} + \sum\_{i=1}^{N\_{surfacc}}h\_{i}A\_{i}(T\_{si} - T\_{z}) + \sum\_{i=1}^{N\_{surfacc}}\dot{m}\_{i}\mathbf{C}\_{p}(T\_{zi} - T\_{z}) + \dot{m}\_{\text{inf}}\mathbf{C}\_{p}(T\_{zi} - T\_{z}) + \dot{Q}\_{\text{sys}} \tag{11}$$

where

∑*Nsl <sup>i</sup>*=<sup>1</sup> *<sup>Q</sup>*˙ *<sup>i</sup>* is the sum of convective heat transfer from the zone surfaces <sup>∑</sup>*Nsur f aces <sup>i</sup>*=<sup>1</sup> *hiAi*(*Tsi* − *Tz*) is the convective heat transfer from the zone surfaces *m*˙ *inf Cp*(*Tzi* − *Tz*) is the heat transfer due to infiltration of outside air <sup>∑</sup>*Nsur f aces <sup>i</sup>*=<sup>1</sup> *m*˙ *iCp*(*Tzi* − *Tz*) is the heat transfer due to interzone air mixing *Q*˙ *sys* is the air systems output *Cz dTz dt* is the energy stored in zone air, and *Cz* = *ρairCpCT*

Infiltration is outdoor air unintentionally entering the building due to the opening of doors as well as air leakage through windows and other openings. Infiltrated air is mixed with air in the various thermal zones of the building. Determining infiltration (or air tightness) values contains significant uncertainty, as it requires a complex and elaborate procedure often referred to as blower door test. Infiltrated air is commonly modeled as the number of air changes per hour (ACH) and taken into account in the air heat balance at temperature equal to that of ambient air. In EnergyPlus, infiltration is modeled based on Equation (12).

$$\text{Infibration} = \left(I\_{\text{design}}\right)\left(F\_{\text{schedule}}\right)\left[A + B\left|\left(T\_{\text{zone}} - T\_{\text{odd}}\right)\right| + C\left(\text{Windspeed}\right) + D\left(\text{Windspeed}\right)^2\right] \tag{12}$$

where

*Idesign* is the user defined infiltration value (ACH−1) *Tzone* is the zone air temperature at current conditions (deg C) *Todb* is the outdoor air dry-bulb temperature (deg C) *Fschedule* is a user defined schedule value between 0 and 1 *A* is the constant term coefficient *B* is the temperature term coefficient *C* is the velocity term coefficient *D* is the velocity squared coefficient

Similarly, ventilation can be modeled using a schedule, maximum and minimum values, as well as delta temperature values, and is determined by the equation

$$\text{Ventilation} = (V\_{\text{design}}) (F\_{\text{schedule}}) [A + B | (T\_{\text{zone}} - T\_{\text{adv}}) | + \mathcal{C} (\text{Windspred}) + D (\text{Vindspred})^2] \tag{13}$$

where

*Vdesign* is the user defined ventilation value (ACH−1) *Tzone* is the zone air temperature at current conditions (deg C) *Todb* is the outdoor air dry-bulb temperature (deg C) *Fschedule* is a user defined schedule value between 0 and 1 *A* is the constant term coefficient *B* is the temperature term coefficient *C* is the velocity term coefficient *D* is the velocity squared coefficient

Furthermore, the energy provided to each thermal zone by the HVAC system, *Q*˙ *sys*, is given by Equation (14).

$$
\dot{Q}\_{\text{sys}} = \dot{m}\_{\text{sys}} \mathcal{C}\_p (T\_{\text{sys}} - T\_z) \tag{14}
$$

Equations (11) and (14) can be transformed to yield zone air temperature as shown in Equation (15) below.

$$T\_z^t = \frac{\sum\_{i=1}^{N\_{\text{sl}}} Q\_i^i + \dot{m}\_{\text{sys}} \mathbb{C}\_p T\_{\text{supply}}^t + \left(\mathbb{C}\_z \frac{T\_z}{t} + \sum\_{i=1}^{N\_{\text{surfactor}}} h\_i A\_i T\_{\text{si}} + \sum\_{i=1}^{N\_{\text{resins}}} \dot{m}\_i \mathbb{C}\_p T\_{\text{zi}} + \dot{m}\_{\text{inf}} \mathbb{C}\_p T\_{\text{so}}\right)^{t-\delta t}}{\frac{\mathbb{C}\_z}{\delta t} + \left(\sum\_{i=1}^{N\_{\text{surfactor}}} h\_i A\_i + \sum\_{i=1}^{N\_{\text{resins}}} \dot{m}\_i \mathbb{C}\_p + \dot{m}\_{\text{inf}} \mathbb{C}\_p + \dot{m}\_{\text{sys}} \mathbb{C}\_p\right)} \tag{15}$$

*2.2. Shape Modeling and Systematic Inefficiencies Correction of the Prediction Model: Presentation of the Properties and Capabilities of the Shape Invariant Model and Implementation in the Current Study*

At the first stage, the energy simulation model was carefully tuned to reduce prediction errors. However, there are some sources of error that cannot be effectively treated only through parameter tuning. Therefore, advanced statistical modeling approaches are considered to further reduce the prediction model's deviance from the true situation. At this postprocess stage, the actual phenomena of energy and mass transfer between the environment and the building as well as within the building itself are not explicitly modeled. Instead, this is a functional modeling approach, wherein the shape of the prediction for the energy demand is appropriately modeled and rearranged to better approximate the actual (measured) energy shape. As a result, the proposed shape modeling approach obsoletes deficiencies caused on the difference of the shape between prediction–reality while the remaining biases are further treated through appropriate Kalman filtering procedures discussed in Section 2.3.

#### 2.2.1. The Shape Model Approach

Let us denote by (*t*, *Xj*(*t*)) the observations representing the energy demand at a specified day *j* where *t* ∈ [0, 24) represents the hour of the day and *Xj*(*t*) the observed energy demand at time instant *t*. In general, these measurements are available only on certain time segments (e.g., per hour); therefore, the accurate daily shape is not known. That means that although the true shape of the energy demand is a continuous function with respect to time parameter *t*, in practice only some points of this shape are known at specified time segments *ti*; therefore, the available data are of the form {*ti*, *Xj*(*ti*)}*<sup>n</sup> <sup>i</sup>*=1, which can be considered as landmarks. As we are interested in working with the shapes of the daily energy demands, we need to consider the shape of the day *j* as a function with respect to time, i.e., *fj* : [0, 24) <sup>→</sup> <sup>R</sup>. Then, using the available data {*ti*, *Xj*(*ti*)}*<sup>n</sup> <sup>i</sup>*=<sup>1</sup> from day *j* we are able to estimate a smoothed version of the energy shape employing any typical interpolation method or nonparametric filters (e.g., spline smoothers, kernel-based smoothing methods, etc. [22–24]) by choosing appropriately the mollification parameters in order not to lose important aspects of the information. In this manner, the shape function of the intra-day power demand is sufficiently recovered with the advantage that we can get estimates for the demand even in time instants that no data are available and allowing to treat data with functional statistics techniques.

The daily shape of the energy demand is not expected to change dramatically between two days given that we consider typical working days (i.e., not weekends or public holidays). As a result, a standard energy-demand picture is expected to be observed with small fluctuations from one day to another, where these fluctuations can be efficiently calibrated by appropriate shape model considerations. Consider, for example, that for an arbitrary day *j*, *fj*(*t*) denotes the observed energy demand and *f <sup>j</sup>*(*t*) denotes the prediction obtained by the simulation model discussed in Section 2.1. Clearly, as the simulation model prediction is not expected to coincide with the true state of the energy demand, if systematic inefficiencies are presented between the prediction and reality, then a shape correction procedure could remarkably reduce errors caused to daily energy shape deviances. We discuss an approach under which we expect to provide corrections to the simulation model prediction by properly "reshaping" the simulation output in order to better match the observed energy shapes based on previous data. Such an approach is possible under the framework of deformation models (see, e.g., [25–28]) where the observed function *fj* (i.e., observed energy shape) is considered as a deformation of the prediction model *f <sup>j</sup>* (predicted energy shape), and this relation is mathematically expressed through the model

$$f\_j(t) = \mathcal{R}\_j(f\_j(t)) + \epsilon\_j(t) \tag{16}$$

where *Rj* : <sup>R</sup> <sup>→</sup> <sup>R</sup> is called a deformation function and *<sup>j</sup>*(*t*) is considered as a white noise process. Although several models can be proposed to parameterize the deformation function (e.g., shape-invariant model [25–27]), for the particular nature of the data we consider in this paper, we may propose a simpler model which consists of modeling the reshape function *α<sup>j</sup>* defined by

$$a\_{\hat{f}}(t) := [\widetilde{f}\_{\hat{f}}(t)]^{-1} f\_{\hat{f}}(t). \tag{17}$$

If the modeler had knowledge of the reshape function, then he/she could perfectly adjust the initial prediction *f <sup>j</sup>*(*t*), provided by the simulation model discussed in Section 2.1, to obtain exactly the true energy demand *fj*(*t*). In particular, knowledge of the exact functional form for *α*(·) would allow predictions even at intermediate time instants for which observations are not available. In this section, we propose a functional statistical model to estimate the reshape function from past data of initial simulation model discrepancies from the actual measured energy shapes, so that it can be used for future predictions.

As the exact knowledge of the reshape function is not an option, we can estimate it using data from the previous days (we should choose a small time window ~5–10 days) to model the "typical" reshape function that is observed in the near past. Clearly, using the information provided from the last *N* days, let us define the set of reshape functions A := {*α*1, *α*2, ..., *αN*}. For the case under consideration, it is expected that there are certain aspects of the daily energy shape which the prediction model cannot efficiently calibrate and systematically does not capture, and as a result, the reshape functions must be

very similar objects as shown for example for a typical set of observations in Figure 2. In such a case, the approach we propose here is valid.

**Figure 2.** Example of the observed and simulated daily energy demand (per hour) and the corresponding reshape functions for the time period 22/01/2018 to 26/01/2018.

Under this consideration, for an appropriately chosen period of time each function, *α<sup>j</sup>* should not present significant fluctuations from the reshape function that is considered as the "typical" one. Therefore, an appropriate notion of mean regarding the set of reshape functions A is needed to properly define/represent the mean element in the set. The latter task requires the calculation of the mean element among a number of functions living in a space which does not necessarily has linear structure, therefore the notion of Fréchet mean needs to be exploited [29,30] for this purpose. In the context we discuss here, each reshape function *αj*(·) is considered as a deformation to the mean element (i.e., Fréchet mean) *α*¯(·) with respect to which an appropriate notion of deviance is defined and its minimization will provide the mean element (please see [31] for technical details on this subject). A general model like the *shape-invariant model* [26] can be used for the parameterization of the functional characteristics of the reshape functions in order to define a consistent parametric form for the Fréchet mean. According to the shape-invariant model, shape deformations like vertical (scale) and horizontal (time) shifts can be efficiently captured. The standard shape-invariant model applied to the energy reshape functions can be written as

$$
\pi\_j(t) = \beta\_j + \kappa\_j \mathbb{1}(t - \mathbb{Z}\_j) + \varepsilon\_j(t), \ \varepsilon\_j(t) \sim \text{WN}(0, \sigma^2) \tag{18}
$$

where *α*¯(*t*) denotes the mean pattern of the energy reshape function, *β<sup>j</sup>* and *κ<sup>j</sup>* introduce vertical shifts parameterization, while *ζ<sup>j</sup>* introduces time-shift parameterization. Recall that any *α<sup>j</sup>* is considered as a deformation of the observed mean pattern (Fréchet mean) of the set A. As a result, the incurred mean energy reshape function estimated by the data from the previous *N* days can be calculated through the equation

$$\mathbb{E}\left(t; \boldsymbol{\theta}^\*\right) = \frac{1}{N} \sum\_{j=1}^{N} \left(\frac{1}{\kappa\_j} \alpha\_j(t + \zeta\_j) - \beta\_j\right) \tag{19}$$

where vector *θ*∗ = (*κ*, *ζ*, *β*) = (*κ*1, ..., *κN*, *ζ*1, ..., *ζN*, *β*1, ..., *βN*) contains all the deformation parameters. Clearly, as these parameters uniquely define the mean reshape function must be selected to minimize the mean model variance from the set A, i.e., the vector *θ*<sup>∗</sup> is chosen as

$$\theta^\* := \arg\min\_{\theta \in \Theta} \Lambda(\theta) = \arg\min\_{\theta \in \Theta} \sum\_{j=1}^N \int\_{\mathcal{T}} \left( a\_j(t) - \mathbb{I}(t; \theta) \right)^2 dt \tag{20}$$

where T represents the time period within a day and Θ represents the space of the deformation parameters which are subjected to some normalization constraints, i.e.,

$$\Theta := \left\{ \theta \, : \, \sum\_{j=1}^{N} \beta\_j = 0, \, \sum\_{j=1}^{N} \zeta\_j = 0, \, \prod\_{j=1}^{N} \kappa\_j = 1, \, \kappa\_j \ge 0, \, \text{for all } j = 1, 2, \dots, N \right\}. \tag{21}$$

Clearly, the estimation of the reshape function for the day *j* = *N* + 1 will be used to improve the initial prediction, provided by the simulation model, for the energy demand of this day through the model

$$f\_{\dot{j}}(t) = \dot{f}\_{\dot{j}}(t) + \eta\_{\dot{j}}(t) := \mathbb{R}(t; \theta)\dot{f}\_{\dot{j}}(t) + \eta\_{\dot{j}}(t) \tag{22}$$

where the error term *ηj*(*t*) is considered as a white noise process. After the initial mean reshape function estimation from the first batch of data corresponding to the initial *N* days has been obtained, an exponentially weighted scheme can be used to update appropriately the reshape functions taking into account both the effect of the initial reshape function and more recent observations on the reshape function. In particular, the proposed scheme can be described through the following steps.

**Initial Step** Set *k* = 0 and provide a choice for *λ* ∈ (0, 1). Given the simulation model predictions { *f j*} and the corresponding measurements { *fj*} for *j* = 1, 2, ..., *N* estimate the reshape functions {*αj*} from (17). Then, set as *<sup>α</sup>*ˆ(*k*)(*t*) the Fréchet mean of {*αj*} (calculated by (19)–(20)) and provide the prediction ˆ *fj*(*t*) = *α*ˆ(*k*)(*t*)*f <sup>j</sup>*(*t*) for *j* = *N* + 1. For every new prediction task, repeat steps 1–3.

**Step 1** Given the new measurement *fj*(*t*) set *k* = *k* + 1, *α*0(*t*) := *α*ˆ(*k*−<sup>1</sup>)(*t*) and *α*1(*t*) := *f* −1 *<sup>j</sup>* (*t*)*fj*(*t*). **Step 2** Set as *<sup>α</sup>*ˆ(*k*)(*t*) the Fréchet mean of *<sup>α</sup>*0(*t*) and *<sup>α</sup>*1(*t*) with weights 1 <sup>−</sup> *<sup>λ</sup>* and *<sup>λ</sup>*, respectively. **Step 3** Given the simulation model prediction *f <sup>j</sup>*+1(*t*) provide the improved (reshaped) prediction

 ˆ *fj*+1(*t*) = *α*ˆ(*k*)(*t*)*f <sup>j</sup>*+1(*t*).

The exponential weighting is used to reduce the effect of older observations to the new predictions, i.e., for the prediction on the day *N* + 2 using the weight parameter *λ* ∈ (0, 1) we could weight the most recent information (last observed reshape function *αN*+1(·)) by *λ* and the older observations by (1 − *λ*). In this manner, choosing *λ* close to 1 we allocate more weight to the most recent realizations and reduce the effect of the older observations. Otherwise, choosing *λ* close to zero, we forget the older observations with a very slow rate allowing the past information to contribute more to the prediction. Clearly, the choice of this parameter is critical to the quality of the prediction, and the final choice of this parameter depends strongly to the nature of the application that the prediction model is employed to. Note also that at each step the Fréchet mean of previous reshape functions is taken into account as an observation through the term *α*0(*t*) although it is not. However, this modification allows to simultaneously condense and appropriately weight (according to our preferences provided by the choice of *λ*) the past information into a single term.

#### 2.2.2. The Weighted Shape Model Approach

In practice, it has been proved that there are periods of time that the prediction models do not provide reliable predictions and they may significantly deviate from the true situation. An example of such a situation is the energy demand prediction of the building during the period of summer holidays discussed in Section 3. In such cases, it is very important to quickly perceive when this happens and rapidly adjust the prediction to an acceptable level of deviance from the reality. For such purposes, we present here a small variation of the scheme presented in Section 2.2.1 where a weighted version of the reshape model is used.

The main idea is to divide the prediction into two parts: (a) the prediction provided by the reshape model and weight it by a proportion *w* ∈ (0, 1) and (b) the prediction provided by previous observed energy shapes (measurements only) weighted by the proportion 1 − *w*. For the second part, the same procedure that used for the estimation of the mean reshape function is used, i.e., given the past measurements of the daily energy shapes { *<sup>f</sup>*1, *<sup>f</sup>*2, ..., *fN*} their Fréchet mean ¯ *f*(*t*) is estimated by (19)–(20) substituting *α<sup>j</sup>* with*fj*. The Fréchet mean in this case is interpreted as the most typical energy shape observed in the previous days without taking into account any information provided by the energy simulation model (predictive model). Then, one could shape the prediction either by constantly setting the weighting parameter *w* to a specific value or by changing this value each time new data become available to adjust the weighted shape model as close as it is possible to the true situation as observed until that time instant. Such a weight allocation criterion could be constructed as follows. Define by *gj*+1(*t*; *w*) := *w* ˆ *fj*(*t*)+(<sup>1</sup> <sup>−</sup> *<sup>w</sup>*) ¯ *fj*(*t*) the weighted prediction, where ¯ *fj* denotes the mean energy shape as estimated until day *j* and ˆ *fj* denotes the shape model's prediction for the day *j* + 1 derived from the approach discussed in Section 2.2.1. Then, given the measurement of the day *j* + 1, *fj*+1(*t*), the weight *w* for the next prediction (i.e., the day *j* + 2) is chosen as the minimizer to the criterion

$$\min\_{w \in (0,1)} \int \left( g\_{j+1}(t;w) - f\_{j+1}(t) \right)^2 dt. \tag{23}$$

Clearly, if the prediction provided by the simulation model is completely misplaced, then this criterion will act rapidly as a safety filter and will provide a prediction that is based more on the empirical data (previously observed energy shapes); otherwise, the prediction will be based on the simulation model. Therefore, criterion (23) will act as a detection mechanism of significant inconsistencies of the estimates provided by the simulation model and if such a significant shift occurs, will immediately properly adjust the prediction and improve its accuracy. Moreover, monitoring the value of *w* for a reasonable period of time, provides a measure of efficiency for the prediction model that is used. Ideally, we expect the value of *w* to remain in high levels (near to 1) except of some periods that major changes happen where the simulation model has not yet re-adjusted.

#### *2.3. Postprocessing Using Kalman Filtering: The General Algorithm, Capabilities and Areas of Application, and the Filter Proposed for the Present Work*

Numerical simulation models in several applications, including energy prediction systems, are exposed to systematic or non-systematic biases, an issue in which a wide number of internal or external parameters are involved: inability of simulating sub-scale phenomena, limitations in the parameterization of all the engaged parameters, and numerical schemes problems can be listed among them.

Towards the limitation of the problems above and in the framework of integrated forecasting systems, statistical postprocesses and bias removal techniques have a critical role. In particular, Kalman filters [32–34] provide a very popular approach for systematic bias removal, combining recursively available records with direct modeled outputs, by using weights that minimize the corresponding biases, requiring limited CPU resources and data backlog. Kalman filters have been successfully applied to a wide number of numerical models including atmospheric parameters (see, for example, [35–39]), sea state (see, e.g., [40–42]), extreme events (see, e.g., [43,44]), as well as renewable power resources (see, e.g., [45–47]). A general description of the basic Kalman filter algorithm is summarized below.

The main target is the simulation of two parameters that evolve in time in parallel. The state vector *x* and the observation corresponding *y*. The change of the two processes in time are described by the *system* and the *observation* equations, respectively:

$$\mathbf{x}\_{t} = F\_{t}\mathbf{x}\_{t-1} + w\_{t} \tag{24}$$

$$y\_t = H\_t \mathbf{x}\_t + v\_t.\tag{25}$$

The coefficient matrices, *Ft* and *Ht*, are defined as the *system* and the *observation* matrix, respectively. The corresponding covariance matrices *Wt*, *Vt* of the random vectors *wt* and *vt*, respectively, should be determined before the application of the filter while *wt* and *vt* need to be independently distributed according to Gaussian probability laws. The Kalman filter theory provides a method for the recursive estimation of the unknown state *xt* utilizing all the observation values *y* up to time *t*. The following equations describe a full integration step of the Kalman algorithm.

$$\mathbf{x}\_{t} = F\_{t}\mathbf{x}\_{t-1} + K\_{t}(y\_{t} - H\_{t}\mathbf{x}\_{t|t-1}),\tag{26}$$

where

$$K\_t = P\_{t|t-1} H\_t^T (H\_t P\_{t|t-1} H\_t^T + V\_t)^{-1} \tag{27}$$

and the covariance matrix of the unknown state *x* is given by

$$P\_{t|t-1} = F\_t P\_{t-1} F\_t^T + \mathcal{W}\_t \tag{28}$$

where

$$P\_t = (I - K\_t H\_t) P\_{t|t-1} \,. \tag{29}$$

In the present work, a linear filter has been adopted for the reduction of possible systematic biases in previous simulations steps. In particular, the estimation of the bias *yt* in time is estimated by means of the direct model output *mt* (where *mt* ):

#### *yt* = *x*0,*<sup>t</sup>* + *x*1,*tmt* + *vt*

The above provides the observation equation of the filter with observation matrix *Ht* = [1 *mt*] and state vector *xt* = [*x*0,*<sup>t</sup> x*1,*t*]. The covariance matrices *Vt* and *Wt* of the state and observation errors are estimated by utilizing a training window for the observed and modeled data (see [38,41,42]). The estimated by the filter bias values are utilized for the improvement of model outputs in the next time steps. It is worth noticing that the optimum training windows as well as initial values are case sensitive and should be estimated separately for every application under study.

#### **3. Test Cases and Results**

In this section, the proposed integrated model presented in Section 2 is implemented to the prediction of the daily and hourly energy demand of Leaf Lab for the time period 2017–2018.

The Leaf Lab is an industrial Near-Zero Energy Building (NZEB) of 6000 m<sup>2</sup> total floor area, integrating energy efficiency measures, advanced automations, renewable energy generation, and storage. The building envelope is highly insulated and consists of walls and double glazed windows with U values of 0.226 W/m2K and 1.793–3.194W/m2K, respectively. The HVAC system of the Leaf Lab is composed of ground water source heat pumps with a nominal heating COP of 4.8 and a cooling EER of 6.2–7. The HVAC is coupled to a thermal storage water tank which is heated or cooled using excess PV power. The roof of the Leaf Lab is covered by a PV system of 236.5 kWp. Energy systems are integrated by MyLeaf platform, which allows monitoring of measurements and advanced control functions [48]. Leaf lab is part of the Leaf Community, a microgrid integrating industrial and office buildings with various renewable energy systems (biPVs, tracker PVs, and micro-hydro power system), storage (electrochemical and thermal), and electric vehicles [49]. A complete description of the Leaf Lab and systems installed along with details of the building modeling and validation procedure are available in [50]. The 3D representation of the Leaf Lab simulation model is presented in Figure 3.

**Figure 3.** The Leaf Lab 3D simulation model.

First, an indicative analysis is performed to the initial simulation model outputs comparing to the energy demand measurements to reveal and discuss the weak points of the non-optimized simulation model (SM) and what improvements are obtained from its optimized version (OSM) discussed in Section 2.1. Next, we present the improvements in energy prediction that are obtained by adopting the shape modeling approach (Section 2.2) and Kalman filtering (Section 2.3) in the postprocessing stage. For the model performance assessment, standard statistical indices are incorporated that are widely used for measuring the performance of prediction models. Let us denote by {*X*ˆ*t*}*<sup>T</sup> <sup>t</sup>*=<sup>1</sup> the under study model predictions for the energy demand at certain time instants *<sup>t</sup>* <sup>=</sup> 1, 2, ..., *<sup>T</sup>* and as {*Xt*}*<sup>T</sup> <sup>t</sup>*=<sup>1</sup> the true (measured/observed) energy demand stature at the same time instants. In particular, the following statistical indices are used.


$$NSE = 1 - \frac{\sum\_{t=1}^{T} (X\_t - \hat{X}\_t)^2}{\sum\_{t=1}^{T} (X\_t - \tilde{X}\_t)^2}.$$

taking values on (−∞, 1] where an index equal to zero corresponds to a perfect prediction whereas values below zero corresponds to the case where the prediction model was outperformed by a reference model *X*.

#### *3.1. Indicative Analysis of the Initial Simulation Results: Revealing the Weak Points*

A snapshot of the simulated versus measured HVAC electric power for working days of the week from 7/8/17 to 11/8/17 is presented in Figure 4. Simulated versus actual measured total HVAC electric energy consumption for the year 2017 is presented in Figure 5a. It is noted that there are significant deviations on a monthly basis which become more evident for months of unstable environmental conditions such as March, May, and September, as well as in February. The total annual HVAC electrical energy consumption in 2017 for Leaf Lab was 285,719 kWh (47.61 kWh/m2), whereas the respective simulated value is 197,305 kWh (32.88 kWh/m2). The corresponding values excluding weekends is 232,669 kWh and 180,522 kWh, which correspond to an unacceptably high level of monthly CVRMSE equal to 40.7%.

**Figure 4.** Simulated vs. measured HVAC electric power for the week from 7/8/17 to 11/8/17.

Following parameterization of the model monthly deviation in HVAC, electric energy consumption are minimized to a percentage difference ranging from 1.78% in February to 30.28% in September and an acceptable CVRMSE of 13.89% (Figure 5b). Total simulated HVAC electric energy in this optimized case is 252,689 kWh or 221,946 kWh excluding weekends associated to a percentage difference of 11.5% and 4.6%, respectively.

**Figure 5.** Measured total HVAC electrical energy consumption for the year 2017 versus: (**a**) initial model (baseline) (**left plot**) and (**b**) parameterized model (optimized) (**right plot**).

Table 1 contains monthly measured HVAC electrical energy consumption values and the corresponding simulated ones obtained from the baseline and optimized models for 2017.


**Table 1.** Monthly measured HVAC electrical energy consumption values and the corresponding simulated ones obtained from the baseline and optimized models for 2017.

Similarly, results from the initial simulation in 2018 are presented in Figure 6a. Notably, there are significant levels of underprediction especially in April, May, and September. The total annual HVAC electrical energy consumption in 2018 for Leaf Lab was 290,760 kWh (48.46 kWh/m2), whereas the respective simulated value is 221,136 kWh (36.85 kWh/m2), which equals a percentage difference of 23.9%. The corresponding values excluding weekends is 250,042 kWh and 200,461 kWh, which correspond to an unacceptably high level of monthly CVRMSE equal to 31.42%.

**Figure 6.** Measured total HVAC electrical energy consumption for the year 2018 versus: (**a**) initial model (baseline) (**left plot**) and (**b**) parameterized model (optimized) (**right plot**).

Table 2 contains monthly measured HVAC electrical energy consumption values and the corresponding simulated ones obtained from the baseline and optimized models for 2018.


**Table 2.** Monthly measured HVAC electrical energy consumption values and the corresponding simulated ones obtained from the baseline and optimized models for 2018.

An intensive search of the parameters that affect the daily, monthly, and annual power distribution profiles is presented in the following to improve the initial results of the model based by minimizing deviation from HVAC power consumption data. Through trial and error various combinations and fine-tuning of the all of key parameters are carried out to reach the optimum results both when assessing deviation at intra-day, monthly, and annual timescale. During this approach, the two key parameters identified as critical to decreasing the model's deviation from actual metered energy consumption are the COP and infiltration rates. COP nominal values are initially used but as they are representative of standard conditions, validation is required to better estimate their performance in dynamic conditions. Furthermore, the energy model uses performance curves to simulate the dynamic behavior of heat pump systems which is hard to define in the absence of very elaborate measuring equipment especially for systems customized to provide heating and cooling for large facilities. Besides, fine-tuning model parameters to match the actual performance in such systems is further justified by the fact that user behavior is not recorded and even if it was, modeling of such

a complex activity is not feasible provided the current features of energy building software tools. Indicatively, the opening of external windows, the entry (or exit) of people or industrial equipment and the manual control of the temperature set point in the various thermal zones are some factors of uncertainty as to the actual thermal energy exchanges of a building such as the Leaf Lab.

Following careful fine-tuning of the model, simulation results are improved as shown in Figure 6b. Specifically, parameters related to the infiltration, ventilation, internal loads (electric equipment and lighting), and occupancy of the various zones were varied until the best possible correlation was established. Although some noticeable at monthly level differences remain especially in April and May, overall, the simulated results approach actual measured values as revealed by CVRMSE of 14.33%. In Figure 7, simulated versus measured HVAC electric power for a week in November 2018 is presented. In this case, the simulated pattern follows the actual measured values; however, there clear deviations especially with respect to the peaks early in the morning and late in the afternoon are observed.

**Figure 7.** Simulated versus measured HVAC electric power for the week from 12/11/18 to 17/11/18.

In Figures 8 and 9, the evolution of statistical indices (Bias, MAE, and RMSE) regarding the deviance of the initial simulation model is illustrated, and its optimized parameterized version comparing to the measured energy demand for the years 2017–2018. With regards to Bias, the optimized model significantly outperforms the initial energy model of the building which underestimates energy consumption for most months in 2017. Concerning MAE and RMSE, it is demonstrated that the optimized model generally performs better except for months March, April, October, and November for which the initial model performs slightly better. Similarly, in 2018, the results of the optimized model are associated with a Bias error of lower magnitude for months from April to October. MAE and RMSE, in this case, are lower in the case of the optimized model compared to the baseline model for months from February to November. In contrast, the initial model performs marginally better for estimating energy consumption in January, October, and November.

**Figure 8.** Statistical error indices for the energy prediction through the simulation model (SM) and the optimized simulation model (OSM) for the year 2017 (per month).

**Figure 9.** Statistical error indices for the energy prediction through the simulation model (SM) and the optimized simulation model (OSM) for the year 2018 (per month).

#### *3.2. Diagnostic Results for the Complete Model Outputs*

In this section, we discuss the improvements obtained to the energy consumption prediction for the Leaf Lab building for the time period 2017–2018 by implementing the postprocess part of the integrated model, i.e., the shape modeling approach and Kalman filtering. For convenience, let us introduce the abbreviations (SM) for the initial (benchmark) simulation model, (OSM) for the optimized/parameterized simulation model, (RS) for the reshaped simulation model, (w-RS) for the weighted reshaped model, and (KF-RS) for the integrated model by implementing Kalman filtering at the later stage of the reshape model output. For illustration purposes, we divide our analysis into two parts: (a) regarding the intra-day (hourly) energy demand prediction and (b) regarding the prediction of daily summary of the energy demand (total energy demand on the whole day).

#### 3.2.1. Evaluating Intra-Day Energy Demand Predictions

In Table 3, the yearly diagnostic results of the intra-day energy demand predictions (hourly), provided by all incorporated prediction models for the years 2017 and 2018, are presented. Note that the NSE coefficient is calculated using as reference model the initial simulation model (SM). Based on the various indices, it is clearly demonstrated that the reshape modeling provides much higher levels of accuracy compared to the model's simulated outputs. Specifically, it is observed that RMSE in both years is reduced ~50%. The Kalman filtering procedure does not further improve the intra-day prediction performance of the reshaped model, indicating that the utilized shape model filters out efficiently most systematic sources of error affecting the intra-day (hourly) energy demand prediction.

**Table 3.** Model diagnostic results for intra-day energy demand predictions (hourly) for the years 2017 and 2018.


Figures 10 and 11 illustrate the monthly Bias, MAE, and RMSE values for both reshaped models RS and w-RS computed for the intra-day performance of the models compared to the results obtained by the optimized simulation model (OSM). Both reshaped models perform better as all statistical indices in the majority of cases are significantly improved indicating the superior predictability of the reshaped approach.

**Figure 10.** Statistical error indices for the Reshaped Simulation Model (RS) and the Weighted Reshaped Simulation Model (w-RS) for the year 2017 (per month).

**Figure 11.** Statistical error indices for the Reshaped Simulation Model (RS) and the Weighted Reshaped Simulation Model (w-RS) for the year 2018 (per month).

#### 3.2.2. Evaluating Daily Summaries Energy Demand Predictions

At a second stage, we are interested in the performance of the incorporated models in predicting the total energy demand per day. The same model approaches are tested, and their results are illustrated in Table 4. In this case, the Kalman filtering step further improves the prediction since the model biases are reduced significantly and the best NSE values are obtained through the KF-RS model. This happens because although the RS and w-RS models efficiently capture intra-day effects in the shape of the energy demand, when we are interested in the daily summary, these effects are mutually canceled and become obsolete. However, the KF-RS model, in general, improves the average biases, as it detects these inefficiencies as systematic ones from the perspective of daily summaries, and therefore it further improves the model outputs.

**Table 4.** Model diagnostic results for the predictions on the daily energy demand (summaries) for the years 2017 and 2018.


In Figures 12 and 13, the performance of the KF-RS model in predicting the daily energy demand (summaries) is illustrated, comparing again to the OSM model. It is evident that the later model is significantly improved by correcting major inconsistencies in certain time periods throughout the year and more accurately approximating the true levels of energy needed. In particular, for special cases where systematic biases are present (see, e.g., highlighted periods in Figures 12 and 13), the KF-RS model proves able to eliminate systematic over- or underestimations.

**Figure 12.** Prediction of the daily energy demand (summary) obtained by optimized simulation model (OSM) and Kalman filter-enhanced reshaped simulation model (KF-RS) for the year 2017.

**Figure 13.** Prediction of the daily energy demand (summary) obtained by optimized simulation model (OSM) and Kalman filter-enhanced reshaped simulation model (KF-RS) for the year 2018.

#### **4. Conclusions**

In this paper, validation of the building energy model of an industrial near-zero energy building is investigated over a two-year period. It is demonstrated that manual extensive fine-tuning of key model parameters is valuable to improve initial overall error levels using trial and error of key parameters. However, this process requires a high level of expertise, deep knowledge of the facility under study, and it is very time-consuming. Furthermore, it is demonstrated that systematic deficiencies continue to occur even after careful fine-tuning of key model parameters. On the one hand, this is due to the complicated issue of modeling the behavior of users as well as of equipment use in a large industrial facility. On the other hand, the dynamic performance of custom-built HVAC systems requires specific measurements if it is not to be solely defined at a high level and modeled based on nominal performance coefficients. To address this problem two modeling approaches are integrated: (a) the shape and weighted shape model and (b) Kalman filtering postprocessing. The above methods are applied in a postprocessing stage, which is tested and evaluated. Initial, optimized simulation results and results from postprocessing are analyzed compared with the aid of bias error, mean absolute error, root mean squared error, and the Nash–Sutcliffe model efficiency coefficient. Results are explored to demonstrate the effectiveness of the method in capturing and reducing intra-day systematic deviations as well as the overall performance gap. Overall, the proposed integrated approach proves to be very effective in eliminating systematic over or under estimations and critically reduces the magnitude of deviations and as a result significantly reduces the performance gap compared to the optimized simulation model. In particular, the integrated prediction model succeeded in reducing the RMSE (in mean values) approximately 39% and 43% for the hourly predictions in years 2017 and 2018, respectively, and approximately 32% and 25% for the daily energy demand for the same years which are rather impressive improvements. Note also that in all cases, the integrated prediction model substantially increased model efficiency coefficient (NS) at least 70% comparing to the efficiency of the

optimized simulation model (OSM). The proposed approach is applicable in several types of buildings (i.e., residential, commercial, public, etc.) provided that reliable data of energy consumption and of human activity are available over a significant period of time (ideally more than one year). As a future step, the presented approach could be deployed and tested across several types of buildings that fall into different energy efficiency categories to assess its performance over a wide spectrum of applications.

**Author Contributions:** Conceptualization, N.K., G.I.P., D.K., G.N.G. and A.N.Y.; Data curation, N.K., G.I.P. and D.I.; Formal analysis, N.K., G.I.P., G.N.G. and A.N.Y.; Funding acquisition, N.K., D.K. and C.C.; Investigation, N.K., G.I.P., D.K., G.N.G. and A.N.Y.; Methodology, N.K., G.I.P., D.K., G.N.G. and A.N.Y.; Project administration, C.C.; Resources, D.K., G.N.G., D.I., C.C. and A.N.Y.; Supervision, D.K. and C.C.; Validation, G.I.P.; Visualization, N.K. and G.I.P.; Writing – original draft, N.K., G.I.P., G.N.G. and A.N.Y.; Writing – review & editing, N.K., G.I.P., D.K., G.N.G., D.I. and A.N.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skeodowska-Curie grant agreement No 645677.

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

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
