A New Bi-Objective Approach for Optimal Sizing of Electrical and Thermal Devices in Zero Energy Buildings Considering Environmental Impacts

This paper proposes a new bi-objective optimization model, trading-off cost and environmental impacts, for sizing the key electrical and thermal devices in a zero energy building (ZEB), i.e., a building that roughly generates as much renewable energy as it consumes annually. A salient novel feature is the consideration of the environmental impacts, computed through a rigorous life cycle assessment approach, of buying electricity from the grid and manufacturing devices. Furthermore, an enhancement of the proposed model, as compared to the existing models, is to prioritize storing the ZEB excess of energy rather than selling it to the grid. The proposed solution approach of the initial mixed-integer nonlinear programming model relies on McCormick relaxation linearization to obtain a more tractable mixed-integer linear model. An augmented ϵ-constraint method is applied to solve the obtained bi-objective model. Finally, considering the building owners’ willingness-to-pay for environmental impacts, a decision-making criterion is proposed to select the optimal size of the devices among all non-dominated solutions of the Pareto front.

I e+ (t)/ I e− (t) Binary variables indicating that PV production is higher/less than electricity demand at time interval t I th+ (t)/ I th− (t) Binary variables indicating that solar thermal production is higher/less than heat demand at time interval t P e+ (t)/ P e− (t) Excess/Deficit of electricity at time interval t [kW] P th+ (t)/ P th− (t) Excess/Deficit of heat at time interval t [kW] P b (t)/ P s (t) Active power buy from/sell to grid at time interval t [kW] P c (t)/ P d (t) Charging T HE consequence of using fossil fuels to produce electricity is the increase in greenhouse gases emission. Buildings, as one of the major contributors to CO 2 emissions, are in charge of 40% of the final energy consumptions and 36% of CO 2 emissions in the European Union. To mitigate this, the idea of zero energy buildings (ZEBs) was imposed in the recast of the Energy Performance of Building's Directive in 2010 [1]. By definition, a ZEB is a building that roughly generates as much energy from renewable sources as it consumes annually [2].
Multiple technologies can be employed to fulfill the criterion of a ZEB for both renewable energy production side, e.g. rooftop photovoltaic (PV) systems or solar thermal collectors, and improved energy efficiency side (e.g., heat pump). Albeit the annual net energy consumption of a ZEB is almost zero, it still exchanges a substantial amount of electricity with the grid due to the mismatch between generation and consumption patterns. Accordingly, storage systems, e.g. building-integrated battery energy storage (BES) or heat storage tank, can be economically beneficial for the ZEB owner as well as for the grid (e.g. in congestion relieving manner) [3]- [6]. This paper focuses on finding the optimum size of each energy production, storage and energy-efficient device at ZEB planning stage that satisfies the zero energy constraint.
A BES sizing problem for residential ZEBs, in which real data of a household in Portugal is used, concludes that the designed BES system can reduce the energy export/import to/from the grid by 76% and 78.3%, respectively [7]. The compensation of the mismatch between the production and consumption patterns of ZEBs is addressed in [8]. A BES sizing problem for residential ZEBs is proposed in [9]. The model prioritized to store the excess electric power in the battery rather than selling it to the grid. This is an appropriate tactic since the electricity price is usually fixed and it is higher than the feed-in-tariff for residential buildings. However, the model suffers from not having a tractable mathematical formulation with the guaranteed optimality gap.
To meet the ZEB criterion, it is also important to consider the size of the devices on the thermal side of the building [10]. A cost-optimal ZEB design is studied for a typical house in Germany by [11]. Multiple possible technologies from both thermal and electrical standpoints are studied. It is shown that depending on the solar irradiation, energy costs, and building consumptions, some of the technologies might or might not be selected for a cost-optimal ZEB designing model.
All aforementioned technologies have certain impacts on the environment, mostly during their manufacturing. Furthermore, the electricity injected into the grid by generating units is also responsible for CO 2 emissions, which depend on the mix of technologies used to produce electricity [12]. The growing environmental concerns and green conscience of many people make it pertinent to consider further environmental aspects in the ZEB planning problem rather than only satisfying the zero energy constraint. This underpins the proposal of this paper, which is specifically in line with approaches from other areas [13], to consider a bi-objective optimization model trading-off an economic goal and an environmental goal.
The main contributions of the paper are as follows: r A new mathematical model and solution approach for the joint optimal electric and thermal device sizing of a ZEB is proposed. The proposed model is formulated as a mixed-integer nonlinear programming (MINLP) problem. Then, it is reformulated to reduce nonlinearities into only two bilinear constraints. Finally, the bilinear constraint is linearized using McCormick relaxation to obtain a mixedinteger linear programming (MILP) model. Further, the proposed model is enhanced to prioritize storing ZEB excess energy rather than selling it to the grid. This is a more realistic approach since the electricity price is usually fixed and it is higher than the feed-in-tariff for residential buildings.
r The proposed model is bi-objective, trading-off an economic goal and an environmental objective related to the impacts of buying electricity from the grid and manufacturing thermal and electrical devices. To solve the proposed bi-objective model, a tailored solution algorithm based on the augmented ε-constraint method of [14] is proposed. Finally, a decision-making criterion, which takes into account the building owners' willingness-to-pay for environmental impacts, is proposed to select the optimal size of the devices from all non-dominated solutions of the Pareto front. The rest of the paper is organized as follows. The ZEB configuration is described in Section II. The proposed optimal device sizing model and the solution algorithm are presented in Section III. Numerical results are presented in Section IV and conclusions are given in Section V. On the electric side, the building is equipped with PV and BES systems to supply the electrical demand. On the thermal side, a solar thermal collector along with a heat storage tank supply the thermal demand. Note that the heat pump plays an interfacing role between electrical and thermal sides. It consumes electricity and produces heat corresponding to its coefficient of performance (COP).
First, the optimal device sizing model is formulated as an MINLP problem. Then, the MINLP model is effectively reformulated such that the nonlinear terms are reduced. Finally, the bilinear constraint is linearized using McCormick relaxation to obtain a MILP model.
It is worth mentioning that the proposed model is capable of prioritizing supplier of the demand to maximize the efficiency and minimize the curtailment of renewable sources. That is, the model considers that the required electricity demand at time interval t is firstly provided by the PV and BES systems. The electricity demand exceeding the PV production plus the stored electricity in the BES system has to be bought from the grid. On the other hand, the excess electricity produced by the PV is firstly stored in the BES system. If the BES system does not have enough capacity to store the excess electricity, it is sold to the grid. A similar trend is considered on the thermal side.

