Demand-Side Flexibility Impact on Prosumer Energy System Planning

: Apart from numerous technical challenges, the transition towards a carbon-neutral energy supply is greatly hindered by limited economic feasibility of renewable energy sources. This results in their slow and bounded penetration in both commercial and residential sectors that are responsible for over 40% of ﬁnal energy consumption. This paper aims to demonstrate that combined application of sophisticated planning methodologies at building-level and presents incentive mechanisms for renewables that can result in prosumers, featuring hybrid renewable energy systems (HRES), with economic performance comparable to that of conventional energy systems. The presented research enhances existing planning methodologies by integrating appliance-level demand side management into the decision process and investigates its effect on the planning problem. Moreover, the proposed methodology features an innovative and holistic approach that simultaneously assess electrical and thermal domain in both an isolated and grid-connected context. The analyzed hybrid system consists of solar photovoltaic, wind turbine and battery with thermal supply featuring solar thermal collector and a ground-source heat pump. Overall optimization problem is modeled as a mixed-integer linear program, while ranking of all feasible alternatives is made by the multicriteria decision-making algorithm against several technological, economic, and environmental criteria. A real-life scenario of energy system retroﬁt for a building in the United Kingdom was employed to demonstrate overall cost savings of 12% in the present market and regulation context.


Introduction
Current trends and existing policies related to energy supply aim to alleviate global dependency on fossil fuels due to their finite nature and proven devastating effect on the environment. This led to the development of new technologies ranging from alternative and sustainable energy sources to innovative management solutions for available energy resources. Considering exploitation of renewable energy resources as the most prominent of them, many countries adopted suitable regulation and incentive mechanisms to improve the status quo and increase share of renewables in the overall energy mix. For example, the previous target set by the European Union (EU), which envisaged at least 20% of renewable energy sources in the final energy consumption by 2020, was recently updated to reach at least 32% share of renewables by 2030 (raised from 28% set previously) [1].
As part of this transition and in line with contemporary directives [2], buildings are increasingly being equipped with local energy sources and granted with capability to conduct individual energy management strategies, yielding a new entity in energy systems referred to as prosumer. However, ensuring a widespread transition from consumers towards prosumers depends on the economic viability of renewable energy sources and accompanying equipment. Although renewable energy sources are already a relatively mature technology, prosumers are still facing significant issues regarding upfront investments, making it hard to level up with the cost-effectiveness of conventional energy supply from coal and gas. Aiming to overcome this issue and spark the necessary widespread adoption of renewable technologies, national governments implemented various incentive programs and mechanisms such as feed-in or generation tariffs, carbon credits, tax refunds, and procurement subsidies tailored for each type of renewable technology. All these measures share a common goal to alleviate the economic impact of associated capital and operational expenditures. In doing so, they intend to make the overall investment profitable on the long run and render the transition from conventional to renewable sources viable for consideration. Apart from instantaneous effects to present investments, these incentive programs also contribute to a decrease of manufacturing costs of renewable energy technologies in the long-term by scaling up the overall number of installations. Besides various financial incentive mechanisms, economic performance of prosumers can also be improved with application of innovative energy management technologies. They not only increase operational efficiency of renewables but also contribute to increased reliability and availability of energy supply in the context of multisource hybrid renewable energy systems (HRES). Given the intermittent nature of renewable energy sources, the latter becomes increasingly important as the transition to renewables in urban areas also needs to meet very high availability and quality of energy supply standards. Therefore, to establish a cost-effective and reliable prosumer, one would need to asses two fundamental aspects: (a) The planning/dimensioning problem, which represents an optimal rated power split of different renewable sources and storage capacities within prosumer systems, resulting from a multicriteria decision making process against complex objectives combining maximization of economic performance, environmental neutrality, and independence from the power grid. (b) The operation problem, which focuses on optimal energy management strategy for a given prosumer and its energy assets. Moreover, it considers optimization of prosumer's energy imports and exports as well as internal power flows between multiple renewable/conventional energy sources and storages against multiple technological, economic, and environmental criteria.
As reviewed by [3], there are many commercial software solutions that aim to solve different individual aspects of these two problems while offering different types of outputs (high-level efficiency or economical parameters, dynamic operational values, etc.). However, although seemingly independent from each other, the two are intrinsically correlated as the decision on optimal prosumer energy assets is impacted by their day-to-day management and, ultimately, reached through a long-term simulation of its operation.
The considered planning problem was extensively investigated starting from standalone HRES in isolated rural areas, where there was a lack of conventional energy supply to those grid-connected, as the penetration of local energy sources in urban areas became more significant. Following is a brief overview of more recent research efforts dealing with such planning problems and related optimization approaches. An approach for both stand-alone and grid-connected modes using the energy filter was discussed in [4], while [5] proposes a multicriteria decision analysis for PV-WT grid connected systems. A research effort conducted on HRES in the form of a microgrid system in [6] incorporates the addition of a battery energy storage systems (BESS) and employs two procedures, a source sizing and battery sizing algorithms in sequence. Scalfati et al. [7] proposes a mixed-integer linear programming (MILP) based solution for sizing that can, in its general form, be used for different microgrid architectures and storage technologies. Sizing with sensitivity of a microgrid structure specific for a university campus is discussed in [8], while optimal sizing with implemented DR strategies is discussed in [9] where a HRES configuration is considered. HRES optimizations with regards to energy, economic, and environmental indicators can also be found in [10] where a multiobjective optimization was employed to determine the best system configuration. The authors note that no single system configuration can simultaneously satisfy the selected three criteria and proposes the selection of a trade-off Pareto optimal solution. Economic parameters and capital investment benefit analyses can also be found in [11] where optimal sizing and power man-agement of prosumers equipped with PV was considered using a two-step approach, with the first step analyzing the technical model (short-term assessment resulting in outputs such as component sizing and battery lifetime), and the second one tackling the long-term assessment through economic modeling. The authors report that an increase of profitability by up to 14% was achieved. Finally, a stochastic approach using MILP for determining optimal sizes of prosumer assets (PV generation capacities and batteries) is proposed in [12]. The authors describe a methodology that minimizes the joint combination of investment, maintenance, and operational costs in different scenarios that result in varying energy consumption levels from the grid.
Regarding the optimization problem, various different techniques can be found in the related literature. As the proposed methodology includes a proactive approach based on utilization of demand-side management (DSM), and in line with findings from a systematic review from [13], the most prominent technique for this type of problem is linear programming (LP) with its mixed-integer variant being the one most commonly used. For example, ref. [14] proposes a multiobjective mixed integer linear programming (MOLIP) technique to facilitate residential DSM in a system where effects of storage systems are specifically analyzed. Also, ref. [15] formulates a MILP model to be used for optimizing profit on the electricity market of a system with photovoltaic (PV) panels and BESS. These optimizations are performed for a horizon of 24 h with hourly varying prices. This model does not consider load to be appliance-based but rather views it as an aggregate value. Also, since the simulation is performed for a short amount of time, the monetary investments and maintenance costs of running such a system are not considered. Paper [16] continues with a MILP model also employed for a 24 h simulation horizon but with a shorter, 15 min-long sample period. The modeled system considers optimal appliance scheduling with a PV source present, with load scheduled on a per appliance basis. The results were obtained and discussed for both single and multiuser scenarios. However, the investment and maintenance costs related to running a renewable energy source are also not considered. The framework laid out in [17] also implements a MILP model for optimal appliance scheduling during a 24 h horizon with a 15 min sampling period. The chain rule, defining that a given appliance can only be started after another one finishes its operation, is introduced. The model output is discussed for three scenarios in a specific use case: a fixed price tariff, a variable price tariff with ripple control (devices that switch on or off appliances based on the current tariff), and a variable price tariff with optimal scheduling. Concluding that, because of the insignificant difference in the applied price tariffs, it would not be viable for an average domestic consumer to look for a solution more sophisticated than ripple control. The authors also state integrating distributed generation and storage into the model as a future research point. Finally, ref. [18] focuses on optimizing energy management of a residential microgrid with the goal of analyzing the relation between the level of demand flexibility and cost savings. This paper also models investment, maintenance and replacement costs of BESS as well as distributed PV and wind turbine (WT) generators, and it introduces a discreetly operated appliance whose operation can be split into multiple nonconsecutive time periods. Considering that this system uses a relatively lengthy, one-year-long horizon with a sample rate of 1 h, the model is designed on an efficient window-based concept. Nonetheless, this approach does not considered load to be appliance-based, i.e., specific appliance activations cannot be traced in the final results. Also, a notable addition with this paper is that the sizing problem is solved simultaneously with scheduling using the appropriate variables implemented in the model. Although this is a very efficient way to solve such a problem, the linear programming paradigm constrains those variables in a linear way, and thus, limits the modeling potential to a certain extent.
Finally, building upon previous research and aiming to improve the current state of the art, this paper proposes the introduction of a two-step optimization process that jointly considers the planning and operation problems. As will be described in greater detail in the following section, the inner loop optimizes operation of individual configurations using the available load shifting mechanisms while the outer loop explores key performance indicators for each of these configurations and selects the one that adheres best to the desired preferences. The remaining part of the paper is organized as follows. Section 2 presents the proposed planning methodology and description of the optimization process as a whole. Section 3 unveils the mathematical implementation of both operation and sizing models supported by Appendix A which outlines the models used for simulation of renewable technologies (RET). Section 4 describes the employed real-life use case scenario for methodology testing and verification, while Section 5 discusses the obtained results depending on the desired criteria. Lastly, Section 6 provides concluding remarks and conclusions. The paper is also supported by a nomenclature table in Appendix B that contains the list of all abbreviations and variables used.

