Uncertainty quantification for combined building performance and cost-benefit analyses

Building and Environment 62 (2013) 143e154

Contents lists available at SciVerse ScienceDirect

Building and Environment journal homepage: www.elsevier.com/locate/buildenv

Uncertainty quantification for combined building performance and cost-benefit analyses Sebastian Burhenne a, *, Olga Tsvetkova a, b,1, Dirk Jacob a, Gregor P. Henze c, Andreas Wagner d a

Fraunhofer Institute for Solar Energy Systems, Freiburg, Germany Masdar Institute of Science and Technology, Abu Dhabi, United Arab Emirates University of Colorado, Boulder, USA d Karlsruhe Institute of Technology, Karlsruhe, Germany b c

a r t i c l e i n f o

a b s t r a c t

Article history: Received 11 August 2012 Received in revised form 26 November 2012 Accepted 12 January 2013

Building performance simulation is most often used to improve the design and at times the operation of buildings. Within a building model, the thermal characteristics of the envelope and the HVAC (heating, ventilation, and air conditioning) equipment are described by parameters that often cannot be estimated with high accuracy (e.g., occupant behavior, building envelope and HVAC equipment performance). Another common part of the design process of a building is a cost-benefit analysis to compare design options and different scenarios. The results are also heavily dependent on assumptions about uncertain economic parameters (e.g., future inflation rates and energy costs). In this paper a Monte Carlo based methodology for uncertainty quantification that combines the building simulation and the cost-benefit calculation is developed and demonstrated. Furthermore, Monte Carlo filtering is applied to determine the model inputs (e.g., design specifications and boundary conditions) that lead to the desired model output (e.g., a positive net present value of the investment). The aim is to propose a methodology that helps to enhance the design process or building operation and supports related decision-making.
2013 Elsevier Ltd. All rights reserved.

Keywords: Uncertainty quantification Uncertainty analysis Sensitivity analysis Building performance simulation Cost-benefit analysis Monte Carlo simulation

1. Introduction Uncertainty analysis (UA), or alternatively uncertainty quantification (UQ), when applied to building performance simulation (BPS) can improve the design process by considering the uncertainty in important design parameters that impact decision making. Several studies have revealed this very fact. Lomas and Eppel [1] conducted one of the first studies involving sensitivity analysis (SA) using Monte Carlo (MC) techniques for building simulations in 1992. They compared different SA methods and analyzed over 70 uncertain inputs. Macdonald et al. compared different techniques for uncertainty quantification and implemented some methods in a simulation program [2e4]. De Wit et al. [5,6] state, that it is essential to communicate uncertainties to decision makers. But in building simulation practice, an explicit assessment of uncertainty is more the exception than the rule. Therefore, most of the

* Corresponding author. E-mail address: [email protected] (S. Burhenne). 1 Olga Tsvetkova has been associated with Masdar Institute of Science and Technology, Abu Dhabi, United Arab Emirates since September 2011. 0360-1323/$ e see front matter
2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.buildenv.2013.01.013

decisions are based on point estimates. If a building simulation is considered as a decision support instrument, it seems to be necessary to assess and communicate the problem of uncertainties properly. Otherwise decisions might be made, which are based on faulty assumptions. This can lead to disappointments when the actual building does not perform as communicated through the design process. With a classical building simulation the answer to design questions is often either yes or no depending on the assumptions. An example for such a design question could be: Will the renewable energy system supply 20% of the total yearly energy demand? With a UA it is possible to answer design questions with probabilities, e.g.: The probability of more than 20% energy supplied by renewable energies is 80%. Furthermore, it is possible to compare design options and choose the option with the highest probability of reaching a specific goal. De Wit and Augenbroe [6] analyzed uncertainties in building design mainly related to thermal comfort. They identified uncertainties in ventilation rates and indoor air temperature distribution as highly influential on the simulation result. Furthermore, they proposed a method for decision making in building design. Mara and Tarantola [7] performed a sensitivity analysis for a test cell model and focused on SA to aid the model development. They

144

S. Burhenne et al. / Building and Environment 62 (2013) 143e154

calculated the sensitivity indices of simulation inputs using a variance-based method. Burhenne et al. [8e10] performed UA and SA for BPS applied to the design as well as to the operational stage of a building. Furthermore, different ways of generating random numbers were analyzed with respect to their performance (i.e., convergence and robustness). Eisenhower et al. [11] performed an UA and SA investigating the influence of about 1000 parameters using quasi random sampling and a meta-model (response surface). Hopfe and Hensen [12] distinguish between three types of uncertainties (physical, design and scenario uncertainties). They report results of a MC based UA and SA implemented with Latin hypercube sampling (LHS). In their work, the physical UA was conducted by varying physical properties of materials, while the design UA was accomplished by adjusting geometry and glass surfaces. The scenario uncertainty was implemented by varying the infiltration rate and internal loads. Booth et al. [13] performed an UA of housing stock models using Bayesian techniques. Their aim is to analyze different retrofit options at an urban scale. In practice many decisions are based on their monetary value. It seems to be a natural way to talk about different options in the design phase in terms of money. Particularly, for people without engineering background it is appealing to make a decision on the basis of a monetary value. Cost-benefit analysis (CBA) is a common approach to making economic decisions. The main idea underlying CBA is that total expected costs are compared with the total expected benefits of a project in order to choose the most profitable of the compared options. Benefits and costs are usually expressed in terms of money [14]. The duration of the project considered plays an important role. Usually the cash flow of costs and benefits occur at different times and the analyzed time frame for buildings is usually 5e50 years long. The longer the time frame, the greater the difference in monetary value becomes. Due to this fact, costs and benefits need to be adjusted for the time value of money. All cost and benefit flows in the course of the project (which usually take place at different points in time) are expressed in a common term, their present discounted value (PDV) [15]. Rysanek and Choudhary [16] conducted an analysis where they combined building simulation and economic calculations. The economic uncertainties were implemented with a best and a worst-case scenario and different refurbishment options were analyzed by changing the technical system in the building simulation. In this paper we combine UA, BPS and CBA. An uncertainty quantification for the uncertain BPS and CBA inputs is performed. Sampling based on Sobol’ sequences is used because quasi random sampling proved to be more effective in terms of convergence and robustness compared to pseudo-random sampling, stratified sampling and Latin hypercube sampling [10,11]. A method for quantifying the uncertain inputs of the CBA is outlined. This paper employs a methodology for performing the analysis that is scalable depending on individual project requirements and can be applied to answer several questions common in building design and operation. The methodology consists of the characterization of the uncertain input, employing an effective sampling method, i.e. sampling based on Sobol’ sequences, and recommendations concerning the coupling of the BPS and the CBA. Practical advice on how to implement the methodology is given and appropriate ways of displaying the results of the analysis are shown. MC filtering is applied to determine which input parameter drives the result into a certain region (e.g., efficient design, positive net present value). This application is of high value to supplement the typical building design process. It can lead to more efficient systems and an improved design. The ultimate goal is to support decision makers and implicitly take salient uncertainties into account in the decision making process. Furthermore, insights are