A. MINLP Model
To model the economical objective function, first we should take into account the financial source. Assuming that the initial investment costs of the devices are financed with an interest rate of i for a period of n years, the capital recovery factor CRF can be calculated by (1) [15].
The economical objective function (2) of the proposed optimal sizing model includes costs (initial investment, operation and maintenance) of PV, BES, solar thermal, heat storage, and heat pump, as well as the cost of buying/selling electricity from/to the grid. Note that the proposed objective function (2) is also valid in the case of time-varying electricity price and feed-in-tariff (i.e., EP(t) and FIT(t) can be different for diverse intervals). The second objective function, which minimizes the environmental impacts, is modeled by (3). It includes the total environmental impacts of buying electricity from the grid, which depends on the prime sources of electricity generation, and the environmental impacts of manufacturing the devices.
The management of excess and deficit of electricity in the ZEB (see Table I) is modeled by a particular set of constraints as follows. Two binary decision variables I e+ (t) and I e− (t) are defined to respectively represent the case of having excess and deficit of electricity at time interval t. If the PV production at time interval t is higher than the total electricity demand at that time interval, then I e+ (t) = 1. Similarly, if the PV production at time interval t is less than the total electricity demand at that time interval, then I e− (t) = 1. Constraint (4) imposes that at each time interval one of the binary variables I e+ (t) and I e− (t) is equal to one, and the other one must be zero. In the case of having excess of electricity at time interval t, the excess of electric power can be calculated by (5). Similarly, for intervals that deficit of electricity exists, the deficit of electric power can be obtained by (6). As imposed by (7) and (8), the continuous variables P e+ (t) and P e− (t) are always non-negative.
Furthermore, two binary decision variables I e max and I e min are introduced to model possible statuses of the BES system at time interval t. If the battery has enough capacity to store the excess of electric power without reaching to its maximum state of charge (SOC) at time interval t, then I e max is equal to zero, and it is equal to one otherwise. Similarly, the binary variable I e min is equal to zero if the battery has enough energy to provide the deficit of electric power without reaching its minimum SOC at time interval t, and it is equal to one otherwise.
Introducing four binary decision variables for each time interval t, 2 4 = 16 possibilities exist. However, constraint (4) prunes eight infeasible possibilities. Constraints (9) and (10) must be considered to prune the other four infeasible possibilities and set the feasibility region according to the desired statuses presented in Table I.
The amount of electric power that can be sold to the grid is modeled by (11), where the battery charging power P c (t) can be calculated by (12). Likewise, the amount of electric power that has to be imported from the grid is modeled by (13), where the battery discharging power P d (t) can be obtained from (14).
The battery SOC at time interval t with respect to the desired statuses mentioned in Table I is modeled by (15). According to (15), the battery SOC is fixed to SOC max and SOC min for status 1 and status 3, respectively. Furthermore, constraint (15) models status 2 of Table I in which the battery SOC is equal to the summation of the excess of electric power stored in the battery and SOC(t − 1). Similarly, to model status 4, the battery SOC is equal to SOC(t − 1) minus the deficit of electric power provided by the battery. The energy preservation constraint (16) imposes that the battery SOC in the last interval should be equal to the initial SOC of the battery. Constraint (17) enforces the limitation of the battery SOC.
The same logic is proposed to model the solar thermal collector and heat storage tank on the thermal side of the building. Two binary decision variables I th+ (t) and I th− (t), which are constrained by (18), are defined to respectively represent the case of having excess and deficit of heat at time interval t. The excess and deficit of heat at time interval t are modeled by (19) and (20). As represented by (21) and (22), the continuous variables P th+ (t) and P th− (t) are non-negative.
Constraints (23) and (24) prune four infeasible possibilities and set the feasibility region similar to the desired statuses mentioned in Table I (similar logic as of (9) and (10)). With a similar interpretation as of (15), the heat storage SOC is modeled by (25). The energy preservation constraint and the limitation of the heat storage SOC are modeled by (26) and (27).
The heat pump is modeled by (28)-(31). Constraint (28) imposes that in case of not having enough stored heat to supply the deficit of heat demand (I HS min (t) = 1), heat pump must be committed. The relation between the heat pump generation and its electricity consumption is modeled by (29). The linear relation between the heat pump heat production (Q HP ) and the heat pump electric power consumption (E HP ) can be defined by the COP of the heat pump [16]. Constraint (30) limits the heat pump generation to be equal or less than the size of the heat pump. Constraint (31) is considered to select the size of the heat pump in a way that it can robustly supply the annual peak of the heat demand profile. According to the zero energy constraint of (32), the annual electric power drawn from the grid must be less or equal than the annual electric power injected into the grid.