Proposed Planning Methodology
The overarching objective of this paper is to demonstrate that by combining sophisticated methodologies for planning and operation of prosumers leveraging small-scale residential HRES and existing governmental incentives for renewable technology, one can retrofit or even completely replace conventional energy supply without making a compromise regarding economic viability of such transition. In other words, the presented research aims to mitigate the economic barrier associated with global uptake of renewables by demonstrating cost-effectiveness comparable to that of conventional energy sources. Moreover, the considered real-life use case scenario exhibits the possibility to even exceed it in the context of existing long-term incentive programs and country-specific energy regulations.
To reach this objective, a novel planning methodology for future prosumers was established by leveraging the benefits of increasingly utilized mechanisms of DSM. In particular, the underlying approach exploits load elasticity, both in time and intensity, to reduce capital investment in renewable energy technologies, improve their economic performance and increasing overall penetration in energy supply portfolio. In addition, the proposed decision making process simultaneously assesses multiple consumer-defined criteria presented in Sections 2 and 5. The techno-economic performance of viable configurations is discussed in the latter section, where the effects of the proposed planning methodology are evaluated in the most conspicuous way. Moreover, the proposed methodology assumes a holistic approach for the planning process, which simultaneously considers both electrical and thermal domains. Although previous research mainly assessed these two domains separately, they are inevitably cross-correlated, especially in cases when thermal demand is satisfied through a heat-pump, which combines any available thermal source (e.g., ground, solar, or air) and a proportional amount of electricity, as described in Appendix A. Moreover, such consideration is even more relevant in cases where thermal demand is satisfied from several different sources, (e.g., gas, electricity, or solid fuels). Following a theoretical elaboration, the proposed methodology is demonstrated through its practical application in a real-life scenario featuring actual technical, economic, and environmental constraints.
The proposed HRES planning methodology, developed to devise an optimal system topology as well as sizing of its individual components, aims at fundamentally enhancing existing planning tools and algorithms by combining different approaches and adding new design aspects. In short, it introduces and brings together the following aspects: 1st. The overall HRES planning process considers simultaneously both electric and thermal energy demand, while current approaches typically consider electric or thermal domain, exclusively, with the methods for such optimizations previously discussed in [19]. Employed methodologies therein are focused on balancing the selected demand type with available energy sources, conversion elements, and storages. However, increased utilization of devices like heat pumps, which satisfy the thermal demand while contributing to electricity demand, requires a holistic assessment approach. The differences between the traditional approach and the one proposed by this paper is illustrated in Figure 1. 2nd. The most utilized approach to consider isolated (island) HRES deployment scenarios is extended towards consideration of grid-tied deployment, which brings a more dynamic context where varying import and export energy prices are applied and unlimited energy exchange with power grid is enabled.
3rd. The increasing application of DSM programs and, more specifically, DR schemes in day-to-day operation is considered on an appliance level and corresponding implications on the planning of HRES and dimensioning of individual components are evaluated.
4th. MCDMA is employed to rank feasible HRES topologies with capabilities of simultaneously evaluating a wide range of technical, economic, environmental, and societal design criteria.
In the following elaboration, each design alternative will be referred to as HRES configuration, which is defined with a set of discrete sizes for each RET components. Hence, the configuration will consider both the HRES topology and sizing of the energy assets within.  The optimization process can be split into use case evaluation and two main distinct sections: operation optimization and sizing optimization, as presented in Figure 2. Firstly, in the evaluation stage, the values like demand profiles, financial information and RET parameters are collected and fed into the model. A set of predefined HRES configurations deemed fit for the selected use case is defined, and when it comes to the definition of the search space for the optimal HRES configuration, a set of context-defined and user-defined constraints is established. The following list summarizes the most influential design constraints into several categories, which are simultaneously assessed by the proposed methodology to deliver optimal HRES topology and sizing: The listed constraints, in fact, define a set of boundaries for the space in which the optimal design solution is searched for. The operation optimization stage is initiated by equipping the model iteratively with one of the predefined alternatives.