gained how technical systems can be improved in order to reach performance goals or requirements. The proposed methodology is applied to an example case. Main results are that the probability of more than 20% of the total energy demand supplied by renewable energies is 60%, whereas the probability that the investment in the solar thermal equipment is profitable is 9%. The most influential parameters, which lead to these results, are identified and it is analyzed under which conditions higher probabilities can be reached. This article is organized as follows: the methodology section covers the theoretical aspects of the work, i.e. modeling approach, CBA, ways of quantifying input uncertainties, sampling related issues and MC filtering. The simulation section describes the BPS model used in the case study, the CBA of the case study and the used tools. The results and discussion section presents the results obtained in the case study: BPS results, CBA results and MC filtering results. Finally, the conclusion and outlook section outlines the conclusions of the work carried out and suggests further directions for research. 2. Methodology 2.1. Modeling approach The literature introduced in the introduction shows that the most commonly applied approach for UA in BPS is based on MC techniques. However, it is also possible to reformulate models in a way that they can take probability density functions (PDF) as input and compute a PDF as output rather than producing point estimates of the simulation output. A study conducted by Jacob et al. [17] showed that this method requires extensive mathematical modeling efforts and more computational resources than MC techniques for a given accuracy and building problem. The analysis described in this paper is conducted using a 2-stepprocedure. The first part is the MC simulation of the performance of the building and its HVAC equipment. The output of the BPS (e.g., the amount of energy supplied by a solar thermal system) is used as input for the cost-benefit MC-analysis. Fig. 1 shows an overview of the simulation process. All CBA inputs, which are considered uncertain, are sampled except for one input of the CBA that is the result of the BPS. The histograms represent the inputs that are considered as uncertain and the uncertain output, respectively. Besides the uncertain inputs there are many parameters and variables that are considered as known. Hence, point estimates are used for these inputs. The decision if a parameter or variable is considered as known or as uncertain can be based on several criteria. These are the project scope and associated questions to be answered, the available information, the resources and time that can be spend to conduct the analysis and the specific method that is applied. In most cases it is desirable to consider the most uncertain and the potentially most influential parameters as uncertain. A good example for such a variable is the air change rate of a natural ventilated building since it is known to be influential on the energy demand and highly uncertain because it depends on the occupant’s behavior. Furthermore, the inputs that can be easily changed in the design or operation phase of a building are of great interest. The mass flow rate of the solar thermal collector is an example for such a variable. 2.2. Building performance simulation model One part of the proposed methodology is a dynamic BPS. The level of detail might vary depending on the analyzed case. However, the computational cost of the model plays an important role in the MC setting. The possibility for parallelization can help

S. Burhenne et al. / Building and Environment 62 (2013) 143e154

145

BPS

CBA point estimates

point estimates Fig. 1. Overview of the MC simulation with consecutively conducted BPS and CBA.

to deal with computational expensive models. Details on the model that is used in the case study are outlined in Chapter 3 Simulation. 2.3. Cost-benefit analysis In a CBA expected costs are compared with the expected benefits of a project. The different times at which cash flows take place require an analysis that takes interest and compound interest as well as inflation rates and energy price escalation into account. This is done by calculating the present discounted value (PDV) of every cost and benefit over the analyzed time span [15]. The PDV can be calculated according to

PDV ¼

n X

CFt

t ¼1

1 þ IR

!

(1)

1 þ Inflt where CFt is a future cash flow at time t, t is the year at which the cash flow takes place, IR is the nominal interest rate, Inflt is the average inflation rate from year 1 until year t and n is the number of years taken into account in the analysis. Inflt can be calculated according to

Inflt ¼

t Y

!1 t

1 þ Inflz

1

(2)

z¼1

where Inflz are the inflation rates from year 1 till year t. The final decision criterion used in this paper is the net present value (NPV) over the life cycle (n) of the analyzed technical system. The NPV is defined as the difference between the PDV of the benefits and the PDV of the costs [15] (see Equation (3)).

NPV ¼ PDVbenefits  PDVcosts

compared, the one with the highest NPV is the most beneficial option. The international standard EN 15459 “Energy performance of buildings e Economic evaluation procedure for energy systems in buildings” contains specific information on calculation methods, cost types (e.g., maintenance cost) and typical life spans for HVAC equipment [18]. However, it does not contain information concerning UA for economic evaluations. 2.4. Quantification of input uncertainty A correct quantification of the input uncertainties is essential to obtain reliable results by performing these analyses, i.e. requiring appropriate probability distributions and associated parameters. Depending on the project and the question to answer the knowledge about the uncertainty varies, and hence the approach to quantify the input uncertainties may vary as well. Saltelli and Tarantola point out that the quantification of the input uncertainty can be based on measurements, estimates, expert judgment, physical bounds, output from simulations and analogies to similar simulation input [19]. Several studies exist that deal with the quantification of input uncertainty in BPS (e.g., [3,5,20,21]). In the methodology section of this paper we focus on predicting the future trends for economic data and the corresponding uncertainty. In econometrics time series models such as ARIMA (auto regressive integrated moving average) models are often applied to predict future values based on time series data. In the following the basics of ARIMA models are briefly introduced. Univariate time series of economic variables can be used to fit an ARIMA model. An ARIMA model without differencing (it would be an ARMA model) can be described by

(3)

A positive NPV indicates that the analyzed option is beneficial for the given set of assumptions. When different design options are

yt ¼

3t

þ

p X i¼1

ai yti þ

q X j¼1

bj 3 tj

(4)

146

S. Burhenne et al. / Building and Environment 62 (2013) 143e154