B. Reformulated MINLP Model
To reduce the nonlinearity of the proposed MINLP model, auxiliary variables ω B (t) and ω HS (t) are introduced as (33) and (34), respectively. They represent the amount of energy that is stored in the optimal size of the battery and heat storage tank at time interval t.
To linearize the nonlinear constraint of (15), new linear constraints (45)-(47) are proposed. In status 1 when I e+ = I e max = 1 and I e− = I e min = 0, constraints (46) and (47) are inactive while (45) imposes that SOC (t) = SOC max . This is the same conclusion obtained from (15). Checking other statuses of Table I, it can be concluded that linear constraints (45)-(47) are totally equivalent to (15).
Finally, the proposed reformulated MINLP model for optimal device sizing problem is the minimization of objective functions (2) and (3)  Clearly, the proposed reformulation reduces the number of nonlinear constraints. Now, bilinear terms of (33) and (34) are the only nonlinearity of the model. The fact that current versions of standard solvers, e.g. CPLEX [17] and Gurobi [18], can solve MILP models very efficiently motivates us to propose a MILP model for the optimal device sizing problem.

C. MILP Model Using McCormick Relaxation
The MILP model of the optimal device sizing problem is derived using a McCormick relaxation technique [19]. Mc-Cormick underestimators (57) and (58) along with McCormick overestimators (59) and (60) can be added to the model to relax the bilinear constraint (33). Similarly, the bilinear constraint of (34) can be replaced by McCormick envelopes (61)-(64). The advantage of McCormick relaxation over other approximation techniques is that the relaxed MILP model contains the feasibility region of the original MINLP model completely. Ideally, the McCormick envelops represent a convex hull for the MINLP feasibility region [19].

