Stochastic energy-demand analyses with random input parameters for the single-family house

In the analyses of the uncertainty propagation of buildings’ energy-demand, the Monte Carlo method is commonly used. In this study we present two alternative approaches: the stochastic perturbation method and the transformed random variable method. The energy-demand analysis is performed for the representative single-family house in Poland. The investigation is focused on two independent variables, considered as uncertain, the expanded polystyrene thermal conductivity and external temperature; however the generalization on any countable number of parameters is possible. Afterwards, the propagation of the uncertainty in the calculations of the energy consumption has been investigated using two aforementioned approaches. The stochastic perturbation method is used to determine the expected value and central moments of the energy consumption, while the transformed random variable method allows to obtain the explicit form of energy consumption probability density function and further characteristic parameters like quantiles of energy consumption. The calculated data evinces a high accordance with the results obtained by means of the Monte Carlo method. The most important conclusions are related to the computational cost reduction, simplicity of the application and the appropriateness of the proposed approaches for the buildings’ energy-demand calculations.


Introduction
The building energy demand is usually calculated based on the strictly defined values describing boundary conditions, material parameters or profiles of the occupants' behaviour. Such a deterministic approach is often applied in the engineering calculations, in which the worst-case scenario may be considered, i.e. to determine the maximum possible demand for the heating system sizing. Nevertheless, those calculations might cause substantial discrepancies between the predicted and the measured values, or misleading conclusions, whenever applied to the calculations of the annual energy demand or investigation of the occupants' thermal comfort. Numerous input parameters of such models are uncertain and should not be considered solely as strictly determined values, but rather as a set of possibilities with a certain probability. Otherwise, the discrepancies between the actual energy demand and the model-predicted values may be very significant, what may induce lack of credibility of such calculations (Coakley et al. 2014). It may be of particular importance in the energy retrofit analyses (Heo et al. 2012), in which verification and calibration of the energy demand should be performed with the reference to the observations. The energy savings and payback time of certain modernizations are strongly influenced by many parameters which are unknown (material parameters, predicted cost of energy in the lifetime of the building, etc.).
Several works indicate that impact of the input parameters uncertainty has to be investigated in order to improve the robustness of the mathematical models. Aude et al. (2000) presented a sensitivity analysis of the internal air temperature. Afterwards, the authors also analysed the uncertainty propagation of such parameters as solar irradiance, outdoor temperature, heating power or wind velocity. A comparison between the air temperature observed in the building and the values obtained using the calculations concerning the lower and upper uncertainty limits was described. Even though the differences between these series were significant, the observed temperature was within the uncertainty boundaries. On the other hand, Allard et al. (2018) considered differences between the calculated and measured energy performance of a single-family building.
Some of the most significant inputs in the analyses focused on buildings' energy-demand are the material parameters, such as heat conductivity, as well as external boundary conditions, especially exterior climate conditions (Allouhi 2015). The uncertainty of thermal conductivity has been thoroughly studied in the past. In Domínguez-Muñoz et al. (2010) the formulas describing the average conductivity and its standard deviation were proposed for numerous insulating materials, based on the available research data. Furthermore, the minimum and maximum values observed were provided for each type of material. Propagation of the thermal properties' uncertainty in the energy demand calculations has been also extensively investigated (de Wit and Augenbroe 2002;Catalina et al. 2008;Heiselberg et al. 2009;Prada et al. 2014;Belazi et al. 2018). In such calculations, one usually applied the typical meteorological year (TMY), which is defined based on the statistical elaboration of the long-time weather observations. However, it is often raised that such an approach is not representative with regard to the average energy demand (Kershaw et al. 2011) and may be further influenced by the climate changes, which are usually not considered in the standard calculations. Therefore, the uncertainty of the boundary conditions, especially external temperature, should be accounted for. Zhai and Helman (2019) studied a set of climate models in terms of the building energy demand. Impact of the climate conditions on the buildings' energy demand has been frequently studied over the past years (Heeren et al. 2015).
Among probabilistic methods, we can distinguish sampling and non-sampling approaches. In the energydemand related problems, the Monte Carlo (MC) or quasi-Monte Carlo methods, which belong to the former ones, are usually applied. The Monte Carlo method contains several sampling techniques. The basic random sampling method (BRS) requires many samples in order to obtain a satisfactory convergence; it has been used in, i.e., Prada et al. (2014) and Prada et al. (2018). In the Latin hyper-cube samplings (LHS) the domain is divided into subdomains of equal probabilities. An applicable number of simulations usually ranges from 100 up to 5 000 (Hopfe and Henses 2011;Eisenhower et al. 2012;Kim et al. 2014;Silva and Ghisi 2014;Yi and Braham 2015). However, this method has some drawbacks, e.g. to improve the accuracy of the LHS one has to increase the number of subdomains or the domain has to be divided again, which will cause overlapping of the subdivisions. Furthermore, the same algorithm of randomization applied to different parameters may result in a correlation between the subdomains. In the perturbation method, which is one of the non-sampling techniques, the random field is expanded by means of the Taylor series at the given order of expansion. In this case the number of simulations, which have to be performed, is much smaller than for sampling methods. The second-order perturbation analysis has been applied to the linear and nonlinear heat transfer problems Kleiber 1997, 1998) and afterwards employed by other researchers (Kamiński and Hien 1999;Wu and Zhong 2016;Kamiński and Strąkowski 2017;Yang and Cui 2017). Generalized perturbation method is applied by expanding the input variables and the state functions using the n-th order expansion of the Taylor series (Kamiński 2010). It has been also recently applied to analyses of the heat transfer (Kamiński and Strąkowski 2017;Ding et al. 2018).
One can also define another type of the probabilistic methods' categorization, depending on the data we can obtain from the analyses. Namely, the distribution, dispersion and reliability uncertainty analyses might be distinguished (Rivalin et al. 2018). Using the distribution analysis, one determines the full distribution of the random variable, the dispersion analysis provides mean, variance and central moments, whereas applying reliability analysis gives the probability of exceeding a threshold. In the energy-demand related problems, the first two types of investigation are mainly applied. The perturbation method, for general response function, is a dispersion analysis. The Monte Carlo method allows one to determine the distribution function.
This manuscript is aimed at an application of two methods, which may be used for propagation of uncertainties in the buildings' energy-demand problems, i.e. the stochastic perturbation method and the transformed random variable method. Statistical study of two parameters, the Styrofoam thermal conductivity and the external temperature, have been performed, in order to determine their probability density functions. The energy-demand calculations have been performed by means of the EnergyPlus software, for a reference single-family house of Poland. Afterwards, the response function of the energy demand has been calculated as a function of each of the examined parameters. The analyses have been performed for the two parameters independently (i.e. one-dimensional analyses) and for simultaneous variations of both parameters (two-dimensional analysis). The response function has been determined using the polynomial regression analysis. Uncertainty of the energy demand has been obtained applying the general perturbation stochastic method and the transformed random variables method. The general perturbation stochastic method has been applied to the heat transfer problems in some previous works. However, in the analyses of the buildings' energy demand, the Monte Carlo or the quasi-Monte Carlo methods are usually used. Meanwhile, the perturbation method allows to obtain the results with a much smaller computational cost. The transformed random variable theorem allows to determine not only some main parameters characterizing the energy demand distribution, but it provides its exact probability distribution function as well. According to the Authors' best knowledge, this method has not been previously used in the energy-demand related problems, despite the fact that its computational cost is much smaller when compared with the sampling probabilistic methods. Therefore, these two methods may be successfully applied in the analyses of the uncertainty distribution in the energy-demand calculations. In the manuscript, those two methods are compared and discussed. The methods themselves are known, but they have not been previously used in the energy demand-related problems. Meanwhile, using them, we may obtain the probability distribution function with a lower computational cost.