Operation Optimization
The underlining structure of the system that was selected for operation optimization and is implemented in this study is a solution proposed in [20], referred to as the Energy Hub, with its general structure illustrated in Figure 3. It represents how the input power vector (P in ) is transformed in several stages (denoted by matrices F −1 in for input transformation, C for conversion and F out for output transformation, as discussed in Section 3.1) while also managing storage charge/discharge rates (Q in for the input storage stage and Q out for the output storage stage), exported power (P exp ) and loads (L). Figure 3. General structure of Energy Hub.
Given the variability of generation from renewable sources, energy demand requirements, and dynamically determined energy costs, each feasible configuration is modeled and its long-term operation is simulated and optimized. When deciding on the simulation time horizon and size of the time step resolution, several factors were taken into account. Firstly, long-term performance of an HRES is heavily influenced by the renewable energy harvesting potential, which implies utilization of historical meteorological conditions for a given location. Moreover, since the performance should be evaluated for the entire HRES lifetime, which is typically around 20 years, the so-called typical meteorological year (TMY) data, which provides hourly data for typical meteorological conditions over a relatively long period and, thus, objectively characterize the specific geographical location, was selected. Secondly, given the objective of long-term assessment of HRES operation, the linear (or mixed integer) models of energy assets running on hourly resolution are sufficiently complex. Consequently, hourly resolution was chosen as the time step while the time horizon is one-year-long. However, long-term performances can be extrapolated by multiplying these results as many times as needed. This is due to the methodology natively taking into account events such as replacement of an asset (e.g., batteries typically last for 5 years) or time-limited governmental incentives (e.g., payment periods of 7-20 years for RET subsidies).
The operation aspect of the described model is finished once all the predefined HRES sizing configurations are optimized and their internal variables saved for later processing in the sizing segment.

Sizing Optimization
Finally, when the iterative simulation procedure of the inner loop is concluded, and the list of evaluation criteria is derived for each configuration, it is necessary to establish appropriate ranking among the configurations, considering multicriteria (multidimensional) evaluation space. To do so, an existing MCDMA algorithm was employed for this segment, previously called the outer loop.
At the initial stage, the MCDMA is fed a list of alternatives A, which ought to be simultaneously ranked across multiple criteria C. Each alternative A, represented by an individual HRES configuration, is first individually evaluated across the whole range of criteria C, referred to as evaluation criteria in the dynamic performance assessment elaboration. Following is the assignment of the weighting factors w to each criterion C. These weight factors allow for non-uniform distribution of the level of importance assigned to each criterion, allowing end user to practically steer the selection process according to desired needs and preferences. Following the findings from [21] noting its common applications for technology evaluation, the acknowledged performance ranking organization method for enrichment evaluation (Promethee) II algorithm [22] was utilized for the purpose of development of MCDMA functionality, as described below.
Firstly, a comprehensive pair-wise comparison is calculated as Afterwards, those differences are passed through a preference degree function π k (a i , a j ) where P k can be defined in a variety of ways, with the most common being a linear variant bounded by q k and p k . Every pair of actions is compared using a multicriteria preference degree where a constant to weights is applied After calculating Finally, ranking all alternatives according to φ(a) gives a complete ranking considered as the output of the Promethee II algorithm. The proposed methodology can easily consider different design criteria, which can be summarized into four general categories. In the following list, these categories are out-lined with detailed elaboration of each category that was implemented and its items following later in the text:

Mathematical Model
The key feature of the proposed methodology is the model implemented to facilitate the optimization process and, just like the process itself, the mathematical model can be represented with two sets of parameters, one depicting operation and one depicting sizing optimization.

Operation Optimization
Generally, a MILP problem is defined as determining a vector of variables that minimizes a certain objective function f which is usually written as Meanwhile, the vector x opt must also adhere to a set of conditions that are split into four categories: equality and inequality constraints, lower and upper bounds and integer constraints. If a subvector x int of x is defined as and the lower bound l b and upper bound u b vector as these constraints can be written as The vector of variables x is formed by arranging a set of all variables required for the model at every instance of the simulated time horizon like If T s is used to denote the sample rate of the model, for any of the variables V f in x, V f (k) may denote the value of V f at the k-th time sample, i.e., V f (k) = V f (kT s ). Since some of these variables should incorporate values for different types of electric energy sources (WT, PV, and electricity from the grid), they are natively multidimensional. However, for computational modeling, x must be a row vector, hence these natively multidimensional variables are flattened out in a predefined order. For example, given a simulation time horizon of T = NT s , the variable P f in (k) holds instantaneous imported power values for WT if 0 < kT s ≤ T, values for PV if T < kT s ≤ 2T and values for grid electricity if 2T < kT s ≤ 3T. Therefore, for modeling purposes, it may be convenient to think of such variables as matrices in a reshaped form. Regarding our example variable P f in (k), we may consider a three-by-one native vector P in (:, k) where P in (1, k) is the corresponding instantaneous imported power from the WT, P in (2, k) is the corresponding instantaneous imported power from the PV array and P in (3, k) is the corresponding instantaneous imported power from the grid, all calculated at t = kT s . Obviously, every matrix equation that employs such variables in their native shape can be easily converted to a vector based expression required for MILP implementation. Thus, for the sake of clarity, variables V f from the vector x will, in the remainder of this paper, usually be written in their native form V(:, k) abbreviated as V(k) with an index k referring to the instance of time at which the variable is being evaluated. For the sake of clarity, the index k is omitted when a variable is referenced in text, but is present in all formulas where that variable is used.