D. Solution Algorithm
The ε-constraint method is well-established and widely used to generate the Pareto front of bi-objective optimization problems. In this work, we use an advanced version of this method, namely the augmented ε-constraint approach [14], which allows to obtain efficiently an even distribution of non-dominated solutions (the Pareto set).
In the ε-constraint method, one of the objectives is optimized while the other one is considered as an inequality constraint. By changing the value of the parameter on the right-hand-side of the constrained objective, efficient solutions of the single objective problem can be obtained. Two main challenges of the ε-constraint method are finding the range of the constrained objective and the domination guarantee of the obtained solution. Finding a valid range for the right-hand-side parameter of the constrained objective is not a trivial task, especially the worst value over the set of non-dominated solutions (i.e., nadir value). Besides, the obtained solution is non-dominated only if the constrained objective function is binding. These challenges are alleviated in the augmented ε-constraint approach by incorporating slack variables and transforming the constrained objective into an equality constraint. More details of the augmented ε-constraint approach are presented in [14].
After obtaining the Pareto front for objective functions (2) and (3), the proposed decision-making criterion of (65) can be used. Parameter WTP is the willingness-to-pay of the building owner to reduce the environmental impacts. According to previous studies and surveys, this value differs around the world [20]. Criterion (65) represents the fact that the difference between the cost of the proposed bi-objective model considering environmental impacts and the single objective model ignoring the environmental impacts should be equal to (or less than) WTP. It will be shown numerically in the next section that the value of WTP can be between zero and a positive upper bound. For the values of WTP larger than the positive upper bound, the model would be infeasible.