where 3 is a noise term, a is the coefficient of the autoregressive (AR) part of the model, p is the order of the AR part, b is the coefficient of the moving average (MA) part of the model, q is the order of the MA part. Differencing is necessary in the case of nonstationary time series. The order of differencing is called d. Differencing is performed by replacing yt with Dyt ¼ yt  yt1 [22]. Furthermore, it is possible to take seasonal effects into account when fitting the ARIMA model. When seasonal effects are taken into account the model is sometimes called SARIMA model. In the case of a model with seasonal effects, the order of the different model components (AR, MA, I) of the seasonal part have to be set. These are P for the AR part of the seasonal model, Q for the MA part and D for the order of differencing. The interested reader is referred to Ref. [22] for further information concerning ARIMA models. In this paper ARIMA models were used to predict future values of the inflation rate and the gas price. The parameter identification of the coefficients for the models can be done with the help of statistical software. The orders of the different model parts have to be chosen by the analyst. An iterative process can do this (e.g., Box and Jenkins method [23]) or an automated parameter variation can be performed. Another possibility is to use an optimization algorithm for this task. After such a parameter variation the log likelihood or the Akaike information criterion (AIC) can be used to choose the best order of the different model parts. Both serve as relative measures for the model fit. A model setting with higher log likelihood is preferred over a model with lower log likelihood, which is based on the maximum likelihood method [22]. The AIC includes a penalty term for the number of parameters to avoid too many parameters in the models and thus weights model performance against model complexity. This is important because a model with a large number of parameters could produce a near perfect fit to the data but noise would be inherent in some of the model parameters. In this case, these parameters do not contain useful information and should not be used in the model. A model setting with lower AIC is preferred over a model setting with higher value for the AIC [22]. The AIC is calculated with the following equation:

AIC ¼ 2k  2lnðLÞ

(5)

where k is the number of parameters and L is the value of the likelihood function of the model. With the help of an ARIMA model, an instance of the future values of the economic data is predicted. The analyst should check if the prediction looks reasonable for the analyzed case. However, the prediction is subject to uncertainties and therefore the predicted time series serve as base case for the future trend of the economic input. The predicted time series are varied in the MC simulation by means of scaling factors to account for different scenarios. For the determination of the scaling factors the confidence interval or the standard error of the predicted time series are used.

strategies were analyzed with respect to the convergence of the estimate of the mean value and their robustness. Latin hypercube sampling and sampling based on Sobol’ sequences [25,26] had the fastest convergence of the mean estimates and outperformed random sampling and stratified sampling. Moreover, the study showed that the sampling based on Sobol’ sequences produced the most robust results. For the referenced study the same BPS model as in this paper was used. A study by Eisenhower et al. [11] also showed, that the convergence rate of quasi-random sampling techniques (sampling based on Sobol’ sequences is a quasi-random sampling technique) is faster than the one of pseudo-random sampling algorithms. Quasi-random numbers can be applied similar to pseudorandom numbers. The difference is that successive sampled points take the position of previously sampled points into account [27]. Hence, quasi-random sampling techniques generate samples that are distributed as uniformly as possible over the multi-dimensional parameter space. In this paper sampling based on Sobol’ sequences is used. As a measure of convergence, the mean and variance estimates of the model output are calculated and compared for an increasing sample size. 2.6. Monte Carlo filtering In the experiment setting of this paper it is interesting to analyze which model input values drive the model output into specific regions. An example of such a range is the solar fraction above a certain threshold or the range of positive NPV in the case of the CBA output (more details on that will be explained in the case study of this paper). MC filtering isolates the input of the model realizations that produce the desired result. With this method it is possible to infer the necessary design and settings of the technical systems that lead to the desired outcome (e.g., a positive NPV). Hence, the method is used to determine the subset of the model input space that drives the model output into specific regions (Fig. 2). The main advantage that MC filtering provides is insight under which conditions a certain design target is reached. This might lead to design improvements as well as more insights in the building and HVAC characteristics. In Fig. 2 (XjB) is compact, which might not be the case depending on the analyzed simulation. Performing MC filtering consists of the following steps [28]: 1. Separating the MC simulation output and the corresponding simulation input into two subsets. One is the behavioral subset (e.g., where the desired design goal is reached; B) and the remaining is the non behavioral (e.g., where the desired design goal is not reached).

(X | B )

2.5. Sampling Sampling is the process of generating random numbers. A random number is a mathematical definition and “real” random numbers cannot be generated in a computer experiment because computers use algorithms for sampling. Therefore random number generators are often and more correctly called pseudo-random number generators [24]. The sample size, which is required until a MC simulation converges, depends on the number of input parameters (k), the analyzed statistical measures (e.g., mean and variance estimates) and the properties of the analyzed model (e.g., nonlinearity). In a paper by Burhenne et al. [10], different sampling

(X | B )

B B

Simulation input

Simulation output

Fig. 2. Mapping of input and output subsets. B is the desired region of the simulation output (e.g., positive net present value) and (XjB) is the subset of the simulation inputs that results in B (adopted figure from Ref. [28]).

S. Burhenne et al. / Building and Environment 62 (2013) 143e154

2. Plotting the empirical cumulative distribution functions (ECDF) for each input parameter for the full set of simulations, the behavioral subset and the non behavioral subset. 3. Comparing the plots and performing a visual inspection of the difference of the ECDFs. 4. Performing a KolmogoroveSmirnov test (two-sided) to compare the two subsets. If the ECDFs are visually different and the maximum difference between the two cumulative distributions computed by the KolmogoroveSmirnov test is high, the analyzed simulation input is influential in driving the simulation output into certain regions. Furthermore, the EDCF gives insights on the mean, the maximal and the minimal value of the results for the complete set, the behavioral and the non behavioral subset. In this paper, the application MC filtering is shown for the combined BPS and CBA. 3. Simulation The main objective of this paper is the development of a UA methodology for a combined BPS and CBA rather than the presentation of a particular example scenario. The test case is supposed to illustrate the capabilities of the methodology. Furthermore, we wanted to analyze a realistic case, which at the same time is easy to explain and follow. It should allow the reader to focus on the proposed methodology rather than on modeling details. The central question of the presented case study is the cost effectiveness of different building design options. In order to answer this question, two scenarios are analyzed. The base case scenario is a conventional gas boiler that supplies the building with heat and domestic hot water. In a second scenario a solar thermal collector system in conjunction with a gas boiler provides the heat and domestic hot water. We analyze whether the additional investment for the solar thermal collector, the storage tank and the additional necessary equipment is cost effective or not. This requires only one set of MC simulations with the renewable plant equipment. This is possible because the total energy consumption of both design options is the same. The economic efficiency is assessed on the basis of the necessary additional investment for the solar thermal system (costs) and the savings based on the reduced gas consumption when a renewable system supplies a fraction of the total energy demand (benefits). The analyzed period of the CBA is 25 years. This period has to be chosen depending on the durability of the analyzed system and other project boundaries. It is assumed that the performance of the building and its energy demand is similar each year. Given the time frame of 25 years climate change and other variations for the BPS might play a role and could be analyzed with the proposed framework. However, this is not the focus of this paper and was not taken into account. 3.1. Building performance simulation model The simulated building is a typical German building with a net floor area of 436 m2. The model was introduced in Refs. [10,29]. There is no air-conditioning available in the building and the heat is emitted by radiators with thermostatic control valves. The building is equipped with sensors (outside temperature, heat meter, room temperatures etc.) to allow for a validation of the simulation. The envelope area to volume ratio of the building is 0.81 m2/m3, the mean U-value is 0.53 W/(m2 K), the total window area is 106 m2. As a zone model the simple hourly method (SHM) according to the ISO 13790 standard is used [30]. This zone model is based on five resistances and one capacity. The model was calibrated for this building; it showed good agreement with the measured room temperature and the heating demand of the building [29]. The