Energy Balance
The instantaneous imported power can be from either the WT, PV array, or the grid and this power can either be stored at the input level or dispatched to the rest of the system through the appropriate transformers. The law of conservation of input power states that the balance must hold. The power sent to the storage unit is converted into energy via The available energy of the storage unit is determined by an integral expression with an initial condition as E in (1) = E in1 defining the SOC in the first sample, and is usually set to zero if optimizations on consecutive time intervals are not being performed on the same system. The energy not being stored at the input is sent to the conversion stage defined as The output of the conversion stage can then either be exported back to the grid or further dispatched to the output stage. This is simply written as with exporting power which was previously imported from the grid back to said grid being directly prohibited by (∀k) D exp P exp (k) = R exp (17) where D exp is a matrix that determines which carrier is to have a fixed (restricted) export and R exp sets those values. The power not being exported is then sent to the output transformation stage that aggregates the carriers into an arbitrary number of values depending on the number of load types. This operation is performed in the equation Analogous to the input, the output can also feature a storage option. The charge/discharge rate is obtained from Similarly to the input stage, the output storage SOC is calculated by and an initial condition is given as E out (1) = E out1 which concludes the set of equations governing the energy balance of the energy hub system. Furthermore, additional equations must be added for load management mechanisms.

DSM and Load Variables
To properly allow the model to perform necessary optimizations, some auxiliary features must be introduced as constraints. First of all, load L can either be attributed as fixed load L fix which cannot be optimized in any way or flexible load L flex which can used for optimizing by means of time shifting (shiftable load), splitting its operation into multiple time instances (dispersible load) and adjusting its instantaneous power consumption (elastic load). Nevertheless, for every appliance i, we define an on/off state variable that is defined by The flexible load at the k-th time sample can be written as with the total load being expressed as However, the product being summed in (22) would represents a nonlinear operation between two variables and so it must be rewritten. To achieve this, P i is divided into three components: nominal power draw, positive power deviation, and negative power deviation from the nominal value, or in other words Nevertheless, (24) also incorporates a product between variables, however d + i and d − i will be constrained to having nonzero values only when the appliance is turned on, this expression can be reduced to Finally, total load can be expressed as Since P nom i is set beforehand, this expression is actually a linear combination of variables (subvectors) from x , and can therefore be implemented as a MILP constraint. According to the classification laid out in [23], elastic loads can be classified either as being either energy-based meaning that they must consume a predefined amount of energy within a specified time window or comfort-based meaning that they must control an environmental variable within a desired range. With effects of DSM on comfort-based appliance being previously investigated in [19], this paper will focus only on energy-based elastic appliances for DSM with an option to elastically adjust their power within given bounds. Therefore, for each appliance, a set of windows (activation cycles) is defined with one of them, a vector 0, k is not in the window n, k is in the window (27) defining the n-th window of i-th appliance by having nonzero values equal to n at time instances that belong to that window. This implementation allows for windows that need not be continuous, i.e., they can be split into an arbitrary number of segments. Since these windows are also predefined like the nominal power, they can be used for forming an energy constraint stating that a specified appliance i must only be active a given amount of times so that the amount of energy it spends during that activation cycle n is equal to the product between nominal power P nom i and the length ∆t (n) i of nominal activation belonging to that window. Because power deviations also affect the energy consumption, it is also stated that the sum of power deviations must be equal to zero during a given window, or in other words thus finalizing the set of equality constraints required for the model. Nonetheless, these relations are not sufficient for the model and some additional conditions must be applied in form of inequalities.

Auxiliary Constraints
Because a set appliance should not have a nonzero positive and negative power deviation at the same time, two variables are introduced to indicate when these deviations are active by defining with the constraint prohibiting simultaneous positive and negative nonzero deviations. Additionally, defining another constraint in combination with (31) forces the indicators to have a correct sign. Finally, by specifying a link between the deviations and their respective indicators, and thus, forcing these variables to uphold the definition set by (30) was provided. To allow for an easier implementation of the results obtained from this model and to facilitate load dispersion penalization in the objective function, another variable called the device start indicator is introduced by marking that the device i will be turned on in the following time sample. Since this definition describes a nonlinear relation between y i and z i , the proper values for z i are obtained by simultaneously enforcing three inequality constraints This same indicator is also used to allow for modeling both dispersible and nondispersible appliances by specifying and concluding the primary set of equalities and inequalities for the model.

Variable Bounds
To supplement the equality and inequality constraints stated previously, a set of bounds in posed for the variable vector x. All of the instantaneous values of power must be non-negative, and thus, the following bound is imposed On the other hand, the imported power P in must be equal to the available amount provided by the renewable sources P renew at all times if considering those respective elements, or unbound if considering elements describing the import from the grid. Therefore, At both the input and output stages, storage levels must be between the lowest possible (zero) and highest possible (battery capacity) SOC and so with Q in and Q out being limited by where Q max is the highest achievable charge rate, and thus, also bounding q in and q out . The total load L only has a defined lower bound equal to the value of fixed load since the flexible load is non-negative, and thus (∀k)(L fix (k) ≤ L(k)).
As mentioned before, both deviations d + i (k) and d − i (k) also have bounds equal to a predefined upper and lower deviation limit, respectively, applied during specified windows as follows As for the indicator variables, the device start z i has a lower bound of zero and upper bound of one for all time samples, i.e., while the on/off state y i has the same bound within windows and both bounds set to zero outside A similar logic is employed as to limit the deviation indicators Finally, since the indicator variables should only assume a value of either zero or one, thus rendering this problem to be classified as MILP rather than LP, we specify