A. System Description and Simulation Setup
To illustrate the advantages of the proposed model, a typical residential building in Luxembourg is studied. Solar irradiation data and per unit output power of roof-top solar panels in Luxembourg are obtained from [21]. For the annual electricity consumption of the residential building, CREOS (distribution electricity utility in Luxembourg) database is used [22]. Furthermore, the annual heat demand of the building is acquired from [10]. The annual electrical energy consumption and heat demand are 4738 kWh and 4147 kWh, respectively. The environmental impact of each device depends on its energy technology and manufacturing process. Life cycle inventories for all technologies have been extracted from the ecoinvent 3.6 database [23], and adapted to Luxembourgish conditions. In particular, the following inventories were used: "Market for electricity, low voltage, (Luxembourg)", "Photovoltaic slantedroof installation, multi-Si (Global)", and "Battery cell, Li-ion (Global)". The heat storage system consists of a basic chrome steel tank, which is equipped with a heat exchanger and a boiler, made in Europe. The solar thermal system consists of a solar flat-plate collector for a one-family house with a combined system for hot water and heating. The heat pump is an air-water heat exchanger operated in Luxembourg. The environmental scores are calculated in a sequence of distinct steps. First, the lifecycle inventories, which contain thousands of environmental flows (emissions to air, water, soil, or raw material extraction), are characterized into a series of so-called "midpoint" indicators, which convey the environmental stress generated indirectly and directly by each system via various impact mechanisms. For example, the "climate change" impact is accounted in kg CO 2 equivalents (a reference unit to convert the set of greenhouse gas emissions into one single indicator) using their global warming potentials over a 100-year horizon as the characterization factors. The characterization method chosen is the one recommended by the European Union's Joint Research Centre in the context of the Environmental Footprint methodology [24]. Each midpoint indicator is expressed in "points" (Pt), where one point represents the global impact of that category divided by the world population (yielding the global per-capita average for that category). Once expressed in the same unit, results can be aggregated into a single score, as shown in Table II.
The following assumptions are considered: r According to [9], [11], the initial investment costs of PV, BES, solar thermal collector, heat storage tank, and heat pump are assumed to be 1500 €/kWp, 50 €/kWh, 800 €/kWp, 40 €/kWh, and 2500 €/kWh, respectively. The annual operation and maintenance cost of PV and solar thermal collector is assumed to be 1% of their investment costs, while it is assumed to be 2% for other devices [11]. The interest rate of the financial institution is assumed to be i = 5%.
r The warranty period (or the life cycle) of the devices is assumed to be greater than the payback period of the loan. Therefore, the building owner will not pay any replacement cost over the life span of the systems. This is a reasonable assumption since, for example, the life cycle of a typical Lithium-Ion battery can be more than 15 years [25], [26].
r Capacity to power ratio (CPR) of the battery is assumed to be small enough to fully charge or discharge the battery between two consecutive intervals. It is a reasonable assumption since the BES systems for the application of the proposed device sizing problem are not high-storage (i.e., in the range of several kWh) and current technologies of charging systems are able to provide continuous power up to few kWh.
r The electricity price and feed-in-tariff in Luxembourg are 0.17 €/kWh and 0.121 €/kWh, respectively. To solve MILP and MINLP models, ILOG CPLEX 12.9 and Baron 19 are respectively used. All solver options are set to defaults except CPLEX parallel computing options "cpx_param_parallelmode" and "cpx_param_threads", which are respectively set to 1 and 48. We used a full node of a high-performance computer (up to 187-GB-RAM) with 48 Intel Skylake @ 2.1GHz cores per node.