147

building is heated by a gas boiler. For this analysis, it is assumed that it is a multi-tenant residential building with 12 occupants and a solar thermal system. A solar thermal collector with an area of 25 m2 and a 2000 l storage tank are modeled. The collector model is implemented in Modelica using a model description by Isakson and Eriksson [31] The collector flow rate is controlled by an on/off controller, the collector flow rate per unit area is 0.004 kg/(s m2) and the storage tank is modeled as a simple one resistor/one capacitor (R1eC1) network. The solar thermal system is designed for domestic hot water and space heating. When the heat from this solar thermal system is not sufficient, a gas boiler meets the remainder of the load of the building. Further details on the model can be found in Ref. [10]. It is assumed that the air change rate (ACH) and the mass flow rate of the domestic hot water (m_dot_DHW) are uncertain. Furthermore, the number of occupants (occ) present at a particular time cannot be determined exactly. The mass flow of the solar liquid in the collector (m_dot_coll) is also varied during the MC simulation. This was done to determine, which flow rate leads to a good performance. Therefore, these four values (m_dot_DHW, ACH, occ, m_dot_coll) are varied in the MC simulation. The schedule for the domestic hot water flow rates is generated with a program which was developed in the Solar Heating and Cooling Program of the International Energy Agency (IEA-SHC), Task 26: Solar Combisystems [32]. The air change rates and the number of occupants are also implemented using a schedule. The variation of the schedule values is implemented either by multiplying a one-value MC-sampled scaling factor or adding a scaling summand. The distribution parameters of the uncertain inputs are listed in Table 1. The distributions and the standard deviation were chosen on the basis of the authors’ experience and literature. Pietrzyk and Hagentoft also use normal distributions for characterizing air change rates and the mean and standard deviation of ACH implemented in the simulations is comparable with cases illustrated in their study [20]. In the case of m_dot_DHW a standard deviation of 10% of the schedule values and a normal distribution seems appropriate for accounting for different consumption scenarios. For the scaling of m_dot_coll a the standard deviation of 20% around a typical design flow rate was intended to explore possible changes toward lower and higher flow rates. Occupancy is very specific to the analyzed case and the uniform distribution around the assumed schedule is intended to explore some possible scenarios. Note, that the values for the scaling summands for the number of occupants are rounded (integers). Fig. 3 illustrates the use of the scaling factors respectively the scaling summand for ACH and occ, The red line (for interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) represents the base case schedule and the box plots indicate the different values for ACH and occ used in the MC Simulation. In reality, there are causal relationships between some of the uncertain input variables (e.g., the number of occupants present influences the air change rate and the domestic hot water consumption). However, these relationships are neglected in this study. Table 1 Distribution parameters. Parameter

Distribution

m

s

ACH scaling factor m_dot_DHW scaling factor m_dot_coll scaling factor

Normal Normal Normal

1 1 1

0.2 0.1 0.2

Parameter

Distribution

Min.

Max.

occ scaling summand

Uniform

3

þ3

S. Burhenne et al. / Building and Environment 62 (2013) 143e154

8

Occupancy [--]

4

1.2

0

0.0

0.4

ACH [1/h]

12

16

148

0

2

4

6

8 10 12 14 16 18 20 22 24 Time [h]

0

2

4

6

8 10 12 14 16 18 20 22 24 Time [h]

Fig. 3. Variation of schedule values using scaling factors and scaling summands.

3.2. Cost-benefit analysis The quantification of the simulation input uncertainty in the building performance simulation was conducted using the expert knowledge of the authors and literature. For the CBA, a different approach was chosen for the inflation rate and the gas price. An analysis based on literature data was performed to quantify the input for a probabilistic uncertainty quantification [33]. This data was used to fit the ARIMA model and to predict the base case scenario for the development of the gas price and the inflation rate. This base case scenario is varied during the MC simulation. 3.2.1. Interest rate Depending on the specific project boundary conditions, different assumptions for the interest rate used in the CBA can be chosen. The biggest difference might be whether the money for any additional investment has to be borrowed from a bank or not. If the money comes from a loan with a fixed loan period and contract, the interest rate is not a subject to any uncertainty analysis. In this paper a fixed interest rate over the whole period was assumed. However, to quantify the influence of the interest rate on the result its percentage is varied in the MC simulation. The interest rate is sampled according to a uniform distribution in the interval [2,10], hence Xi w U(2,10) with i ¼ 1,2,.,N. Interest rates in this range seem reasonable for investments in renewable energy systems. 3.2.2. Investment cost During the design process prices for the investment are assumed. Depending on the market the final investment cost might

3.2.3. Inflation rate The inflation rate is the rate at which prices in an economy increase. As a basis for the analysis the inflation rates for Germany in the period of 1997e2011 were used [34]. The data of this period was used to fit an ARIMA model with seasonal component. For the ARIMA model parameter identification the statistical software package R was used [35]. The order of the model parts (p, q, d, P, Q, D) was chosen on the basis of a parametric run of all possible

0

2

4

6

Table 2 ARIMA model parameters.

-2

Inflation rate [%]

be different than the costs assumed in the design stage. To take this into account, the investment cost was assumed to be uncertain. Furthermore, the influence of the investment cost compared to the other economic input can be compared in this way. We just take the additional cost for the solar thermal system into account. The specific investment cost for one square meter collector area is sampled according to a normal distribution with m ¼ 700 EUR/m2 and s ¼ 50 EUR/m2, hence Xi w N(700,50) with i ¼ 1,2,.,N. Beside the investment cost, running costs, periodic replacement cost and maintenance costs are part of many CBA (specific information can be obtained from Ref. [18]). However, for simplicity reasons they are not considered in this paper. It is assumed that in the case of the analyzed solar thermal system these costs are relatively low compared to the other considered cash flows (e.g., the order of magnitude of maintenance cost is much smaller than the investment costs). Furthermore, the cash flow for replacement costs takes place 10e15 years after the investment took place and its influence on the NPV is reduced because of the compound interest effect.