Objective Function
Depending on the effect that is desired to be achieved, different objective functions can be formed. Relevant literature most commonly considers cost minimization, discomfort minimization, and maximization of on-site generation use. These criteria can also be combined as to create a mixed objective and so one such possibility is considered

Cost Minimization
One of the most frequently considered parameters when discussing feasibility of renewable generation systems is the monetary cost that ultimately falls on the end consumer. To model the effects that running a HRES has over the simulated horizon, the prices of energy attributed to imports and exports is taken into consideration. The active cost of running such a system can be calculated with Factors α WT and α PV usually represent zero or negative values because the use of renewable generation is generally subsidized by governments and their values, like the values of other factors in (47), vary depending on local regulations and acting price tariffs. Therefore, the parameters of such a cost function are use case dependent, and their exact values will be discussed later on in the paper.

Dispersion Minimization
Sometimes, splitting one appliance's operation cycle into several disjointed segments may be considered as unwanted behavior impacting user's comfort. To combat this behavior for dispersible appliances, a criterion is defined as Minimizing such a function by itself would lead to no appliance utilizing dispersion, hence, a combined criterion J cd = J c + J d may be used and is implemented in the proposed solution to simultaneously balance minimizing cost and penalizing unwanted load dispersions.

Sizing Optimization
As mentioned previously, the second part of the proposed methodology, besides the optimal management of energy resources accomplished by the MILP solver, is determining the proper configuration of the site, also referred to as the sizing problem. Since multiple combinations of available renewable generators and storage options are being considered, a set of criteria is selected to facilitate MCDMA ranking of all available configurations with the first one being optimal with regards to the specified weights associated with each of the given criteria.

Total Cost (EMI)
Since the model application considered in this paper will mainly focus on residential users, the total cost of running a renewable project is one of the most important factors to be considered when ranking different configurations. The yearly cost relating to importing and exporting energy is expressed by (47). However, this relation does not take into account initial investment costs, maintenance and eventual replacement costs associated with using equipment with a finite life span that may malfunction during its operation. To assess all of these expenses in a comprehensive way, equated monthly installments and maintenance costs are calculated. For a piece of equipment labeled i one month's equated installment (equivalent to rent) would be equal to where δ = 0.42%. The corresponding criterion for the MCDMA can now be written as with E = {BESS, WT, PV, STC, GSHP} symbolically denoting all of the devices for which the costs need to be included.

Net-Zero Energy Building Rating
The concept of net-zero energy buildings (NZEBs) gained a lot of prominence lately and its importance was also recognized by several governments with the EU's Energy Performance of Buildings Directive (EPBD) even requiring all new buildings from 2021 to be at least near-zero energy (nZEBs) [2,24]. A NZEB is defined as an energy efficient, self-sufficient structure which roughly consume the same amount of energy as it produces on site over the course of a year. With around 40% of global energy consumption being attributed just to buildings, such a concept of energy management and planning is set to greatly improve the current environmental effects of powering buildings in the near future. To facilitate a balance necessary for a NZEB or nZEB rating, we may consider minimizing the criterion where P out (k) is the amount of power consumed or stored and P in (1, k) + P in (2, k) is the amount of power generated on site by either the WT or PV for each time step.

CO 2 Emissions
Finally, to quantify the impact that running the operation of the considered system has on the environment, the total mass of CO 2 may be considered. This criterion can be simply calculated as

Use Case
To display the capabilities of the proposed methodology, a concrete use case is assumed in form of a simulated residential household property located on the Ravenscliff Road in Motherwell, Scotland, United Kingdom (GPS coordinates 55.80, −3.96) shown in Figure 4, located around 25 km from the center of Glasgow.

Energy Hub Model
The considered system, with its appropriate Energy Hub structure presented in Figure 5, has three different types of input electric energy: WT generated, PV-generated, and grid-imported. No thermal energy is taken into account because the thermal load is considered to be met by the STC and GSHP and is observable trough the fixed electric load generated by these systems as will be described later. Also, the raw imported energy cannot be stored in any way, and therefore the appropriate matrices required for input storage modeling are The imported energy is directly dispatched to the conversion stage so Since the user only pays for the active power measured by the local meter, as not to interfere with the objective function calculation, η grid = 100% is assumed while some losses may incur when converting power from renewables, and thus, η WT = η PV = 95% was set, which is a common inverter efficiency value. Moving on, users with distributed generation such are WT and PV are allowed to export excess energy back to the grid. However, exporting electricity imported from the grid is prohibited with performs an aggregation of all three types of energy carriers at the output stage as only electric loads exist in the model. The output stage features a single storage unit for electric energy, and therefore S out = 1 and S qout = 1.

Tariff Information
The address of the proposed use case is located in an area with the supply code 18 named "South Scotland" meaning that its energy is supplied by ScottishPower. Although several tariffs are provided by ScottishPower with varying prices between them, those differences are mostly negligible, and Online Fixed Saver December 2019 [25] tariff was selected with the consumption costs at 19.63 cEUR/kWh and 4.68 cEUR/kWh for electricity and gas respectively, and daily standing charges at 21.87 cEUR/d for both energy carriers.
Residential users in Great Britain can take advantage of a Feed-in tariff (FIT) scheme administered by the Office of Gas and Electricity Markets (OFGEM). Although feed-in tariffs usually refer to an amount paid just for exporting energy back to the grid, a generation tariff is also offered. With related numbers changing every three months, the generation tariff scheme available for new applicants at the time of writing this paper is 4.40 cEUR/kWh for accredited PV installations smaller than 10 kW and 9.47 cEUR/kWh for WT installations smaller than 100 kW On the other hand, the export tariff equals 5.97 cEUR/kWh for both the PV and WT. The generation tariffs are only paid out during the 20 years following the contract signage.
Furthermore, OFGEM also administers the Domestic Renewable Heat Incentive (RHI) [26] that pays users that are running environmentally friendly heating and cooling solutions like the STC and GSHP. Just like the FIT, RHI is also time limited, with payments limited to 7 years. The current rates are 22.64 cEUR/kWh for heat extracted from the STC and 23.32 cEUR/kWh from the GSHP, but for the sake of simplicity, both are assumed to equal the lower value. Because of these time limitations, it would not be fair for a model with a one-year-long horizon to assume that the RHI payment would be paid out in full for that year, and then take the output of that model as a metric when discussing long-term feasibility. Since the future of similar renewable programs is currently unknown and hard to predict, having in mind the fact that the acquisition costs of renewable sources are sure to decline when eventual replacements are necessary, the model is modified to scale the yearly FIT and RHI incentive payments by the ratio of available payment years and the estimated life span of the device (20/20 years in the case of WT and PV and 7/20 years in the case of GSHP and STC), thus giving a more fair representation of benefits that can be obtained by enrolling in the mentioned programs.