Distribution of the uncertain parameters
Thermal conductivity depends on numerous factors, such as density, local heterogeneities and inclusions, workmanship errors, conditions of production, moisture, temperature, etc. Most of these issues are independent (EN 10456 2007), which is why they are not investigated in this research. The uncertainty of the EPS thermal conductivity may be caused by experimental uncertainties (related to the measurement errors and procedure of determining the declared thermal conductivity), information uncertainties (caused by lack of data concerning the exact EPS type or inaccurate values given by the producers) or intrinsic uncertainties (which may occur due to differences of material properties within the material). In the works concerning the uncertainty of the material thermal conductivity, its distribution is usually assumed as normal, except for the sensitivity analyses, in which uniform distribution of either thermal conductivity or thermal transmittance may be imposed. In the presented research, the data published by the Construction Control Authority in Poland (https://www.gunb.gov.pl.probki) was analysed for extruded and expanded polystyrene for years 2016-2019. The results obtained in 313 measurements, with 4 samples for each test, have been investigated. The data is presented in the form of a histogram in Figure 1.
The observed values were in the range of 0.0292 W/(m·K) up to 0.0539 W/(m·K). Normal and lognormal distributions have been fitted to the data and compared by means of the Pearson's chi-square method. A better correlation has been obtained for the normal distribution with the mean of 0.0368 W/(m·K) and standard deviation of 0.0046 W/(m·K), thus such a distribution was assumed in further analysis.
The latter parameter, the uncertainty of which is analysed in this research, is the external temperature. It is very prone to variations and, simultaneously, affects the results profoundly. As mentioned before, the typical meteorological year (TMY) is usually applied in the calculations of the energy demand. It constitutes a set of hourly values, based on long-term observations, and used as a representative year for a given location. However, it strongly depends on the observation span and the random character of the weather. Therefore, it may be highly inaccurate, especially if the observations have not been sufficiently long. Furthermore, the climate undergoes constant changes, which may result in a variation of the temperature typical for a given localization. In the energy demand simulations, a typical reference year for Warsaw has been used. It was based on the observations performed in years 1971-2000. In case of TMY an average external temperature is equal to 8.26°C. In the calculations performed in the following sections, the temperature was defined for each hour based on the typical meteorological year determined for Warsaw. To investigate the impact of the external temperature uncertainty on the energy demand of analysed building, the following standard deviations were taken: σT = 0.82 °C and σ T = 4.0 °C. The first value reflects the standard deviation from the average annual temperature taken over 1971-2000 for Warsaw (https://danepubliczne.imgw.pl). Due to the global warming, the average temperature increases and the weather becomes more unsteady and unpredictable. It results in the growing number of violent phenomena, such as drought, blizzards, etc. To cover those phenomena the calculations were performed parallelly assuming a much higher value of temperature deviation, namely σT = 4.0 °C.

Mathematical model
In the manuscript, two methods are compared: the perturbation method and the transformed random variables theorem. The schematic flowchart of both methods is presented in Figure 2.

General n-th order stochastic perturbation method
Let us start with the formulation of the stochastic perturbation method assuming the single random parameter. Suppose that b(x) is the random variable and p(b(x)) is its probability density function. The random variable might be any material parameter, initial conditions, boundary condition including the ambient temperature or the heat transfer coefficient, etc. The random variable is defined by its two first probabilistic moments (expected value and variance): The integrals must be bounded. The n-th order perturbation method relies on the assumption that all variables as well as state functions, f(b), might be interpolated by the n-th order Taylor expansion using small perturbation parameter ε > 0 according to the formula: The k-th order variation (k = 1, ..., n) of b about its expected value is: Let us now introduce the Eq.
(2) into the Eq. (1a) giving the expected value of any state function with a small perturbation parameter using Taylor expansion: Placing Eq. (2) into Eq. (1b), and after performing analogous transformation as for expected value, one can obtain an analogous formula for the variance. The expected value of a function f(b) can be expanded using the tenth-order perturbation assuming general non-symmetric probability density function to the formula: denotes the k-th order central probabilistic moment of the random variable b, given by the equation: Applying similar derivation, the tenth-order expansion for the variance can be formulated:   6  2  5  3  4  7  7   2  4  3  5  2  6  1  6  8  8  2   1  8  2  7  3  6  4  5  9  9 2 5 4 6 3 7 2 8 2 1 9 10 10 Similarly, the central probabilistic moments of higher order might be calculated. It must be noticed that all odd moments about the mean (μ k (b) = 0, k = 1 + 2n, n Î N) vanish for the symmetric probability density function. Consequently, the formula defining both expected value and variance must be modified appropriately.
The derivation becomes more complicated when multiple random variables are concerned. Here, we limit the case for two random variables, (b 1 , b 2 ), with symmetric density functions. The second order estimation of the f(b 1 , b 2 ) expected value and covariance matrix S x of 2 random variables is presented in Appendix A (Eq. (A1), Eq. (A2)), which is available in the Electronic Supplementary Material of the online version of this paper.
In the stochastic finite element method, the critical issue is to determine the partial derivatives of the response function concerning random parameters. In the direct differential method, introduced by Hien and Kleiber (1997Kleiber ( , 1998, the partial derivatives are calculated directly from the governing equation. However, in this case, the direct differentiation of the equation describing energy consumption would be difficult. To overcome this obstacle one can calculate multiple solutions of the energy demand around the expected value of a random parameter, i.e. heat conductivity and external temperature. Subsequently, by means of various approximation techniques, one can formulate the response functions. The procedure of the response function determination has been presented in Appendix B, which is available in the Electronic Supplementary Material of the online version of this paper.

Transformed random variables theorem
Let us introduce the technique, which allows one to approximate the probability density distribution of any random variable function. It is based on the well-known transformed random variables theorem (Pugachev 1984;Krysicki et al. 2007). Theorem 1 Let b Î (x 1 , x 2 ) be the random variable with the probability density function p b (b). Suppose that the new random variable y is the function of b, namely y = f (b).
, then the following equation defining the probability density function, p y (y) for y (y 1 , y 2 ), holds (Pugachev 1984;Krysicki et al. 2007): and p y (y) = 0 elsewhere. Usually, in a practical application f Î C ¥ , which allows one to expand the function f in Taylor series. The main technical difficulty, which arises during the application of the equation is finding the exact form of the inverse function h = f −1 . Let us show the application of the above theorem assuming that f depends on a single random variable and that f is the function of two random variables. By introducing the perturbation Eq. (2) and assuming the perturbation parameter equal to unity the random variable function f might be expressed in any given point b: where b is the expectation of a random variable b. One can notice that Eq. (10) resembles the Taylor series derived for function f. Let us discuss some special cases, i.e. by truncating the terms with order higher than two, one gets the equation: The similar approximation might be performed for higher order polynomials. Assuming that C = 0 the transforming function, f, becomes a linear one. Having normally distributed random variable b with the expectation b and the variance var(b) = μ 2 (b), the new random variable y = f(b) is also normally distributed with the expected value y A Bb = + and the variance var(y) = B 2 var(b). However, in most practical cases, the simplification that f is linear becomes too strong. Therefore, polynomials with higher order have to be taken into account. On the other hand, the difficulty of the inverse function h = f −1 determination grows along with the rise of polynomial order. The assumption of ( ) 0 f b * ¢ ¹ for any b * Î (x 1 , x 2 ) might be overcome by splitting the interval (x 1 , x 2 ) into two parts ( Subsequently, the theorem can be straightforwardly applied for both subintervals. Now let us consider the function of two random variables, namely y = f(b 1 , b 2 ). Assuming that the input variables are independent, the probability density function reads: We can formulate two alternative versions of the Theorem 1. By applying the theorem on the marginal probability density function, one can obtain the expression defining the probability density distribution of the random variable, y, being the function of two random variables b 1 , b 2 , which might be expressed as: with respect to the first variable, the second variable is assumed to be the fixed parameter, and y is a Jacobian, which may be generalized for any given number of variables (Pugachev 1984;Krysicki et al. 2007). Alternatively, it can be given as: with respect to the second variable. To complete the task, one has to apply the Taylor expansion of the bivariate function: , ,

Characterization of the analysed building
The national residential sector of Poland consists of single-(SFHs) and multi-family houses (MFHs) with SFHs standing for a majority (approx. 90% of inhabited buildings). Energyrelated analyses of residential buildings can be performed either for a specific object or a group of them, constituting a cluster. It is recommended to use specialized software, which includes dynamic phenomena (especially heat transfer and heat accumulation/release). Those advanced analyses are typically performed using computational simulations with at least an hourly calculation step, allowing an accurate assessment of the investigated issue.
The object of performed analysis was one of the representative buildings of the Polish residential sector. The residential stock of Poland was analysed in a specialized report named TABULA (Typology Approach for Building Stock Energy Assessment) (NAPE 2012). The major assumption of the TABULA report was to develop a residential building typology of Poland. The main purpose of the proposed typology was to present a general overview of the Polish buildings' energy performance by means of the statistical elaboration for the total amount of over 60 thousand existing single-family houses.
In the presented study, the analysis was performed for a representative SFH of Poland, which is composed of the most typical elements for the defined building type. Objects determined as representative SFHs can be also named as "average buildings" and successfully used to assess the results of the application of the modernization for existing houses. The TABULA report specifies 7 types of Polish SFHs, characterized by construction periods. Out of the examined Polish typology, the representative SFH constructed after 2009 was selected for further investigations.
According to NAPE (2012) most of the construction data is known, nevertheless, some assumptions have to be made in order to define a valid computer model of the analysed building. The selected SFH has traditional domestic heating (DH) and domestic hot water (DHW) systems and natural ventilation (NV). The DH system is relatively simple: it includes a heating stove (supplied by gas), supply piping system and heater units (water radiators). The DHW system was not defined in the performed model due to the fact that it was not under consideration for the conducted study. The buildings have two occupied floors and an unoccupied attic with a gable-type roof (30 degrees sloped). On each floor, there is a bathroom, with area equal to approx. 10% of the floor area. Total occupied area is equal to 188.49 m 2 , with a volume of 540.00 m 3 . The construction consists of red-clay brick walls, concrete slabs and wooden roof structure covered with tilling, while windows are double-glazed. The geometry and HVAC installation schema of the analysed building is presented in Figure 3.
The energy profile of the analysed building is defined by the following parameters: As it was already mentioned, some additional assumptions were required to define a valid building model in the EnergyPlus software. Those assumptions concerned parameters which were not included in the TABULA report. The most important ones related to building's shape (square-based, 10.48 m × 10.48 m), numbers of residents (family of four), and their activities (mid-level activityresting), operation timetables (typical working schedulefrom Monday to Friday, from 8 AM to 4 PM), design interior air temperature (22 °C in living zones and 24 °C in bathrooms) as well as lighting system (LEDs usage) and equipment design loads. The expanded polystyrene foam was used as insulation material of exterior walls, slabs and ground floor, while mineral wool was used in roof. Using the above-mentioned assumptions and the data published in the TABULA report it was possible to define a detailed building model in the EnergyPlus software in order to perform energy-related analyses. The EnergyPlus software is probably the most popular and the most widely used whole building energy simulation programs (Crawley et al. 2001). The analysed building is characterized by good thermal parameters of a building enclosure, with efficient DH and DHW systems.
That building is the closest one (out of all 7 representative ones) to fulfil the actual Polish regulations (regarding thermal properties of building envelope). Also, it was chosen because it is a representative object for many existing buildings in Poland, as well as due to its data availability. In Poland, more than 80% of buildings have insufficient thermal insulation of building enclosures (IES 2017). Results can be considered as a reliable output for a selected group of a residential archetype of Poland.
During the performed study, parametric analyses of two factors (thermal conductivity of expanded polystyrene and dry bulb temperature) were made (see more in Chapter 2). Therefore, one more additional assumption was made due to the variability of a dry bulb temperature (DBT). Heating system operation is available for the whole year, 24/7 in order to ensure heating whenever it is necessary due to variability of outside temperatures.

Energy demand analysis
The energy consumption was calculated using the EnergyPlus software for the single uncertain parameter (random thermal conductivity -case 1, random dry bulb temperature assuming σ T = 0.82 °C -case 2 and random dry bulb temperature assuming σ T = 4.0 °C -case 3) and for both parameters considered as uncertain (random thermal conductivity and external temperature with σ T = 0.82 °C -case 4, random thermal conductivity and external temperature with σ T = 4.0 °C -case 5). The range, where the response function was approximated was equal to , [ 3 , 3 ] λ T μ σ μ σ -+  .
Assuming such a wide range it allowed one to conduct more accurate determination of nonlinear approximation function. Unfortunately, on the other side, one might get unrealistically low or high value of uncertain input parameters, i.e. too low thermal conductivity. The relation between the energy consumption and the random parameters has been determined by means of the polynomial regression. Comparison of the EnergyPlus results with the ones obtained using the approximation with a polynomial of various orders is presented in Figures 4 A, B, C for cases 1, 2 and 3, respectively. One can observe that for case 1 the approximating functions are nearly linear. For case 2, increasing the polynomial order induces slight improvement of the curve fitting, while for case 3 strongly nonlinear relation was derived. For bivariate problem (cases 4, 5), Figure 5, the results obtained using the 2 nd order of polynomial for thermal conductivity and the 6 th order for temperature are presented only to improve the clarity of the presentation. Assuming the smaller temperature deviation (case 4) one obtained more linear approximation, whereas for higher temperature deviation (case 5) strong nonlinearity can be observed. Obviously, the energy demand increases along with the increase in thermal insulation conductivity and decreases along with the growth in external temperature, since no building cooling has been accounted for. As one can notice, good accordance was obtained.

Stochastic perturbation method analysis
Let us notice that the input variables are independent, Fig. 4 Comparison of the yearly energy consumption for heating EH obtained using the EnergyPlus simulations and as a result of the polynomial approximation of different orders for random thermal conductivity, case 1 (A), and external temperature, case 2 (B) and case 3 (C) hence cov(b 1 , b 2 ) = cov(b 2 , b 1 ) = 0. It leads straightforwardly to the simplification of Eq. (A2) in the Appendix. The expected value and standard deviation of the final energy demand were calculated using the general n-th order stochastic perturbation method and has been presented assuming different orders of the polynomial approximation in Table 1 for cases 1-3 and in Table 2 for cases 4-5. For each analysis, the values obtained using the polynomial approximation were compared with the EnergyPlus results. The approximation error, ES, of polynomials drawn in Figure 4 and Figure 5 has been calculated as: is the approximation error in point i, E pol,i is the value Fig. 5 Comparison of the yearly energy consumption for heating EH obtained using the EnergyPlus simulations (black orthogonal frame) and the polynomial approximation (green surface) for the bivariate problem, case 4 (A) and case 5 (B), respectively  obtained by polynomial approximation, E EP,i is the value obtained using the EnergyPlus software and k is the number of all points. For random thermal conductivity (case 1) the approximation error decreases from 0.25% to 0.01% by using the 2 nd polynomial order instead of the linear form, while the expected value is changed by 0.2% with reference to the basic value. For random thermal conductivity with lower temperature variance, the approximation error decreases from 0.51% to 0.03% by using the 6 th order polynomial. Application of higher polynomial orders induces only small improvement of the results, while it induces higher complexity of further calculations. For case 3, with higher temperature variance, it can be noticed that only for the 6 th order polynomial sufficiently accurate results have been obtained, with approximation error below 0.1%. Finally, the second-order was assumed for case 1 and sixth-order polynomials for case 2 and case 3 as the best approximation of energy consumption. Similar polynomial orders were analysed for the bivariate problems. One can notice that using the 2 nd order for thermal conductivity and 6 th order for external temperature the best accuracy of the results (approximation errors are equal to 0.15% and 0.08% for case 4 and case 5, respectively) may be obtained. Therefore, it has been decided that such an error is acceptable, while application of higher orders would increase complexity of further calculations (i.e. determination of the inverse functions to the response function). It should be noted that the stochastic perturbation method allows one to determine the expected value and central moments for linear and non-linear problems, regardless of whether one applies direct differentiation of governing equation or the response function as the output variable.

Transformed random variables method
The comparison of the probability density function and the cumulative distribution functions obtained using the proposed transformed random variables method, assuming the linear form of function f with the one obtained for a 6th-order polynomial is presented in Figure 6 for random thermal conductivity (case 1) and random temperature (cases 2 and 3). The thermal conductivity uncertainty (case 1) has much smaller impact even when compared to the case 3 with lower variance of external temperature (case 2). Due to the high approximation error of the linear response function in case 3 (see Table 2), it can be noted that large differences may be observed in the PDF and the CDF between the linear relation and the 6 th order polynomial response function for this analysis. Figure 7 presents the distribution and cumulative distribution functions for bivariate problems -case 4 and case 5, respectively. Again, for case 5, in which larger variance of temperature has been assumed, assumption of the 6 th order polynomial response function instead of the linear one changes the results significantly. It can be noticed that the nonsymmetrical distribution of energy demand, assuming a symmetrical normal distribution of input parameters, was obtained for the cases 3 and 5 - Figure 6 and Figure 7. The expected value and standard deviation obtained using this method for different polynomial order is presented also in Table 3 and Table 4 for cases 1-3 and 4-5, respectively. To assess the impact of the parameters' randomness on the response function, the quantiles of energy demand, defined as: for arguments equal to q = 0.8, 0.9, 0.95 were investigated for first-order and higher-order polynomials. The values of the quantiles are given in Table 5 for the cases 1-3. By application of polynomial response functions of higher order, the Q E (0.95) is modified by 0.04%, −0.11% and −0.05%   for cases 1, 2 and 3, respectively, with reference to the results obtained by means of the linear approximation (see Table 5). For bivariate problem, a change by −0.11% for case 4 and 1.71% for case 5 can be noted (see Tables 6 and 7). The results were compared with the results obtained using the Monte Carlo method. The analysis was performed for every single random parameter, in which either 1 000, 10 000 or 100 000 simulations were conducted. The distributions of input parameters were taken from Section 2: N(μ, σ) = N(0.0368 W/(m·K), 0.0046 W/(m·K)) for the random thermal conductivity (case 1), and N(8.26 °C, 4.0 °C) for random mean ambient temperature (case 3). The histograms of the energy demand are given in Figure 8A and Figure 8B for the case 1 and the case 3, respectively. The Monte Carlo method is a technique, which allows us to determine not only central moments, but the distribution itself as well. However, one can notice that even though a large number of simulations have been used in both analyses, the final results may differ - Figures 8A, B. The expected values and the standard deviations are presented in Table 8.
The proposed transformed variable and stochastic perturbation methods allow one to determine the uncertainty of energy demand with a much smaller number of runs than in case of the Monte Carlo method. For the MC method, in order to obtain robust results, the number of the simulations has to be very high -at least 10 000 simulations for the presented cases of single uncertainties, and even at this number accuracy may not be sufficient. Simultaneously, the proposed transformed variable method provides explicit distribution function with only 11 runs of the EnergyPlus software for a single random variable and 121 runs for a bivariate problem. This number depends on the shape of the response function. For strongly non-linear relation or non-continuous functions, this number should    It should be noted that it is possible to apply polynomial fitting formula also in the Monte Carlo method. If the response function fitting formula allows to obtain accurate results, it is possible to determine the result for each input parameter value according to the response function instead of performing time-consuming calculations for a large set of the input parameter values. If the response function fitting formula is defined based on a wide range of the input parameter, the probability of getting a value from outside this range is very low in a draw of a single value. Assuming that the range of [μ − 3σ, μ + 3σ] was used in the polynomial regression analysis, this probability would be equal to appr. 0.27%. However, depending on the number of the Monte Carlo individuals, this probability increases accordingly, what causes the possibility of extrapolating the values of the response function, what should be avoided. Therefore, additional conditions, which would perform additional calculations in such a case, would be required, or truncated normal distribution should be assumed for the Monte Carlo sampling. In the perturbation method and the transformed random variables method, such transformations are not required.
What is more important, the Monte Carlo method approximates the PDF. The results we obtain are strongly dependent on the sample size. Even though the perturbation method and the transformed random variable theorem also depend on the grid of points used to obtain the response function, they allow to obtain the exact PDF. Further, the number of the bins in the histogram will not impact the results -which may be problematic in the analyses of the Monte Carlo method results. In Figure 8 in the manuscript it is visible that even for a large sample size in the Monte Carlo simulation, the irregularities in the histogram were observed, which was partly related to the preciseness of the input data in the EnergyPlus software. In case of the transformed random variable theorem, we may also determine the PDF and the CDF.
In both methods analysed, the computational cost is dependent on determination of the response function. Determination of the input uncertainty distribution and the response function is performed by means of the polynomial regression analysis (least square method). Using higher polynomial order, the approximation error is limited. Similarly, by application of a denser grid of points in the range of the uncertain input parameters, we may increase the accuracy. These processes also induce larger computational cost. Number of the points in the grid is not of paramount importance for the calculation time when two-dimensional analyses are performed. However, if this method was applied to more complicated problems, exceeding this number could induce severe computational requirements. Therefore, it should be thoroughly investigated before analysing any particular problem. It has been decided after initial calculations to limit the grid of points to 11 for onedimensional and to 121 for two-dimensional analyses. In the presented research, increasing this number did not change the results by more than 1%. Higher polynomial order also induces larger requirements in both methods, especially for determination of the inverse function. Furthermore, even though increase of the polynomial order could cause further reduction error, it was decided that if this error was about 0.1%, it was considered as negligible when compared to other uncertainties in the calculations, therefore higher orders of the polynomial have not been considered. Determination of the response function took appr. 3 seconds to find the optimal polynomial coefficients. The time necessary to perform calculations according to the perturbation method and the transformed random variables theorem was equal to appr. 35 seconds and 60 seconds, accordingly. Meanwhile, the calculations according to the Monte Carlo method, with 100 000 runs, took approximately 5 hours.

Conclusions
The impact of the material conductivity and external temperature uncertainty on the building's energy consumption was analysed. The input variables were assumed to be independent. Distribution of the uncertain parameters was based on statistical analyses of the material and weather databases. Calculations of the energy demand were performed using the EnergyPlus software.
Two techniques belonging to the disperse and the distribution methods of uncertainty analyses were comparedthe general n-th order perturbation method and the transformed random variables method, respectively. The statistical analysis was based on the polynomial approximation of energy consumption determined using EnergyPlus program considering random input parameters. 1) The Monte Carlo method is often used for its simplicity.
In case of intractable problems, it may allow to obtain results in an easy manner when compared with other approaches. It is also easy to apply in majority of the cases. It is a distribution approach -it provides data concerning not only central probabilistic moments, but also approximation of the probability distribution function. However, this method has some drawbacks as well. The computational cost of this method is considerable, when large number of individuals is analysed. If number of simulations is too small, it may, from the other hand, provide inaccurate results. 2) The disperse method -stochastic perturbation methodallows one to calculate the expected value and central moments of any desired order. The entire probability distribution may be obtained only for the linear form of the response function. In case of applying low polynomial order, it may be applied to problems characterized by low coefficient of variation. For higher coefficients of variation, higher order expansion of the equation should be applied, what also increases computational cost, so order of the polynomial should be determined carefully. This method requires calculation of the derivatives of the state function with respect to the random variable. It may be obtained by means of the direct differential method or the response function. The first approach may be applied in case of uncomplicated mathematical models, in which we may obtain derivatives of the equation. For more complicated models, the response function is determined by polynomial regression analysis for a set of points in a range of the input parameter and afterwards derivatives are obtained. This approach requires only a fraction of the computational cost necessary to perform robust Monte Carlo simulations. 3) Application of the transformed random variable method allows to obtain explicit form of the probability distribution function. The output functions have to be monotonic and of C1 class in order to enable determination of the inverse functions, which are necessary in this method. If requirement of the invertible function is not fulfilled in some isolated points, the domain must be split into subdomains, where the output functions are invertible. 4) It is shown that approximation error diminishes along with the rise of polynomial order. The expected value and standard deviation change less than 0.1% by increasing the polynomial order from 1 to 2 assuming thermal conductivity as the single random variable. For the second considered random parameter (external temperature) the statistical parameters change by around 9% and 3%, appropriately, after rising the polynomial order from 1 to 6. Assuming bivariate problem (thermal conductivity and external temperature together) the increase in polynomial order from 3 to 8 causes change of the expected value and standard deviation by about 0.2% and 2%. 5) The distribution method -transformed random variable method -provides the exact formula for probability density function for the single and bivariate cases, appropriately. As shown in Figures 5, 6 and 7, the accuracy of the method strongly relies on the applied approximation technique. Having the normal (symmetric) distribution function one might obtain the nonsymmetrical distribution of the output variable, depending on the shape of the transformation function. Assuming a normal distribution of the input variable and a linear transformation function one gets the normal distribution of the output parameter. Having the exact form of the distribution function, one can calculate any q-quantile of the output variable. 6) The stochastic perturbation method is very effective for estimating the parameters of the probability distribution function. It might be easily applicable for most of the practical problems. The application of response function allows one to avoid the necessity of calculating the partial derivatives of the output variable directly from the differential equation. The transformed random variable method requires more sophisticated mathematical apparatus but provides the explicit form of output variable distribution function. The second method is more general, but its application demands the profound understanding of statistical analysis. 7) It must be noticed that for dependent input variables the application of the stochastic perturbation method is quite simple. The covariance between input random variables must be introduced into the equation defining the expected value and variation. Unfortunately, such a modification would be much more difficult for the second method. 8) The parameters obtained using both stochastic methods were verified against data calculated using the standard MC method for 1 000, 10 000 and 100 000 runs. The differences between the MC and stochastic methods are acceptable; however, for a nonlinear problem, the difference is larger than for linear one.
Open Access: This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/