2000

2010

2020 Year

2030

2040

Fig. 4. Time series for the inflation rate in Germany. The values from 1997 to 2011 are historical data. The values starting in 2012 are predicted using an ARIMA model that was fitted with the historical data. The dashed lines represent the confidence interval.

Parameter

Value

SE

Order of AR (p) Order of MA (q) Order of differencing (d) Order of seasonal AR (P) Order of seasonal MA (Q) Order of seasonal differencing (D) Coefficient AR (a1) Coefficient AR (a2) Coefficient MA (b1) Coefficient MA (b2)

2 2 0 0 0 1 0.8749 0.9899 0.7435 1.0000

e e e e e e 0.1278 0.0377 0.3058 0.3829

Result

Value

SE

log likelihood Akaike information criterion (AIC) Variance (s2) estimate

12.43 34.86 0.332

e e e

149

-2

Inflation rate [%] 0 2 4

6

S. Burhenne et al. / Building and Environment 62 (2013) 143e154

2012 2014 2016 2018 2020 2022 2024 2026 2028 2030 2032 2034 2036 Year Fig. 5. Box plot of the sampled values for the inflation rate used in the MC simulation. For illustration reasons the ARIMA prediction (red line) and the corresponding confidence interval (dotted red lines) were also plotted. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

combinations with p, q, P and D varied from 0 to 3 and the order of differencing varied from 0 to 2. This results in 2304 different combinations. As selection criterion the AIC value of the different model configurations was used. Not all possible combinations produce reasonable predictions and some result in error or warning messages of the ARIMA fit function in R. The reason for this might be limitations introduced by the time series data since the analyzed period is relatively small and only yearly values were available. However, the configurations with a low AIC were analyzed by visual inspection. The chosen configuration was among the best (i.e., low AIC) 11% of all possible model configurations. Fig. 4 shows the historic data, the prediction of the chosen model configuration and the confidence interval of the predicted time series. The prediction was performed for the next 25 years (2012e2036). The inflation rate has a periodic characteristic over time. The parameters of the model are listed in Table 2. To vary the base case prediction during the MC simulation, the standard error of the predicted time series is used in the sampling procedure. The variation is assumed to be normally distributed. This results in different scenarios within the MC simulation (particular time series of inflation realization to conserve the periodicity). However, the standard error of the prediction is very high and would overestimate the uncertainty in the opinion of the authors. This opinion is based on historical data and the fact that political measures aim to keep the inflation rate at a certain range. Therefore, the standard deviation in the sampling procedure is arbitrarily assumed to be 60% of the standard error of the ARIMA prediction. Fig. 5 shows box plots of the sampled values as well as the ARIMA prediction and its confidence interval. The box represents the data from the lower quartile to the higher quartile and the

3.2.4. Gas price Historic data from 1991 to 2011 [34] was used to fit an ARIMA model. The ARIMA model parameter identification was conducted similar to the one for the inflation rate. Fig. 6 shows the historic data, the prediction and the confidence intervals of the predicted time series. Unlike the inflation rate, the consumer gas price has an increasing trend rather than a periodic one. The configurations with the lowest AIC were analyzed by visual inspection. Similar to the ARIMA fit for the inflation rate not all possible combinations produce reasonable predictions and some result in error or warning messages of the ARIMA function in R. The chosen model configuration is among the best (i.e., low AIC) 34% of all combinations. The confidence interval indicates that the model outcome becomes highly uncertain when the future values are far away from the time when the investment takes place. The parameters of the model are listed in Table 3. It is not likely that the gas price is below zero although the confidence interval in Fig. 6 indicates that this is possible. The

0.05

0.10

0.15

0.20

Table 3 ARIMA model parameters.

0.00

Gas price [EUR / kWh]

median is also indicated. The length of the box is the interquartile range and the whiskers range from the lowest value still within 1.5 interquartile range of the lower quartile and the highest value still within 1.5 interquartile range of the upper quartile. The circles represent data outside this range (outliers). Predictions until 2036 are highly uncertain. The decision makers involved in the analysis would have to make the final decision which uncertainty should be applied in the analysis. The analyst should provide guidance in this process. This is true for every uncertain input. The on the basis of the data presented in Fig. 5 the average inflation rate Inflt is calculated according to Equation (2).

1990

2000

2010 2020 Year

2030

2040

Fig. 6. Time series for the consumer gas price in Germany. The values from 1991 to 2011 are historical data. The values starting in 2012 are predicted using an ARIMA model that was fitted with the historical data. The dashed lines represent the confidence interval.

Parameter

Value

SE

Order of AR (p) Order of MA (q) Order of differencing (d) Order of seasonal AR (P) Order of seasonal MA (Q) Order of seasonal differencing (D) Coefficient AR (a1) Coefficient MA (b1) Coefficient MA (b2) Coefficient seasonal MA (sb1)

1 2 2 0 1 1 0.6190 0.0649 0.8518 0.7051

e e e e e e 0.3238 0.6773 0.5917 0.5028

Result

Value

SE

log likelihood Akaike information criterion (AIC) Variance (s2) estimate

12.07 34.15 0.1966

e e e

150

S. Burhenne et al. / Building and Environment 62 (2013) 143e154

Fig. 7. Box plot of the sampled values for the gas price used in the MC simulation. For illustration reasons the ARIMA prediction (red line) and the corresponding confidence interval (dotted red line) is also plotted. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

4. Results and discussion In this paper a three-step analysis is performed: First, the results from the BPS are analyzed separately to show the advantages of the UA applied to BPS. In a second step, the uncertainties for the CBA are quantified. Finally, the results of the MC filtering are presented. 4.1. Building performance simulation results

Similar to the BPS, the convergence of the MC simulation for the combined BPS and CBA was analyzed. Fig. 10 shows the convergence of the estimate for the mean and the variance of the net present value. The mean estimate already converged at a sample size of 512 whereas the variance estimate needed a sample size of 1024.

0.8 ECDF

10

0.0

0 0.16

0

200 400 600 800 sample size

Fig. 8. Left: Convergence plot of the estimate of the mean value for the solar fraction. Right: Convergence plot of the estimate of the variance for the solar fraction.