Baseline Consumption
According to the British Department for Business, Energy & Industrial Strategy's (BEIS) yearly statistical report titled "Energy Consumption in the UK (ECUK)" [27], the average household in the study spent around 3903 kWh of electric energy, 10,930kWh of thermal energy for space heating and 3017 kWh for domestic hot water (DHW) per year. Since the use case location is in the northern half of the country where colder climate is prevalent, the baseline thermal demand for the model is selected to equal 120% of the national average, meaning Q heat = 13,170 kWh and Q DHW = 3620 kWh were set as baseline loads whilst the total cooling demand is equal to zero. This load is then distributed over the period of a year according to normalized monthly and daily usage profiles presented in Figure 6a creating a time series of predicted thermal demand. On the other hand, as was mentioned previously, the electric load is separated into two classes: fixed and flexible. The flexible load, available for management trough DSM, is chosen to encompass the largest energy consumers. The nominal weekly appliance usage plan with appropriate timeframes for shifting is depicted in Table 1. A notable addition to the plan was made with the inclusion of an electric vehicle (EV) with its effect on DSM signified by [28]. As this load accounts for a total of 17,672 kWh in yearly electric energy consumption, or 3694 kWh when not considering the EV, the remaining portion, attributed to appliances not listed as flexible (i.e., lighting, refrigeration, etc.), was chosen to equal 550 kWh and is distributed similarly as the thermal load using normalized usage profiles for electricity depicted in Figure 6b. In this case, four types of daily profiles are considered: winter work day, winter weekend day, summer work day, and summer weekend day and so repeating these segments appropriately results in a fixed electric demand time series as presented in the Figure 7. A thing to note is that both graphs from Figure 7 include daily and hourly noise in a similar fashion as is contained in popular microgrid software solutions like HOMER. The constructed baseline, having in mind tariff information given in Section 4.2, equates to a nominal spending of 952.9 EUR/a for gas and 3577.1 EUR/a for electricity totaling 4530 EUR/a if considering a general use case where the thermal demand is met with gas. This number, inclusive of standing charges, will be used later on when discussing feasibility of selected configurations.

Renewable Technologies
To perform the operation optimization, the models from Appendix A need to be instantiated with real-word values corresponding to the selected use case site. First of all, data from TMY records [29] for the given location is obtained to be fed into the RES models for estimated production and demand calculations. Starting with the WT, a power curve is assumed after the manufacturer data for the Lely Aircon 10 turbine. The obtained power curve, normalized for one unit of power of installed capacity, is shown in Figure 8a. When the specified model is used to estimate the power production of the same normalized WT over the course of a year, including parameters given at the top of Table 2, the model gives the results depicted in Figure 8a. Moving on, using the PV parameters from Table 2 for the Suntech STP250S-20/Wd monocrystalline solar module, we may estimate the yearly production of the PV array per unit power considering the appropriate weather data for the use case location, as is presented in Figure 8b.
Finally, the values required for STC and GSHP modeling are given in Table 2. When applying aforementioned TMY data in conjunction with this data to the appropriate models, having in mind the thermal demand profiles set in Section 4.3, Figure 7b is obtained showing the part of fixed electric demand required for running the GSHP, and thus, meeting the thermal demand of the system.
Given the fact that out of the flexible appliances that were modeled, only the EV is set to be dispersible with maximum power draw deviations of ±50 per cent of its nominal value, and also having in mind that splitting its charging process into several sections probably does not impact the user's comfort in a significant manner, ζ i (k) = 0 was set for all appliances i and all time steps k. As for the CO 2 emissions criteria, the carbon footprint values for energy sources were set to f C WT = 20 g/kWh, f C PV = 40 g/kWh and f C grid = 310 g/kWh in accordance with the fuel mix data provided by ScottishPower [30] and median values that can be obtained from renewable production statistics.

Available Configurations
Since the multicriteria analysis is conceptualized as an exhaustive search over a preselected domain, that domain must be specified before the optimization is performed. Therefore, a set of renewable generators and storage options is formulated as depicted in Table 3. The pricing for the WT is loosely based on extrapolations of market data from EnergySavingTrust [31], the PV pricing can be found in a detailed statistical sheet compiled by BEIS [32] and BESS parameters are closely related to those that can be found in the technical documentation for the LG CHEM RESU [33] battery series and quotes by renewable technology resellers. The rest of the system is not considered within the sizing optimization process. This includes the installation for the GSHP system (construction work, piping, fittings, etc.), the pump and the STC. The pump and STC are considered to have an expected life span of 20 years while the GSHP installation generally lasts around 100 years before needing to be replaced. The initial acquisition costs are considered to be 4500 EUR for the STC, 9500 EUR for the GSHP installation and groundwork and 7300 EUR for the pump given that a 13 kW system is to be installed which is estimated to fulfill the thermal demand of the considered household. The associated yearly maintenance cost with running all of the aforementioned renewable technologies is set to 2% of initial acquisition cost.