B. Numerical Results
Applying a K-means++ clustering algorithm, the size of the hourly raw dataset is reduced to 365 and 12 to respectively represent average daily and monthly profiles. The daily and monthly clusters of per unit PV output, electricity consumption, and heat demand are depicted in Figs. 2, 3, and 4, respectively. K-means is a model-free method for data clustering and partitioning. Given the number of clusters (K value), it aims to partition the raw dataset into K clusters in which each data  belongs to the cluster with the nearest mean (cluster centroid). After assigning all data to clusters, recalculated centroids are considered for the next iteration. The algorithm converges when the assignments no longer change. To reduce the dependency of the K-means algorithm to the initial centroids, we apply an initialization approach, called K-means++, that is proven to work better than K-means. In addition, the occurrence time of the data is kept unchanged during the clustering process. Therefore, time-dependent behavior is preserved in the final clusters. Details about these algorithms can be found in [27].
As expected, it can be observed that in colder seasons the electricity and heat consumptions are higher, while the solar generation is lower than in warmer seasons.
1) Ignoring Environmental Impacts: Considering first only the economical objective (2), the simulation results of the proposed  Table III. In order to guarantee the optimality of the solution, the tolerance on the optimality gap is set to zero.
As shown in Table III, investing in the solar thermal collector and heat storage tank is not economical for the case of daily and hourly intervals. One explanation could be the fact that the size of the heat pump is usually selected by the robust constraint (31). Thus, its value must be larger than the peak of the heat demand profile. For the data with higher resolution (e.g., 365 and 8760 clusters), the spiky nature of the profile leads to a higher peak. Therefore, the model suggests a larger heat pump for the case of high-resolution data and a smaller heat pump for the case of monthly data resolution. Since none of the thermal devices can help in providing electricity for the building (see Fig. 1), PV is the only production device responsible for satisfying net zero electricity building constraint. Therefore, similar to the heat pump, the presence of PV is necessary. The model suggests that, due to high COP of the heat pump, it is more economical to supply the heat demand using the electricity provided by the electrical side (PV+BES+grid) through the heat pump (which must exist due to (31)) rather than investing in a solar thermal collector and a heat storage tank. However, it was not the case for monthly intervals due to the smaller size of the heat pump. It should be noted that these results are for a specific geographical location (Luxembourg) with the assumed cost of the devices. Assuming other investment costs for the devices or having more solar radiation might lead to a different decision by the proposed model. It can be observed that the solution time of the MILP model is considerably smaller than that of the MINLP one. For the case of 8760 intervals, the MINLP model does not provide any feasible solution within 100 000 seconds. It demonstrates the intractability of the MINLP model for high-resolution data.
Note that the McCormick-based MILP model provides similar solutions for the case of 12 and 365 intervals with the benchmark MINLP reformulation. This shows empirically that the McCormick relaxation is tight and provides the global optimal solution in this case. Furthermore, similar to all planning problems (e.g., generation and transmission expansion planning [28]), the optimal size of each device depends on the data resolution. Therefore, the preprocessing analysis of raw data is a key step for the device sizing model.
After obtaining the optimal sizes, the values should be rounded up to the closest standard size available. For instance, a PV panel with a size of 6.5 kWp can be selected according to the solutions of the cases with the daily and hourly resolution.
The SOC of the battery for the case of 365 intervals is depicted in Fig. 5. Comparing Fig. 5 with Figs. 2-4, one can conclude that in the warmer seasons, when the consumptions are lower, the battery stays fully charged more often. However, in the colder seasons, when the consumption is higher, the battery stays fully discharged regularly.
As the device sizing problem is in the planning stage, due to the computational burden, it is not possible to model the real-time operation precisely. Having a fast-response control system, in real-time operation, leads to so narrow intervals that it is unlikely that the battery is fully-charged or fully-discharged for a long time. Besides, if one wants to avoid that the battery is fully-discharged, one can add a new reserve constraint to the model or simply set the SOC min B to a larger value. Fig. 6(a) shows a 3D visualization of the bilinear term of (33). Figs. 6(b) and 6(c) visualize the McCormick overestimators and underestimators, respectively. To make the maximum deviation observable, a side view from the angle of an observer located on the intersecting line of the underestimator planes is depicted in Fig. 6(d). As visualized by 3D figures and illustrated by case studies, the McCormick relaxation is tight enough to practically results in a high-quality solution of the proposed device sizing problem.
2) Bi-Objective Model: The proposed bi-objective device sizing model is solved using the ε-constraint method presented in Section III. For the case of daily intervals, Pareto optimal fronts of the reformulated MINLP and McCormick MILP models are shown in Fig. 7. As it can be observed, it is possible to control the number of efficient solutions by properly selecting the number of grid points within the range of objective functions in the ε-constraint algorithm. While the non-dominated solutions of  the MINLP model are selected manually, a number of 50 efficient solutions are automatically generated for the MILP model using the augmented ε-constraint method, which is explained in detail in [14]. The similarity between MINLP and MILP Pareto fronts is another illustration of the tightness of the McCormick relaxation.
It can be seen that it is not possible to reduce the total cost of the environmental impacts less than a lower bound (i.e., 3.0512 Pt). In other words, for the values of WTP larger than a positive upper bound (i.e., 1,447.5 − 1,353.97 = 93.53 €/year), the model becomes infeasible because the environmental objective function (3) is conflicting with the zero energy constraint (32).
After selecting the value of WTP, one can refer to the Pareto front to find the solution with respect to environmental impacts. For instance, for an owner who is willing to pay a maximum amount of 70 € per year, the environmental impacts index can be reduced from 4.9381 Pt to 3.279 Pt (i.e., 33.6%). The optimal size of thermal and electrical devices of a ZEB for such an owner is shown in Table IV. As shown in the table, the value of WTP is spent on selecting a larger BES as the storage system and a larger PV as the production source. Several reasons exist for this selection. First, as shown in Table II, the environmental impacts of buying electricity from the grid are very high. Thus, the optimal size of the thermal and electrical devices should be selected in a way that the amount of buying electricity from the grid is reduced. While the cost of buying electricity from the grid in the case of considering only the economical objective is 366.78 €/year, the proposed bi-objective model reduces it now to  247.78 €/year. Second, for the studied building, investing in the solar thermal collector and heat storage tank is not economically desired, as illustrated in daily and hourly results shown in Table III. Finally, the environmental impacts of manufacturing a larger heat pump are higher than those of a larger PV or BES (refer to Table II). Therefore, a larger BES and a much larger PV, which is close to 10 kWp (the assumed maximum limit of PV capacity), is selected to reduce the environmental impact of the ZEB. The results for the maximum possible WTP ( = 93.53 €/year) are also shown in Table IV. In this case, the optimal size of the BES system is the assumed maximum limit (i.e., 5 kWh) in order to further reduce the amount of electricity bought from the grid.
3) Sensitivity to Demand and Solar Variation: Considering up to 10% and 20% variation respectively in electricity demand and solar radiation, a sensitivity analysis is performed. The sensitivity analysis results for the optimal size of the BES and PV in the case of daily data resolution and WTP = 0 are shown in Figs. 8 and 9. Note that the values of optimal size of BES system becomes 10 times larger for the sake of visibility.
As shown in Fig. 8, the optimal size of both BES and PV is almost the same when up to 10% electricity demand reduction happens. However, by increasing the electricity demand of the building, larger PV and BES system are required to supply it. Due to limitations on PV installation (10 kWp for the studied building), by increasing the electricity demand more than 8%, the optimal size of the PV is selected to its maximum possible value. In this case, depending on the cost of each option, either a larger battery must be installed or the amount of electricity bought from the grid must be increased.
By increasing the solar radiation from 80% to 112%, the optimal size of the PV is linearly decreasing. Having more solar radiation requires a smaller PV to satisfy demand. However, after some point (112% in this case study), the amount of solar radiation becomes high enough to not only meet the demand but also produce extra energy for the sake of selling it to the grid (energy arbitrage). Therefore, installing a larger PV becomes desirable after this point. In addition, by increasing the solar radiation from 80% to 104%, the optimal size of the BES system is fixed. However, by increasing the solar radiation to more than 104%, larger storage is selected.

V. CONCLUSION
This paper has proposed a new bi-objective mathematical model for the optimal sizing of key electrical and thermal devices to fulfill the requirements of zero energy buildings. To break down the computational complexity of the initial MINLP model, we proposed a more tractable MILP model, based on Mc-Cormick relaxation linearization. Numerical results on a typical residential building in Luxembourg illustrate that the proposed MILP model is computationally effective and, contrary to the MINLP model, scalable to large problems with high-resolution data. In our numerical experiments McCormick relaxation was tight, likely due to the rather mild model nonlinearity resulted from the proposed MINLP reformulation.
In addition to the cost-optimal objective function, the environmental impact, computed from life cycle assessment, of buying electricity from the grid and manufacturing thermal and electrical devices has been taken into account as a second objective function. An augmented ε-constraint algorithm was applied to obtain the Pareto front of the proposed bi-objective model and has shown good performances in producing evenly distributed solutions. The Pareto front enables building owners' to make informed environmentally friendly decisions in selecting the optimal size of the devices according to their willingness-to-pay for environmental impacts.
As illustrated by the sensitivity analysis, over-and underestimating the amount of building's energy demands and solar radiation affect the optimal size of the devices. Thus, we plan to consider uncertainties in the proposed model in future work.