0.401

0.2

5

PDF

4e-04

15

200 400 600 800 sample size

2e-04

variance (solar fraction) 0

0e+00

0.207 0.205 0.203

mean (solar fraction)

0.209

20

It is important to check the convergence of MC simulations after the simulation is performed. However, it is difficult to estimate the necessary sample size before a MC simulation because the

4.2. Cost-benefit analysis results

1.0

In the context of this paper the language and environment R [35] for statistical computing is used to generate the samples, change the simulation input file, call the executable of the simulation program, analyze, and visualize the results. The BPS models are written in the object-oriented and equation-based modeling language Modelica [36] and the executable that simulates the system is produced with the software Dymola 2012 [37].

0.6

3.3. Tools

convergence rate depends on the model structure, the number of analyzed parameters (dimensionality of the parameter input space) and the applied sampling strategy. Fig. 8 shows the convergence of the mean and the variance estimate for the solar fraction. The properties of the Sobol’ sequences require a sample size of N ¼ 2j where j˛Nþ [27]. The dashed vertical lines indicate these sample sizes. The simulation converged at a sample size of 512 since a higher sample size does not change the estimates of the mean and variance of the solar fraction. However, the results obtained with lower sample sizes might be also acceptable depending on the desired accuracy. The aim of the MC based UA is to obtain insights into the variance of the simulation results caused by the input uncertainty. Fig. 9 shows the histogram and the ECDF of the solar fraction. The solar fraction varies between approximately 15 and 30% and the probability that the solar fraction is above or equal to 20% is approximately 60% ((1  0.40)$100%). The solar fraction is just one example model result that can be analyzed. The analysis can also be based on other results such as the total energy demand of the building or the energy supplied by the solar collector.

0.4

uncertainty of the gas price time series is considered to be lognormally distributed. This results in different gas price scenarios. Similar to the assumptions for the uncertainty of the predicted inflation rates the standard deviation in the sampling procedure is assumed to be 60% of the standard error of the ARIMA prediction. Fig. 7 shows box plots of the sampled values as well as the ARIMA prediction.

0.20 0.24 0.28 solar fraction [--]

0.16

0.20 0.24 0.28 solar fraction [--]

Fig. 9. Left: PDF computed with the use of kernel density estimates on the basis of the simulation output. The gray area represents the probability that the simulation result is equal to or greater than the solar fraction threshold of 0.2. Right: Empirical cumulative distribution function (ECDF) of the result. The y-axis indicates the probability that the solar fraction is lower than a certain value.

S. Burhenne et al. / Building and Environment 62 (2013) 143e154

1.0e+07

2.0e+07

-8000

variance (NPV)

-7000

4.3. Monte Carlo filtering results

-10000 -9000

mean (NPV)

151

0

200 400 600 800 sample size

0

200 400 600 800 sample size

Fig. 10. Left: Convergence plot of the mean estimate for the net present value. Right: Convergence plot of the variance estimate for the net present value.

Fig. 11. Left: PDF computed with the use of kernel density estimates on the basis of the CBA results. The gray area represents the probability that the simulation result is equal or greater than the net present value threshold of 0. Right: Empirical cumulative density function (ECDF) of the result. The y-axis indicates the probability that the NPV is lower than a certain value.

The NPV of the additional investment varies between approximately 15,000 and þ20,000 EUR (Fig. 11). This massive variation is due to the significant influence of some analyzed parameters. The ECDF reveals that the probability that the investment has a positive NPV is approximately 9% ((1  0.91)$100%) under the chosen assumptions.

On the basis of the MC simulation results, a MC filtering analysis was performed. The first analysis is limited to the BPS. For the separation into the behavioral and non behavioral bins a solar fraction threshold of 20% was used. The bigger the differences between the ECDF for the behavioral and non behavioral result subsets, the bigger is the influence of the analyzed parameter in driving the simulation output into a specific region. Fig. 12 reveals that the air change rate (ACH) has a significant influence whether the solar fraction is above or equal to 20% or not. Lower air change rates lead to less heating energy demand and hence to a higher solar fraction. Other interesting pieces of information that can be obtained by analyzing the ECDF are the minimal and maximal values of a variable in the subsets. The highest scaling factor for the air change rate in the behavioral subset is approx. 1.2 whereas the lowest scaling factor in the non behavioral subset is approx. 0.8. The concept of the scaling factors is explained in 3.1 Building performance simulation model. If minimal values of the behavioral and the non behavioral subset are significantly different this indicates a very influential parameter. The same is true for the maximal values for both subsets. The scaling factor for the domestic hot water consumption (m_dot_DHW) is less influential. However, higher domestic hot water consumption leads to a better solar fraction since the solar thermal collector is able to cover the complete energy demand during the summer period. In the case of domestic hot water consumption, the minimal and maximal values are approximately the same. The number of occupants (occ) in the building influences the internal loads in the heated zone. Hence, many occupants result in less heating energy demand and therefore to a higher solar fraction. The ECDF has steps because of the occupancy scaling summand was rounded to integers. The flow rate of the solar liquid through the collector (m_dot_coll) influences the output slightly and the figure provides some possible design improvements, since a higher mass flow rate through the collector leads to a higher solar fraction in this case. Beside the visual inspection of the ECDFs, a Kolmogorove Smirnov test is conducted to compute the maximum distance between the behavioral and the non behavioral ECDF. Table 4 shows the results of the KolmogoroveSmirnov test. The test also

Fig. 12. MC filtering plot for the solar fraction. The threshold for the separation into the behavioral bin is a solar fraction of equal or greater than 0.2. All input parameters that lead to a solar fraction smaller than 0.2 belong to the non behavioral subset.

152

S. Burhenne et al. / Building and Environment 62 (2013) 143e154

Table 4 Results of the KolmogoroveSmirnov test for the solar fraction. Parameter

D

p-value

ACH scaling factor m_dot_DHW scaling factor occ scaling summand m_dot_coll scaling factor

0.6396 0.1619 0.3212 0.0994

<2.2e-16 5.026e-06 <2.2e-16 (approx.) 0.01553

reveals that the air change rate has the highest influence on the solar fraction. Another interesting result of the test is the p-value, which is a measure of the significance of the result. In the analyzed case the null hypothesis is that the results for the behavioral and the non behavioral subset come from the same distribution. In case the p-value is below the chosen significance level the null hypothesis is rejected and the two different subsets do not come from the same distribution. In this paper a significance level of 0.05 was