Results
The operational optimization was performed using the CPLEX® Matlab Toolbox for all 100 selected configurations in both DSM off and on states with a time constraint of five minutes per configuration and was finished in about 10 hours. An illustrative example of an optimized operation of a modeled appliance is given in Figure 9a with an overview of all relevant energy assets for the same time frame presented in Figure 9b. The results in the (C C , C NZEB , C CO 2 ) space of the proposed criteria for all considered configurations are presented in Figure 10. The Pareto frontier formed by the set of nondominated solutions depicting the best compromise-free solution set is accented over the rest of the obtained results.
However, given that the end result of the optimization process is considered to be only one, best solution, in the context of MCDMA ranking, the final outcome is highly dependent on the selected weights of each of the criteria that are imposed. Since each user will have different preferences, the decision space is practically infinite, however a couple of concrete examples are discussed to showcase the capabilities of the proposed planning strategy. Four different use cases will be presented to provide illustrative examples of how the optimum configuration changes in line with the criteria weight selection process. The first two will place a larger focus on monetary costs with the first one only considering them while the second introduces concerns for the supply and demand balance as well as for estimated environmental pollution. On the other hand, the third and fourth use case place a larger focus on the energy balance and emissions, respectively, while also including, but reduced, concern for monetary costs of running the system.

Total Cost as the Only Criterion
For the sake of simplicity, the first considered setup only focuses on minimizing the estimated total cost criterion C C through the appropriate set of weights w C = 100%, w NZEB = 0 and w CO 2 = 0. (58) In the optimization process, the results were first obtained without employing DSM giving a midway baseline, and the five most cost-effective solutions are shown in Table 4. As can be clearly seen, implementing only raw renewable sources without any additional optimization, in this case, does not yield significant monetary savings when compared with the gas and electricity no-renewable baseline, with the best option well below 3%. The third best, and all following combinations without DSM are not even profitable, with the fifth one costing over 3% more than the original baseline.
However, when DSM is turned on, the order of best solutions, depicted in Table 5, changes somewhat, with the best case scenario (using a 5 kW WT) returning savings slightly above 10% allowing the investment to be considered well worthwhile. A thing to note is that no PV production is included in this best-case scenario which should not be surprising since the considered site does not receive a significant enough amount of sunlight for PV to be cost-effective. However, the absence of BESS in the optimal solution is of significantly greater interest since the obtained result can be considered to show that an efficient implementation of a DSM program can lower the requirement for expensive and potentially environmentally unfriendly solutions like lithium-ion batteries that are commonly used for storing off-peak electric energy. Slightly less cost-effective then the best solution is the combination with the 7.5 kW WT, with the fifth best solution finally including the smallest considered PV panel array, and the sixth including a 3 kWh battery with the originally selected 5 kW WT.

Total Cost as the Primary Criterion
On the other hand, cost need not be considered as the sole criterion when selecting the optimal configuration. Choosing different weights allows for the user to specify what criteria he deems relevant, and thus, the software adjusts the optimal configuration rankings. One such mixed case where cost is still the primary focus, but the other two environmental criteria are also taken into consideration can be defined by selecting appropriate weights to equal w C = 60%, w NZEB = 20% and w CO 2 = 20%.
After ranking the optimally preforming configurations in accordance with the abovementioned weights, the best combination and a list of best alternatives is obtained and presented in Table 6. The results show that although monetary savings are still the prevalent considered parameter, this type of selection slightly favors large renewable sources due to the associated decrease in discrepancy between the spent and produced amount of energy and the equivalent amount of CO 2 emitted.

NZEB Rating as the Primary Criterion
Also, environmental criteria can be accented if the project is such that the ecological impact is more of a priority. One possible weight selection that would accomplish this could be w C = 20%, w NZEB = 60% and w CO 2 = 20%.
Adequate ranking resulted in the list presented in Table 7. Since the monetary aspect is no longer the primary concern, the best selected solutions are generally not profitable when compared to baseline costs. However, they are significantly more ecologically friendly, with all of them generating significantly more energy on-site than is being consumed and CO 2 being cut by about 40 to 50% when compared with that of the configuration without WT, PV, and BESS for all the presented alternatives.

CO 2 Emissions as the Primary Criterion
Finally, amongst environmentally friendly use cases, the main stressed criterion could be emissions, as is accomplished by w C = 30%, w NZEB = 10% and w CO 2 = 60%.
The resulting ranking in this case is presented in Table 8. As was expected, the best solutions are, just like in the aforementioned case, not the most profitable, but offer significant environmental features. Nonzero BESS units now appear in all of the selected configurations because of the need to minimize the emissions trough using locally produced electricity as opposed to grid imports.

Use Case without Incentives
In the end, for reference purposes, a use case is assumed where the various incentives that were previously considered are set to zero. This use case is obtained by specifying α WT = 0, α PV = 0, α grid = 19.6 cEUR/kWh, β WT = 19.63 cEUR/kWh, β PV = 19.63 cEUR/kWh, β grid = 0.
and also setting the RHI incentive to zero for both STC and GSHP. Since the process of energy trading on the wholesale market would certainly include multiple middlemen that would reduce the benefit that the user exporting energy can make use of, this use case assumes arguably generous renewable energy export prices at the same value as is the price of importing electricity from the grid. However, when the operational optimization is performed for configurations that include batteries with capacities of 0 and 9 kWh, PV arrays with capacities of 0 and 4 kW and all wind turbine configurations with nonzero capacities and the resulting results are analyzed with a sole focus on minimizing costs, the results from Table 9 are obtained. Results show that the event is the best-case scenario without incentives is more than 35% more expensive than the nonrenewable baseline, thus reaffirming the conclusion that the renewable installation project in the discussed case continues to highly depend on the existences of some of the mentioned incentive programs to be cost-effective.