chosen and all p-values are below this threshold, which means that two distributions are significantly different and the analyzed simulation input has significant influence on the result. The MC filtering was also performed for the combined BPS and CBA. Fig. 13 shows the ECDFs for the different input parameters. Four input parameters for the BPS (ACH, m_dot_DHW, occ, m_dot_coll) and four parameters for the CBA (interest rate (IR), gas price (GP), investment cost (IC) and inflation rate (Infl)). In Fig. 13 the gas price and the inflation rate are represented by the sampled values obtained by sampling from a uniform distribution in the interval [0,1]. These values were the basis of calculating the values for the gas price and the inflation rate which are different in each year of the analyzed period (see 3 Simulation for details). Hence, the numbers for GP and Infl shown in the Figure can be interpreted similar to scaling factors. By comparing the figures it becomes obvious that the technical parameters do not influence the NPV significantly, whereas the economic parameters have a great influence. The

Fig. 13. Monte Carlo filtering plot for the net present value (NPV). The threshold for the separation into the behavioral bin is a NPV of equal or greater than 0. All negative values for NPV are in the non behavioral subset.

S. Burhenne et al. / Building and Environment 62 (2013) 143e154 Table 5 Results of the KolmogoroveSmirnov test for the net present value. Parameter

D

p-value

ACH scaling factor m_dot_DHW scaling factor occ scaling summand m_dot_coll scaling factor IR GP IC Infl

0.0491 0.109 0.0794 0.0898 0.5619 0.6891 0.1419 0.211

0.9901 0.2951 0.6903 (approx.) 0.5348 <2.2e-16 <2.2e-16 0.07838 0.001555

findings are: low interest rates (IR) lead to a positive NPV and for high gas prices a positive NPV is more likely. Both, the interest rate and the gas price are very influential parameters in the analysis. The probability of an interest rate below 6% in the behavioral subset is higher than 90%, whereas the same probability in the non behavioral subset is below 50%. In the case of the gas price the smallest value that leads to a positive NPV is approx. 0.4 (minimum value of the behavioral subset). The investment costs have a visible influence on the NPV even though the influence is less than the one for the gas price or the interest rates. High inflation rates lead to positive NPVs. The influence of the inflation rate is less than the ones for GP and IR. One reason for this is the periodic characteristic of the inflation rate that leads to high and low values in the course of the analyzed 25 years meaning that the influence compensates itself. Table 5 shows the results of the KolmogoroveSmirnov test for the combined BPS and CBA. Also in this case, the test reveals the same results, which can be confirmed by visual inspection of the ECDF plot in a quantitative way. The p-values for ACH, m_dot_DHW, occ and m_dot_coll indicate that given the chosen significance level of 0.05 the null hypothesis that both subsets come from the same distribution cannot be rejected. This is also true for the case of IC. The reason for this is that the number of samples that lead to positive NPVs is small compared to the number of samples that lead to negative NPVs. Another reason is the strong dominance of IR and GP. 5. Conclusion and outlook In the paper we propose a methodology to conduct an uncertainty analysis for a combined building performance simulation and cost-benefit analysis. Given the assumed uncertainties, the solar fraction in the example case varies from 15 to 30% and the probability, that the solar fraction is higher than 20%, is 60%. Furthermore, the Monte Carlo filtering revealed that the air change rate has the strongest influence on the solar fraction and which changes in the system can lead to the desired design outcome. One possible change is an increased flow rate of the solar liquid through the collector. The net present value of the solar thermal system can vary between 15,000 and þ20,000 EUR and the probability for a positive net present value is 9%. This probability will not be satisfactory in cases where the decisions are made on basis of monetary values rather than on environmental motives. Monte Carlo filtering can be used to determine possible changes in the design of the system. The method can be also used to compare different design options. In the example case it was revealed that the technical parameters have much less influence on the net present value than the economic parameters. The future gas price and the expected interest rate are the most influential parameters. Also a significant reduction of the investment cost can lead to positive net present values. However, compared the other economic parameters the investment cost is less influential than it was expected by the authors. This is due to the compound interest effect of the other parameters (gas price and

153

interest rate). This information could lead to the conclusion that policy makers who want to promote renewable energy systems should provide low-interest loans instead of subsidies for the systems. The results also reveal that the best promotion of renewable energy are high prices for conventional energies. The proposed analysis can supplement the typical design process and offer insight into the important aspects of the design and the role of the economic boundary conditions. The methodology can also be applied in supporting policy makers in providing support for renewable energy systems. The orders of the different parts of the ARIMA models used to predict future gas prices and the inflation rate were identified with an exhaustive search. This is computationally expensive since all possible combinations are computed. An optimization algorithm, which determines the optimal orders could improve the analysis. In the Monte Carlo filtering analysis single simulation parameters were analyzed, whereas interactions of parameters were not analyzed. Future research might focus on interactions in a Monte Carlo filtering setting. Furthermore, climate change and other variations for the building performance simulation during the analyzed time frame could be analyzed in a future study. 6. Abbreviations and notation Abbreviations and notation used throughout this article are given in Tables 6 and 7 respectively. Table 6 Abbreviations. Abbreviation

Description

AR ARIMA ARMA BPS CBA CDF ECDF HVAC LHS MA MC PDF SA SARIMA SHM UA UQ

Autoregressive Autoregressive integrated moving average Autoregressive moving average Building performance simulation Cost-benefit analysis Cumulative distribution function Empirical cumulative distribution function Heating, ventilation and air conditioning Latin hypercube sampling Moving average Monte Carlo Probability density function Sensitivity analysis Seasonal autoregressive integrated moving average Simple hourly method Uncertainty analysis Uncertainty quantification

Table 7 Notation. Symbol

Description

ai bi yt

Coefficient of the autoregressive part of ARIMA model Coefficient of the moving average part of ARIMA model Output of ARIMA model Noise term Air change rate Akaike information criterion Future cash flow KolmogoroveSmirnov statistic Gas price Investment cost Inflation rate Average inflation rate from year 1 to year t Nominal interest rate Number of inputs

3t

ACH AIC CF D GP IC Infl Inflt IR k

(continued on next page)

154

S. Burhenne et al. / Building and Environment 62 (2013) 143e154

Table 7 (continued ) Symbol

Description

L m_dot_coll m_dot_DHW N n NPV occ PDV t

Value of the likelihood function The mass flow rate of the solar liquid in the collector The scaling factor for the domestic hot water consumption Sample size Life cycle in years Net present value Occupancy Present discounted value Year at which cash flow takes place Mean value Standard deviation Year indicator

m s z

Acknowledgments This study was funded by the Reiner Lemoine Stiftung and the German Federal Ministry of Economics and Technology (BMWi) under the program ”EnOB” (BMWi 0327893A). References [1] Lomas K, Eppel H. Sensitivity analysis techniques for building thermal simulation programs. Energy Build 1992;19(1):21e44. [2] Macdonald I, Strachan P. Practical application of uncertainty analysis. Energy Build 2001;33(3):219e27. [3] Macdonald IA. Quantifying the effects of uncertainty in building simulation. University of Strathclyde; 2002. [4] Macdonald I, Clarke J. Applying uncertainty considerations to energy conservation equations. Energy Build 2007 Sep;39(9):1019e26. [5] de Wit S. Uncertainty in building simulation. In: Malkawi AM, Augenbroe G, editors. Advanced building simulation. Abingdon, UK: Taylor & Francis; 2003. p. 25e59. [6] de Wit S, Augenbroe G. Analysis of uncertainty in building design evaluations and its implications. Energy Build 2002;34:951e8. [7] Mara T, Tarantola S. Application of global sensitivity analysis of model output to building thermal simulations. Build Simul 2008 Dec;1(4):290e302. [8] Burhenne S, Jacob D, Henze GP. Uncertainty analysis in building simulation with Monte Carlo techniques [Internet]. In: Proceedings of SimBuild 2010, 4th national conference of IBPSA-USA. New York City: 2010 [cited 2011 May 21]. Available from: http://www.ibpsa.us/pub/simbuild2010/papers/SB10-DOCTS08B-01-Burhenne.pdf. [9] Burhenne S, Elci M, Jacob D, Neumann C, Herkel S. Sensitivity analysis with building simulations to support the commissioning process [Internet]. In: Proceedings of ICEBO 2010, 10th international conference for enhanced building operations. Kuwait City: 2010. Available from: http://www. icebo2010.org/. [10] Burhenne S, Jacob D, Henze GP. Sampling based on Sobol’ sequences for Monte Carlo techniques applied to building simulations. In: Proceedings of building simulation 2011: 12th conference of international building performance simulation association. Sydney: 2011. p. 1816e23. [11] Eisenhower B, Neill ZO, Fonoberov VA, Mezi c I. Uncertainty and sensitivity decomposition of building energy models. J Build Perform Simul 2011:1e14.

[12] Hopfe CJ, Hensen JLM. Uncertainty analysis in building performance simulation for design support. Energy Build 2011;43:2798e805. [13] Booth AT, Choudhary R, Spiegelhalter DJ. Handling uncertainty in housing stock models. Build Environ 2011;48:35e47. [14] Nas TF. Cost-benefit analysis: theory and application. Sage Publications, Inc.; 1996. [15] Pindyck RS, Rubinfeld DL. Mikroökonomie. 6th ed. Munich: Pearson Studium; 2005. [16] Rysanek A, Choudhary R. Multi-period optimization of building refurbishment decisions: assessing options and risk under economic uncertainty. In: Proceedings of building simulation 2011: 12th conference of international building performance simulation association. Sydney: 2011. p. 2447e54. [17] Jacob D, Burhenne S, Herkel S, Wagner A, Dodier RH, Henze GP. Comparing two methods of stochastic modeling for buildings. In: Proceedings of building simulation 2011: 12th conference of international building performance simulation association. Sydney: 2011. p. 1784e91. [18] DIN EN 15459. Energieeffizienz von Gebäuden e Wirtschaftlichkeitsberechnungen für Energieanlagen in Gebäuden [Energy performance of buildings e Economic evaluation procedure for energy systems in buildings]. German version of EN 15459:2007; 2008. [19] Saltelli A, Tarantola S. On the relative importance of input factors in mathematical models: safety assessment for nuclear waste disposal. J Am Stat Assoc 2002;97(459):702e9. [20] Pietrzyk K, Hagentoft C-E. Reliability analysis in building physics design. Build Environ 2008 Apr;43(4):558e68. [21] Corrado V, Mechri HE. Uncertainty and sensitivity analysis for building energy rating. J Build Phys 2009 Jun 19;33(2):125e56. [22] Ljung L. System identification: theory for the user. New Jersey: Prentice-Hall; 1987. [23] Box GEP, Jenkins GM, Reinsel GC. Time series analysis: forecasting and control. 4th ed. New Jersey: John Wiley & Sons, Inc.; 2008. [24] Sobol’ IM, Levitan YL. A pseudo-random number generator for personal computers. Comput Math Appl 1999 Mar;37(4e5):33e40. [25] Sobol’ IM. On the distribution of points in a cube and the approximate evaluation of integrals. USSR Comput Math Math Phys 1967;7(4):86e112. [26] Sobol’ IM. Uniformly distributed sequences with an additional uniform property. USSR Comput Math Math Phys 1976;16(5):236e42. [27] Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M, Tarantola S. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Comput Phys Commun 2010 Feb;181(2):259e70. [28] Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, et al. Global sensitivity analysis. The primer: John Wiley & Sons, Ltd; 2008. [29] Burhenne S, Jacob D. Simulation models to optimize the energy consumption of buildings. In: Proceedings of ICEBO 2008, 8th international conference for enhanced building operations; 2008. [30] ISO 13790. Energy performance of buildings e calculation of energy use for space heating and cooling; 2008. [31] Isakson P, Eriksson LO. MFC 1.0b. Matched flow collector model for simulation and testing. User’s manual. Stockholm; 1994. [32] Jordan U, Vajen K. Handbuch. DHWcalc. Werkzeug zur Generierung von Trinkwasser-Zapfprofilen auf statistischer Basis. Version 1.10 (2003). Lyngby/ Kassel; 2003. [33] Tsvetkova O. Uncertainty and sensitivity analysis for decision support in building design. University of Freiburg; 2011. [34] Entwicklung von Energiepreisen und Preisindizes [Internet]. 2012 [cited 2012 Mar 19]; Available from: http://www.bmwi.de/BMWi/Navigation/Energie/ Statistik-und-Prognosen/Energiedaten/energiepreise-energiekosten.html. [35] R Development Core Team. R: a language and environment for statistical computing [Internet]. Available from:, http://www.r-project.org; 2010. [36] Elmqvist H. Modelica e a unified object-oriented language for physical systems modeling. Simulat Pract Theor 1997;5(6). [37] Dassault Systèmes AB. Dymola. In: Dymola release notes. Library. Dynamic Modeling Laboratory; 2011.