Conclusions
The methodology presented in this paper details a comprehensive optimization process with a proposed linear programming-based model for operational optimization and a multicriteria ranking system for sizing optimization. To illustrate the capabilities of that process, a use case is assumed and the optimizer was employed to prove that the investment in renewable sources can be cost-effective when adequate appliance management techniques are used, resulting in savings just over 10% when only costs are considered as a primary criterion, in line with findings that could be found in the discussed related literature. The algorithm also manages to yield configurations with significant environmental impact improvements when compared to a nonrenewable baseline when the focus is shifted from costs to ecological factors. Also, in some of the selected cases, the obtained optimal solutions show an absence of storage solutions, which shows that DSM methods can be employed to reduce the negative environmental impact of costly chemical-based batteries that are often associated with HRES systems. Although the selected demonstration site was not comprehensive in terms of utilization of all energy sources featured by the Energy Hub methodology, additional sources can effortlessly be integrated into a similar analysis should the use case differ from the presented one.
A crucial aspect that should be considered for the demonstrated energy efficiency improvements to be realized in an arbitrary real-world setting are potential policy and legislative limitations. With different governments offering different terms for energy prosumers, local regulations should be, as was the done in the demonstrated use case, carefully analyzed when approaching both the planning and operation stages of the project. However, results presented in this paper and related publications can hopefully serve as a positive use case and aid in policy adjustments for the benefits of energy end users, energy producers, and grid operators.
Further research in this field is planned to include the implementation of more complex stochastic models into the production and load profiles as well as nominal appliance usage to more accurately model day-to-day variances in user demand. Furthermore, the demandside management aspect can be extended by also including demand response events that would facilitate load increases or decreases during predefined periods of time. Also, the sizing optimization paradigm could be extended to all sustainable sources like the heat pumps and thermal storage units. However, this extension may involve more detailed numerical models that would dynamically assess operation of the mentioned components to ensure long-term efficiency. Finally, the effects of different time step lengths and the introduction of different new criteria should also be analyzed to provide a holistic view of solutions that can be obtained in the desired decision space.

Conflicts of Interest:
The authors declare no conflict of interest.
v(z hub ) = v(z anem ) ln(z hub /z 0 ) ln(z anem /z 0 ) This formula is obtained as a nonlinear fit with experimental data where A is the amount of times that the area of the collector is greater than 1.33 m 2 used in the original experiment and parameter values a 1 = 1.637 m 2 , a 2 = 4.492 kJ/K and a 3 = 0.192 MJ are taken from the original experiment. If the produced value is negative, the collector is detached as not to further reduce the amount of energy stored in the tank. Otherwise, the additional energy is injected into the tank, and so the available amount of heat can be expressed as Q avail (k) = Q(k) − Q loss (k) + max Q prod (k), 0 (A13) However, since the user's demand Q demand at the considered time step is be subtracted, but Q avail may not be sufficient to meet all of the thermal needs. Thus, it is stated that the amount stored in the tank at the start of the next time step is where Q min (k) = mc(T amb (k) − T base ).
If the stored heat is limited by (A14) to Q min , the user's demand was not fully met, and the unmet heat is passed on to an auxiliary heating system whose role is played by the GSHP in the considered configuration.
Appropriate values for parameters k and λ used to fit these expressions are summarized in Table A1. Finally, having in mind that the amount of thermal energy that can be obtained by this system is limited by the designed power, and thus, may only be less than P max T s , we set the consumed amount of electric energy to where in the latter case auxiliary thermal sources must be used to fulfill the remaining part of the demand. k-th criterion out of q for a k d k (a i , a j ) pair-wise comparison value for criterion c k π k ((a i , a j )) preference degree function P k (a i , a j ) preference function p k , q k lower and upper preference boundaries w k weight associated with criterion c k φ + , φ − , φ positive, negative and net preference flow WT wind turbine P test WT (v v ) nominal power-wind speed WT power curve Y WT rated power (power capacity) of the WT system P WT (v v ) capacity and air density adjusted power-wind speed WT power curve ρ, ρ 0 air density at site and at test conditions B laps rate T 0 standard temperature g gravitational acceleration R universal gas constant divided by the molar mass of air z abs hub WT hub height above sea level z anem/hub anemometer and WT hub height above ground v anem/hub wind speed at anemometer and WT hub height z 0 surface roughness level PV photovoltaic P PV instantaneous power generated by the PV system Y PV rated power (power capacity) of the PV system f PV derating factor of the PV system G T solar irradiance incident on the PV array averaged at the current timestamp a P temperature coefficient of power T c PV cell temperature STC (at) standard test conditions G d/b , G diffuse, beam and global horizontal irradiation A i anisotropy index R b ratio between beam radiation on a tilted and horizontal surface f horizontal brightening factor β PV surface slope ρ g ground reflectance (albedo) T a ambient temperature NOCT (at) nominal operating cell temperature η mp maximum power point efficiency A PV surface area of the PV array a solar absorbance of the PV array cover τ solar transmittance of the PV array cover STC solat thermal collector Q tank/loss available amount of heat in the tank and lost amount of heat T tank/amb water temperature in the tank and ambient temperature m mass of water in the tank c specific heat capacity of water η l storage heat loss coefficient per unit of volume T s sample rate Q prod/avail produced amount of heat and available amount of heat in the tank GSHP ground source heat pump T g ground temperature T g mean annual soil temperature A s annual surface temperature amplitude X s soil depth α soil thermal diffusivity k soil thermal conductivity ρ soil density c p soil specific heat t day of the year T ewt entering water temperature T d cooling/heating averaged design temperature T bin binned ambient temperature χ capacity multiplier q c/d cooling/heating demand η, η base instantaneous coefficient of performance (COP) and at baseline P max maximum design power E elec,GSHP electric energy consumed by GSHP Appendix B.2. Mathematical model P in input (imported) power P cin/cout input and output power to and from the converters Q in/out power flow to/from the storage at the input and output stage L loads P out output power to loads P exp exported power q in/out charge/discharge rate for storage at the input and output stages E in/out available energy (state of charge, SOC) for storage at the input and output stages y, z device on/off and device start indicators d +/− positive and negative power deviations I(d +/− ) indicators for positive and negative power deviations S in/qin input stage energy storage conversion matrices F in/qin input and output transformation matrices C energy conversion matrix D exp , R exp export energy restriction matrices S out/qout output stage energy storage conversion matrices P i , P nom i i-th appliance's instantaneous and nominal power draw w nom i , ∆t nom i i-th appliance's n-th activation window indicator and length P max dev+/− i maximum absolute positive and negative power deviations P renew power produced from renewable sources SOC min/max in/out SOC capaciti limits at input and output stages J c/d/cd cost, dispersion and mixed criterion α, β, σ energy import and export prices and standing charge ζ i i-th appliance dispersion penalization factor X EMI/maint i equated monthly installments and maintenance costs γ i estimated lifetime in years for i-th device δ monthly discount rate B i initial acquisition cost for i-th device C C/NZEB/CO 2 cost, net-zero energy building (NZEB) and CO 2 emission criterion f C WT/PV/grid carbon footprint values for WT, PV and grid imported energy