A robust optimization model for designing the building cooling source under cooling load uncertainty

Applied Energy 241 (2019) 390–403

Contents lists available at ScienceDirect

Applied Energy journal homepage: www.elsevier.com/locate/apenergy

A robust optimization model for designing the building cooling source under cooling load uncertainty

T

Jide Niua, Zhe Tiana, , Yakai Lua, Hongfang Zhaoa, Bo Lanb ⁎

a

School of Environmental Science and Engineering, Key Laboratory of Efficient Utilization of Low and Medium Grade Energy, MOE, Tianjin University, Tianjin 300072, China b School of Architecture, Tianjin University, Tianjin 300072, China

HIGHLIGHTS

range of uncertain parameters is investigated by Monte Carlo simulation. • AThewide method can be conducted with computational efficiency. • The robust proposed robust method guarantees the global optimization. • Effectiveness, and accuracy of the proposed robust optimization method are demonstrated based on a case study. • A comparison efficiency between robust and deterministic method is conducted. • ARTICLE INFO

ABSTRACT

Keywords: Robust optimal model Cooling load uncertainty Monte Carlo simulation Cooling source design

The cooling load is the basis of cooling source design, while the inherent uncertainties of building information, weather condition, and internal load make it impossible to obtain a determinate load. Probability method can represent characteristics of uncertain cooling loads well, but a large number of Monte Carlo (MC) simulations should be carried out. The optimal cooling source design can be formulated as a mixed-integer-linear-programming model (MILP), which can be solved efficiently to obtain the global optimality using a state-of-the-art MILP solver. However, if all the MC simulations are used in the optimization problem, the size of MILP model would lead to computational issues or even be unsolvable. This study, therefore, proposes a robust optimization design model for sizing the cooling source when there is cooling load uncertainty. A method named in this paper cooling load bin is proposed and implemented by converting time series cooling loads obtained by MC simulations to those in the frequency domain. Therefore, the cooling load frequency and mean value in each cooling interval are obtained and used in the optimal design model, which can be solved efficiently by the General Algebraic Modeling System (GAMS). The robust model is applied to the optimal design of a cooling source to minimize the cost. A case study on the cooling source of a hospital in Tianjin is conducted to demonstrate the effectiveness of the proposed robust model. Furthermore, the accuracy of the solution and the computational efficiency used to evaluate the proposed robust model are systematically investigated and compared with the deterministic model.

1. Introduction An effective design of building cooling source could improve the energy performance of the building significantly. However, balancing the demand and supply is a challenging task because of the instantaneously varying load and limited operational flexibility of system. In the conventional method, the peak cooling load is used to configure the cooling system [1]. Safety factor or redundancy design are often used to enhance the reliability of the cooling system by considering the ⁎

uncertainties of weather data, building information, and the internal load [2]. These approaches can make the cooling system fulfill any uncertain cooling load demand, but at the same time make the cooling system oversized and also lead to an efficiency reduction if the operational conditions deviate from the design conditions [3]. Although hourly cooling loads for representative days are developed and implemented in a mathematical optimization model for designing an energy system [4–6], this method can only obtain a sub-optimal cooling system design because it does not quantify demand uncertainty.

Corresponding author. E-mail address: [email protected] (Z. Tian).

https://doi.org/10.1016/j.apenergy.2019.03.062 Received 8 August 2018; Received in revised form 17 February 2019; Accepted 7 March 2019 0306-2619/ © 2019 Published by Elsevier Ltd.

Applied Energy 241 (2019) 390–403

J. Niu, et al.

Nomenclature

MARD r L

Tcost In Op N Ca Up CRF f Ele Ep c

annual total cost (yuan) investment (yuan) operational cost (yuan/year) the number of chiller the capacity of a chiller (kW) unit price (yuan/kW) capital recovery factor frequency of cooling load chiller power (kW) electric price (yuan/kWh) a symbol used to count how many cooling loads fall into bin i Y times of MC simulations q the cooling load of a chiller (kW) cload cooling load (kW) bin binary variable a, b, c, d coefficient of chiller COP COP coefficient of performance of the chiller M Big-M, a constant q The cooling load of a chiller (kW) seq SOS2 variable ARD absolute relative deviation of the objective k the width of CLB, kW

moving absolute relative deviation of the objective the discount rate the lifetime of the chiller

Superscripts T D P s ¯ i cload

type of chiller, T = A, B, C deterministic model probability optimization model s-th MC simulation the mean of cooling loads in bin i

Subscripts i j k I

i-th CLB j-th chiller initial value of moving total cooling load interval

Greek symbols α β

Moreover, accurate demand forecasting is almost impossible because of the inherent uncertainty caused by the stochastic nature of weather conditions, occupant behaviors, etc. Therefore, it is necessary to develop a better mathematical optimization model for designing building cooling source with uncertain cooling loads [5,7]. A number of studies have been performed to investigate uncertainties in design of building energy system. Two crucial problems were widely investigated: (1) how to quantify the uncertainty and (2) how to establish the mathematical optimization model for the design of cooling source while considering uncertainty.

correction factor of the investment correction factor of operational cost time (h)

enough accuracy for both of them was achieved after 780 MC simulations. Besides these external parameters, i.e., the load, solar radiation, and wind velocity, the models used to configure energy systems is also have an inherent uncertainty about the physical process. For instance, Li et al. quantified the uncertainties in pressure head calculations and chilling water system construction using normal and uniform distribution functions [15]. Similarly, Cheng et al. performed MC simulations to generate the hydraulic resistance distribution by quantifying the uncertainty [16]. Overall, because the MC simulation based on the probability distribution function is intuitive and reliable, it is widely implemented to treat these uncertainties occurring in the building load calculation and energy system design. However, Tian et al. argued that it is hard to obtain detailed information for some parameters, especially at the initial stage of a building project [17]. Moreover, the total number of uncertain parameters treated by MC simulations would lead to overly large problem sizes and computational tractability issues [8]. Therefore, some non-probability approaches were proposed to account for the uncertainties. For instance, a non-probabilistic Dempster-Shafer theory of evidence approach is applied for energy consumption uncertainty analysis in a green building [17]. Fu et al. proposed an information entropy approach to quantify uncertainty in an integrated energy system; the dependence relationships were further assessed by the mutual information and correlation coefficient [18]. In [19], the delta technique was applied for constructing prediction intervals for future loads instead of forecasting their exact values. Furthermore, Khosravi et al. proposed a new technique, the lower upper bound estimation (LUBE) method, for constructing prediction intervals and it was demonstrated that the LUBE method is effective and reliable [20]. The above literature review on approaches to uncertainty quantification approaches can be summarized as follows: the probability distribution function combined with MC simulations is a suitable method to describe a variable’s uncertainty, and yet computational issues arise. A non-probability approach is easily implemented, as less variable information is needed, but it results in reduced accuracy. In this study, the probability distribution function is used to quantify variable uncertainty and a new measure to treat the computational issues is introduced, which is detailed in the following section.

1.1. Methods for uncertainty qualification As stated in [8] “The mode-based design process of a distributed energy system (DES) is irrevocably affected by uncertainty”. Under uncertain circumstances, it is of primary importance to identify and quantify uncertain parameters [9]. The Monte-Carlo (MC) simulation combined with the probability distribution function is a common method for addressing uncertainty. Zhou et al. assumed that each type of load followed a normal distribution and that 1.96 equaled 0.2µ [7]. Similarly, in [10], the uncertain parameters, i.e., solar radiation, wind velocity, cooling load, and other loads, were all formulated using uniform distributions. The uncertainty of the cooling load can be traced to many aspects, such as the stochastic nature of weather data, the lack of building information, or the unpredictability of internal disturbances. Therefore, in [2,11–12], many empirical probability distribution functions, like the normal distribution for indoor temperatures or the triangular distribution for occupancy, were presented and lots of MC simulations were performed to account for the uncertain cooling load in the cooling source design. The results indicated that the peak cooling load and annual average cooling load in each MC simulation vary greatly given the uncertainties. As a result, a question arises: how many iterations of the MC simulation are sufficient to describe load uncertainty comprehensively? Generally, the researchers assumed that 1000 simulations are sufficient to generate the uncertain load, but no quantized indexes were presented to evaluate whether the results converge or not [13–14]. In [3], the minimum MC simulation number was studied by discussing the peak load and load distribution, and 391

Applied Energy 241 (2019) 390–403

J. Niu, et al.

1.2. Mathematical models for cooling source design

operational variables in the second stage [7]. However, Mavromatidis et al. pointed out that a large number of iterations in the genetic algorithm and an adequate number of samples with the MC method must result in computational issues; otherwise the genetic algorithm would not converge and the solution of model would not be accurate enough [29]. Therefore, there is a trade-off between the solution’s accuracy and the computational time. In [29], a feature-based clustering approach was applied to reduce the time series load scenarios that consider multiple uncertainty parameters, and then a two-stage stochastic program was performed efficiently for an optimal district energy system design. Obviously, there is a compromise in the accuracy of the model’s solution because some of the uncertain load scenarios are overlooked. Reliability-based design optimization (RBDO) is usually used for engineering designs when uncertainty is being considered. Li et al. proposed a single-loop deterministic method for RBDO that converts the probabilistic constraints to approximate deterministic constraints [30]. For optimizing the dispatch of a combined cooling, heating, and power (CCHP) system while considering demand uncertainty, Hu et al. developed a stochastic optimization model in which the probability constraints were added and converted into deterministic constraints according to [30] in order to guarantee that the optimized CCHP strategy was reliable [31]. However, the results in [30] indicated that the accuracy could be drastically reduced when the computational efficiency was improved. Even though the probabilistic method has been widely applied to deal with the optimization of an energy system while considering uncertainty, the different approaches discussed in the above paragraphs also indicate that some shortcomings still exit. Given the probability distribution function, an empirical function may lead to a deviation from the real situation; yet it is difficult to collect huge samples of these uncertainty parameters to determine of the probability distribution function. Moreover, it is impossible to consider all the uncertain scenarios in the optimization model because it is prohibitive in terms of computational expense. Therefore, some other models or methods have been developed to treat the uncertainty in the optimization of an energy system. For instance, Li et al. proposed a dual-interval mixed-integer linear model in which the uncertain parameters are expressed as intervals and dual intervals [32]. In [33], the worst-case realization of uncertainties associated with energy demands and renewable energy resources are considered for optimally sizing and siting tri-generation equipment. However, the worst-case scenario is the energy demand taking the upper bound of uncertainty while the renewable generations taking the lower bound of uncertainty. Undoubtedly, the solution obtained in [33] may be overly conservative. Hence, Akbari et al. introduced a new model for the design of a district energy system [34]. In the proposed model in [34], the uncertainty parameter is represented by the mean value; the constant perturbation and the conservatism level could be determined according to the decision-maker’s preference. Similarly, Jiao et al. proposed a two-stage robust model for scheduling a building energy system wherein three parameters, called the “budget of uncertainty” are used to adjust the level of conservatism in the optimal solution [35]. Even though the interval of an uncertainty parameter is relatively easy to determine compared to its probability distribution function, the solution obtained by the interval optimization method is still overly conservative. Wang et al. pointed out that the dependency problem and wrapping effect would make the bounds of the solution much wider than the true solution region [36]. Therefore, the affine arithmetic algorithm was introduced to represent the interval because it had been proven to provide more compact bounds [36].

Generally, mathematical optimization models are developed to assist the design of an energy system [4–6]. Binary and Special-Order-Set (SOS) variables were introduced in [4] and [6] to address the nonlinear problems and to further formulate a mixed integer linear programming (MILP) model for the design of a distributed energy system [4]. Guo et al. presented a two-stage optimal planning and design method for a combined cooling, heat, and power microgrid system, in which the non-dominated genetic algorithm was applied to solve the optimal design problem in the first stage and MILP model was used to optimize the dispatch problem in the second stage [5]. However, most of the models are formulated as deterministic models in which all the parameters are deemed to be known. As Walker said [21] ‘any deviation from the unachievable ideal of completely deterministic knowledge of the relevant system is uncertainty.’ Hence, the energy system based on a deterministic model is no longer reliable or accurate as lots of uncertainties are overlooked. So, another critical problem is how to combine the uncertain parameters and mathematical optimization model together to get a robust energy system. Sensitivity analysis, the probability optimization method, and the interval optimization method are commonly used to investigate uncertainty in the optimal design of an energy system. In sensitivity analysis, the influence of uncertain parameters on the objective value is examined by varying each uncertain parameter range and then ranking its degree. In [22], a multi-stage sensitivity analysis approach was proposed to identify the main uncertainty parameters for designing a zero/low building. Li et al. adopted average, uncertainty and peaks of energy demands to describe the complex energy demands; they performed a sensitivity analysis to investigate the influence of the energy demands on optimal facility scheme [23]. In order to alleviate the drawback of a local sensitivity analysis, [24] performed a one-way sensitivity analysis, a two-way sensitivity analysis, and a multiway sensitivity analysis to investigate the impacts of uncertain input parameters on the building performance. Similarly, Mavromatidis et al. proposed a two-step global sensitivity analysis combining the Morris and variance-based Sobol methods to examine the most influential parameters driving the total cost of the district energy system [25]. Overall, despite the fact that the sensitivity analysis approach has solved the problem of ‘which of the uncertain input factors is more important in determining the uncertainty in the output of interest?’ [26], it still cannot obtain a single global optimization energy system configuration. Consequently, a model-based probabilistic method was developed for decision-making under uncertainty [10,15–16,27–28]. For instance, Lu et al. assumed that the uncertainty of solar radiation, wind velocity, cooling load, and other loads all followed a uniform distribution. Based on the above assumption, an exhaustive search and MC simulations were performed to a design renewable energy system [10,27]. In [15], the probability distribution of the actual pressure head is introduced to estimate the pressure loss. So, the water pump system could be designed to take uncertainties into account [15]. Similarly, in [16] and [28], some input factors, i.e., weather conditions, occupants, and heat transfer coefficients, were described by probability distributions of different types. Then, a competitive analysis was performed for some typical configurations of chiller plants to achieve an optimal design. The probability distribution function accurately describes the characteristics of uncertain parameters, but an exhaustive search method and typical scenario analysis method are difficult to obtain for the global optimal configuration of the energy system. Therefore, a heuristic algorithm, such as the genetic algorithm or a numerical algorithm, such as the simplex method should be implemented to solve the energy system optimization model with uncertain parameters. Zhou et al. developed a two-stage stochastic programming model for the optimal design of a district energy system, in which the genetic algorithm was performed to search the design variables in the first stage and a MILP model combined with the MC method was used to deal with the

1.3. Literature summary and main objective of this paper Many studies have thus conducted an optimization of the energy system, but some limitations should be addressed and are summarized as following: (1) The main purpose of performing a sensitivity analysis is to assess the relative importance of a model’s uncertainty input 392

Applied Energy 241 (2019) 390–403

J. Niu, et al.

factors. As a result, a global optimization solution cannot be determined by this method. (2) The characteristics of uncertainty parameters can be described exactly by the probability distribution function combined with a large number of MC simulations, which cause prohibitive computational expenses. Moreover, heuristic algorithms used in the first stage of the two-stage programming model may only achieve a local optimum. (3) Even though the interval of an uncertainty parameter is easily obtained, the solution using the interval optimization method usually deviates from the true solution region. Consequently, this paper develops a new robust model for the optimal design of a building cooling source considering its load uncertainty. The proposed model has some of the following features: i. The computational expense is dramatically reduced along with the time series cooling load reduction, which is completed using a new load reduction method introduced in this paper. ii. The proposed model can be easily solved by a numerical algorithm implemented in some optimal solver like GAMS/Gurobi [37] and thus, a global optimization solution can be obtained. iii. the computational efficiency is not limited by the number of MC simulations. Therefore, numerous MC simulations can be performed to achieve an exact presentation of the uncertainty parameters. The paper is structured as follows: Section 2 presents the formulation of a robust model based on load reduction for the optimal design of a building’s cooling source. In Section 3, a case study is introduced to investigate the effectiveness of the proposed model. Section 4 discusses the solution’s accuracy by comparing it to the deterministic model (presented in Appendix A). Section 5 concludes.

implemented to reduce the computing challenge, which is described in detail in Section 2.2. 2.1. Load scenario reduction based on the cooling load bin (CLB) method In this paper, MC sampling is used to generate multiple scenarios for the uncertain parameters which are described by probability distribution functions; then EnergyPlus [38] is launched to generate the cooling load. Generally, the researchers assume that 780 [3], 1000 [11], or 1500 [39] simulation trials are sufficient to guarantee convergence of the MC sampling. Obviously, the more MC samplings, the better accuracy of the simulation results. However, if all the cooling load scenarios are applied in an optimization problem, the deterministic model will be difficult or even impossible to be dealt with directly in the optimal solvers (as illustrated in Appendix A). Therefore, this study proposes a cooling load bin (CLB) method for the load scenarios reduction that is based on the following hypothesis: when the time series cooling loads at different points in time are equal or approximately equal, the cooling source should launch the same cooling strategy. The idea of CLB method can be illustrated by considering the following raw time series cooling data shown in Fig. 2(a). To apply the CLB method, the first step is to bin the range of time series cooling loads—that is, divide the entire range of cooling loads into a series of intervals, defined as the CLBs, as shown by the x-axis in Fig. 2(b) and Fig. 2(c). Then count how many cooling loads fall into each bin using Eq. (1), to determine the frequency of the cooling loads in each bin, as shown in Fig. 2(b).

fi =

2. Methodology

c

(1)

where fi refers to the frequency of the cooling loads in bin i ; c equals 1 if cload is grouped in CLB i ; otherwise, it equals 0; and cload is the cooling load at time . ¯ i , of the mean of the cooling loads Finally, an estimate, cload grouped in each bin i is calculated by Eq. (2) and shown in Fig. 2(c).

The process and steps of the robust optimal design method are illustrated in Fig. 1. In the first step, the uncertainties of the input variables are quantified by probability distributions. MC sampling is then performed to generate a combination of all the variables and write them into EnergyPlus files to get a time series cooling load. A basic hypothesis is that the uncertainties and their combined effect can be propagated to the cooling load equivalently and completely as long as enough simulations are performed. Details of the stochastic simulation are available in [2]. In the second step, the cooling load in the time domain is converted to that in frequency domain. Detailed technology is presented in Section 2.1. In the third step, the cooling source design model is proposed, in which the cooling load in the frequency domain is

¯ i = cload

cload (coolingloadgroupedinbini ) fi

(2)

¯ i is the mean value of the cooling loads grouped in bin i . where cload Thus, the time series cooling load is replaced by the cooling load frequency and the mean of the cooling loads. When considering the cooling load as reduced by the CLB method, the optimal model for the task of optimal cooling source design, shown in Appendix. A is

Fig. 1. Process of the robust optimal design method of the cooling source. 393

Applied Energy 241 (2019) 390–403

J. Niu, et al.

Fig. 2. (a). Time series cooling load in one MC simulation; (b) The frequency distribution of CLBs; (c) The mean value of grouped cooling loads.

reformulated and described in Section 2.2 in detail.

CRF T =

2.2. Robust optimization model for cooling source design

I

fi ·EleiT, j · Ep·

T

Y

(6)

• Thermal balance

The system investment cost over the project lifetime is converted to an annual value by the capital recovery factor (CRF). As shown in Eq. (4), the system investment cost includes the chiller investment and auxiliary equipment investment, which are calculated by a correction factor ( = 1.2 ). Eq. (5) describes the CRF, in which r is the discount rate and L is the lifetime of the equipment. Eq. (6) describes the operational cost expenditure, expressed as follows:

T = A, B, C

(5)

Commonly, the total energy consumption is obtained by accumulating the hourly energy consumption, while in this proposed model, it is calculated by accumulating the product of the cooling load frequency f and the average energy consumption Ele at each CLB i . Y denotes the times of the MC simulations.

(3)

minimize Tcost = In + Op

LT

NT

T = A, B , C i = 1 j = 1

In this paper, the equivalent annual cost (Tcost ) is used as the objective function, consisting of the investment expenditure (In ) and the operational cost expenditure (Op ). The most suitable configuration of the cooling source is obtained by minimizing the total annual cost, as follows:

N T · CaT ·UpT ·CRF T ·

(1 + r )

Op =

• Objective function

In =

r 1

¯ i ) at each bin i is satisfied by the The mean of the cooling load (cload chillers as shown in Eq. (7). The thermal balance constraint (Eq. [8]) is robust, as all the time series cooling load scenarios could be considered in this model. NT T = A, B , C j = 1

qiT, j

¯ i cload

• Energy conversion and other constraints

T

(4) 394

(7)

Applied Energy 241 (2019) 390–403

J. Niu, et al.

For the convenient control and maintenance in late operation, no more than 5 of each type of chiller could be configured in this study:

NT

3. Case study 3.1. Description

(8)

5T = A, B, C

As a case study for this paper, the optimal design of a cooling source for a hospital building in Tianjin, China, is investigated. This hospital has fifteen floors with around 10800 m2 area and is conditioned with Make-up Air Units (MAUs) from June 1st to September 15th, a total of 2568 h (H = 2568). The building outline is shown in Fig. 4(a). The fifteen-story building mainly includes two functional departments; one is the outpatient department, which is around 4320 m2, and the other is the inpatient department, which is around 6480 m2. The inpatient department is occupied all the time and the outpatient department is mainly occupied from 07:00 to 19:00. Other operational schemes in the building planning stage (equipment scheme, light scheme, etc.) are listed in detail in [41]. The overall structural representation of the cooling source for the hospital is illustrated in Fig. 4(b). The task of the optimal cooling source design is to determine what type and how many of these chillers should be installed in order to minimize the total cost under uncertain cooling load scenarios. The cooling water pump, chilled water pump, chiller, and cooling tower are pre-designed on a one-to-one matching basis that is used in [15]. Considering the chillers could be customized by suppliers, the chillers are regarded as discrete equipment, namely the rated capacity is known. In this paper, there are three types of chillers available for optimal selection and their information is shown in Table 1. The optimal design model described in Section 2.2 is automatically encoded in GAMS using Python script. In this study, co-simulation is easy to implement in Python by linking GAMS and EnergyPlus. In this paper, all the work is done on an Intel(R) Core (TM) i5-6200U CPU of 2.40 GHz with a memory size of 8.00G.

Binary variables (bin ) are introduced to describe the status of chillers; bin = 1 means the chiller is on and vice versa. The number of running chiller is constrained by Eq. (9). 5

bin iT, j

NT

(9)

j=1

As for the electricity consumption, it is mainly determined by the cooling load distribution and the energy efficiency. Generally, the coefficient of performance (COP) of an electric chiller varies mainly depending on the partial cooling load [2,40], as shown in Fig. 3(a). The chiller power (Ele ) can also be calculated by the cooling load (q ) directly. In this paper, the electricity consumption of a chiller is estimated by a cubic function equivalently, as shown in Eq. (10) and Fig. 3(b), which simplifies the linearization process as expressed in Eqs. (15)–(18). A binary variable (bin ) is introduced in Eq. (10) to ensure that there will be no electricity consumption when the chiller is off.

EleiT, j = aT ·(qiT, j )3 + bT ·(qiT, j ) 2 + cT · qiT, j + dT ·biniT, j

(10)

The cooling load of each type of chiller is constrained within its lower and upper bounds of capacity, which can ensure that the chiller is running within a reasonable range.

0

qiT, j

CaT

(11)

As described in [40], if a chiller is on, the outputs of the same type of chiller, which are installed in parallel in this paper, should be equal. In this paper, that constraint is defined by introducing an auxiliary variable qiT , as follows:

bin iT, j ·M

qiT

(1

qiT

CaT

(12)

biniT, j )· CaT

qiT, j

qiT

5.00

(13) (14)

When the frequency of the cooling load at some CLBs is equal to 0, the chillers should be off, which is achieved by adding a constrain, as expressed in Eq. (15). This operation reduces the search space and speeds up the calculation.

bin iT, j

(1

seqiT, s = 1

s=0

Ele iT, j

biniT, j ·M

biniT, j )· CaT

qiT, j

Ele

seqiT, s ·qsT T i, j

seqiT, s ·ElesT

3.00

20

30

40

50

60

70

80

90

100

350.0

400.0

part load ratio (%) 90.0 80.0

chiller power (kW)

seqiT, s · ElesT S

biniT, j )· CaT

(1

y = -0.0004x2 + 0.0646x + 2.0705 R² = 0.98

3.50

2.00 10

In this paper, the SOS of Type 2 (SOS2) variables are introduced to linearize Eq. (10). SOS2 variables are a set of variables, such that at most two variables within the set may have nonzero values and these variables have to be adjacent [2,40]. The detailed linearization process is presented as follows:

seqiT, s · qsT

4.00

2.50

(15)

fi

Chiller (Type C)

4.50

chiller COP

qiT, j

(16) (17)

(18)

70.0 60.0

Chiller (Type C) y = 2E-06x3 - 0.001x2 + 0.3348x + 2.4904 R² = 0.99

50.0 40.0 30.0 20.0 10.0 0.0 0.0

(19)

50.0

100.0

150.0

200.0

250.0

300.0

cooling power (kW)

Where seq is the SOS2 variable, qsT and ElesT are these cooling load and chiller power of set S . M refers to a large number as an auxiliary constant associated with the variable bin iT, j .

Fig. 3. (a). Performance curve of chillers at different load ratios, (b) Chiller power at different cooling loads. 395

Applied Energy 241 (2019) 390–403

J. Niu, et al.

Fig. 4. (a) Illustration of the hospital; (b) Set of candidate chillers and schematic diagram of a cooling source.

from 1981 to 2010 (30 years) have been collected and used to represent the uncertainties. Other uncertainty parameters are defined by different probability distribution functions and the detailed uncertainty parameters used as inputs in the cooling load uncertainty study are listed in Table 2 [3,41–42]. Uniform distributions are adopted in Table 2 when there are no accurate distributions [12]. The variable’s mean value, variance or range are obtained from the local design guidelines and assumptions are based on the designer’s experience. To generate uncertain cooling load scenarios, EnergyPlus [3] is used as the building performance simulation tool and the probability distributions shown in Table 2 are assigned to its uncertain input parameters. These uncertain input parameters are drawn with MC simulations and are automatically encoded into the EnergyPlus file using the Python script. Finally, 1000 cooling load profiles (Y = 1000) are generated by launching EnergyPlus via Python script.

Table 1 Information about the chillers. Type

A B C

Capacity (kW)

Unit price (/kW)

Rated COP

1055 528 352

939 945 1002

4.82 4.48 4.24

Coefficient a

b

c

d

2.00E−07 4.00E−07 2.00E−06

−0.0003 −0.0003 −0.001

0.3226 0.2695 0.3348

4.4007 6.9761 2.4904

3.2. Generating the cooling load with uncertainties Many uncertainty parameters that have a significant impact on the cooling load are considered in the cooling load calculation. The weather data is a time series and the dry temperature, solar radiation, and other meteorological parameters are not independent. Hence, weather data Table 2 Information on the inputs used in the cooling load uncertainty study. Items

Parameters

Weather data

Reference value

Distribution

—

1981–2010

Building construction

Wall conductivity W/m2 K Window conductivity W/m2 K Roof conductivity W/m2 K

0.45 2.85 0.25

T (0.36, 0.60, 0.45) T (2.50, 3.24, 2.85) T (0.20, 0.30, 0.25)

Indoor condition

Occupancy density m2/person Power density W/m2 Lighting density W/m2 Ventilation rate 1/hour Indoor set temperature ℃

6 9.6 6.4 0.7 24

T (3, 9, 6) T (3.6, 14.4, 9.6) T (2.4, 9.6, 6.4) U (0.4, 0.9) N (24, 0.1)

Note: T denotes the triangle distribution, U denotes the uniform distribution, N denotes the normal distribution. 396

Applied Energy 241 (2019) 390–403

J. Niu, et al.

3.2.1. Uncertainty analysis of the cooling load After performing 1000 MC simulations, the time series cooling loads are obtained, as shown in Fig. 5(a). It can be seen that these cooling loads vary diurnally mainly owing to the schedule of the indoor temperature. A cooling load from one day is selected and shown in Fig. 5(a) (grey background) in which significant variations of the quantity of the cooling loads are observed, indicating a high dependence of the building cooling load on the realization of the uncertain input parameters. The design cooling load usually corresponds to an unmet hour [43]. In this study, the unmet hour is 50 and the corresponding load at that point in time is defined as the design cooling load. The generated design cooling loads of each MC simulation are put into a histogram, as shown in Fig. 5(b). The design cooling load distribution is located between 704 kW and 1707 kW. The average is about 1207 kW and the standard variance is 138 kW. Besides the design load, Fig. 6(a) illustrates the variation in the cumulative hourly cooling load. A violin plot in Fig. 6(b) is depicted for a better understanding of the uncertainty of the annual cooling energy demand. Overall, the uncertain input parameters result in significant variation of the cooling load for both the design load and the energy demand. Therefore, all uncertain cooling loads (1000 MC simulations) should be considered comprehensively in the cooling source design.

3.3. Results

3.2.2. Cooling load reduction scenario It is a well-known disadvantage of a deterministic programming model that is, scales in size with the number of cooling load scenarios considered. In this work, the cooling load scenario reduction is performed first based on the method in Section 2.1, and then the proposed robust model is applied to the design of cooling source. Considering all the cooling load scenarios, the design cooling load is first determined as 1259 kW (average unmet 50 h in each simulation) and then all the hourly cooling loads are limited to the design cooling load. Fig. 7(a) shows the frequency distribution and Fig. 7(b) shows the mean value of the cooling load in each CLB for the 1000 cooling load scenarios. The red dotted line presents the mean of the cooling load without being limited by the design value and the bar chart is the mean of the cooling load limited by the design value that will be used to design the cooling source.

3.3.1. Effectiveness verification of the robust optimization model The cooling source is optimally configured by the robust design model (C0 shown in Table 3); the eight cooling source configurations (C1–C8 shown in Table 3) are obtained based on the traditional design method outlined in the above rules. The results are shown together in Table 3 and Fig. 8 for comparation. First, it is interesting that the cooling source configuration C5 is the same as the configuration obtained by the robust optimization design model, while the annual total cost (Tcost) of the former is slightly less. That is because different cooling loads are used to calculate the operational cost: the time series cooling load as shown in Fig. 5(a) is used in C5, while the reduced cooling load as shown in Fig. 7 is used in C0. This results also demonstrates that the robust model could obtain the mean optimization from 1000 cooling load scenarios because only a small difference exists between C0 and C5 (−0.004%). Moreover, it can be observed that the minimum objective value is obtained by the

In this section, the results of the robust optimization model are presented and analyzed. To demonstrate the validation of the robust model, the results obtained by the robust model are compared with those obtained using the traditional design method, which follows three rules: Rule 1: The number of chillers is determined only from the design cooling load, as shown in Fig. 7(a) (1259 kW). Rule 2: The total capacity of all configured chillers should be just more than or equal to the design cooling load. For instance, the total capacity of three chiller Bs is 1584 kW, just more than 1259 kW, while it reaches to 2112 kW when four chiller Bs are configured, far more than 1259 kW. Rule 3: When the cooling source configuration is determined based on Rule 2, the operational cost is optimally calculated by the deterministic model (Appendix A) for each cooling load scenario, in which the number of chillers is constant. Finally, the average operational cost for the 1000 cooling load scenarios is used to represent the annual operational cost of the corresponding cooling source.

Fig. 5. (a) Time series cooling load from 1000 simulations; (b) The frequency distribution of the design cooling load of each simulation. 397

Applied Energy 241 (2019) 390–403

J. Niu, et al.

Fig. 6. (a) The cumulative hourly cooling load from 1000 simulations; (b) Violin plot for the annual cooling energy consumption from 1000 simulations.

proposed robust model compared with C1–C4 and C6–C8. Therefore, the proposed model is effective in optimizing the cooling source design with uncertain cooling loads. Regarding the operational cost (Op), Fig. 8 presents the variation of the operational cost obtained by the deterministic model (Appendix A) using violin plots and the operational cost obtained by the robust optimization design model using the dotted line. The orange dots represent the mean value of the operational cost of each configuration (C1–C8). Even though the mean operational cost of C8 is less than that of C0 (2.0%), the investment of the former is far more. The other op-

eration costs are all more than (C1–C4, C6, C7) or nearly equal to (C5) that of C0. These results reveal that the robust model can obtain the optimal solution when considering the investment expenditure (Tcost) and the operational cost (Op) comprehensively, but in the traditional design method, the configuration is first pre-designed and then the operation is optimized by the deterministic model. In addition, the robust design model can guarantee the global optimization because it treats the large number of uncertain scenarios synchronously, while the traditional design method can only obtain get some feasible configurations, which are also known as local optimal solutions.

Fig. 7. (a) The frequency of the cooling load in each bin; (b) The mean of the cooling load in each bin.

398

Applied Energy 241 (2019) 390–403

J. Niu, et al.

Table 3 Results comparison between the robust model and the traditional design method. Configuration

The number of chillers Chiller A

Chiller B

In (yuan)

Op (yuan)

Tcost (yuan)

Difference of Tcost (%)

Chiller C

Robust design

C0

1

0

1

151529.8

341357.3

492887.1

—

Traditional design method

C1 C2 C3 C4 C5 C6 C7 C8

2 0 0 1 1 0 0 1

0 3 0 1 0 1 2 1

0 0 4 0 1 3 1 1

223489.5 168848.1 159140.0 168027.4 151529.8 175637.7 152359.4 207812.5

363765.4 351343.7 347558.1 348028.2 341338.7 348243.3 347995.3 334763.2

587254.9 520191.8 506698.1 516055.6 492868.5 523881.0 500354.7 542575.7

19.146 5.540 2.802 4.701 −0.004 6.288 1.515 10.081

3.3.2. Test of efficiency of the robust optimization model The convergence rate is a crucial index for testing the efficiency of the robust model. The relative gap between the objective bound and objective value is a crucial factor in Gurobi (an optimization solver), which is used to measure whether the optimal solution is obtained or not. Fig. 9 shows the objective bound, the objective value, and the relative gap of each iteration in Gurobi. The relative gap drops gradually and finally converges to 0.001%, while the objective value (Tcost) also converges to 492887.1, which is considered to be an optimal solution. Another important observation is that the solution time is 3.10 s, illustrating that the robust optimization design model could be solved effectively even with numerous cooling load scenarios.

in Section 3, the accuracy of the proposed robust model should be further discussed in this section. The width of the CLB is the key factor in determining the accuracy of the robust model. If the width becomes infinitesimal, the robust model developed in Section 3 becomes identical to the deterministic model formulated in Appendix A, since the number of CLBs will be close to the total hours in Y times of MC simulations (Y equals 1000 and the total hours is 2,568,000 in this paper). Thus, a computational tractability issue arises in the robust optimization model because its size is dramatically increased. CLBs with a large width can make simplify the calculations, but results from the cooling load distribution deviate from reality immoderately. Therefore, two cases are illustrated to analyze the impact of the width of the CLBs on the accuracy of the robust optimization model.

4. Further discussion of the robust model

Case one: Accuracy verification for one cooling load scenario; Case two: Accuracy verification after many cooling load scenarios.

Besides the effectiveness of the robust optimization model verified

Fig. 8. Operational cost (Op) calculated by the robust model and the deterministic model (Appendix A). 399

Applied Energy 241 (2019) 390–403

J. Niu, et al.

540000

14.0%

520000

12.0%

objective

8.0%

480000

6.0%

460000

4.0%

440000 420000

gap

10.0%

500000

2.0% 0.00 0.20 0.40 0.80 1.00 1.20 1.40 2.60 1.80 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10

0.0%

solution time (s)

Objective value

Objective bound

Gap

Fig. 9. The convergence process of the robust optimization design model.

4.1. Accuracy verification for one cooling load scenario

robust optimization model exponentially decreases as the width of the CLBs increases. This result can be primarily attributed to the fact that the number of CLBs also decreases exponentially as the width of the CLBs increases, as shown in Fig. 11.

In case one, one cooling load profile is selected from 1000 MC simulations as the boundary condition for the robust optimization model and the deterministic model. In the robust optimization model, the time series cooling load is reduced by the method in Section 2.1 with 1 kW, 2 kW, 3 kW, … 100 kW as the widths of the CLBs, respectively. In the deterministic model, the time series cooling load is employed directly because the size of the model is computational. The results obtained by the two models separately are compared quantitively by the absolute relative deviation (ARD) of their objectives ( ARDkTcost ), defined as:

ARDkTcost =

|Tcost D TcostkR | k = 1 kW, 2 kW, 3 kW 100 kW Tcost D

4.2. Accuracy verification for many cooling load scenarios To further investigate the accuracy of the robust optimization model, 1, 10, 50, etc., of cooling load scenarios are firstly reduced with 1 kW, 2 kW, 3 kW, …, 100 kW as the widths of the CLBs and then applied to the robust optimization model, respectively. Since using many cooling load scenarios in the deterministic model is not possible due to its computational limitations, only the robust optimization model is performed in this section. Here, an index named the ‘moving absolute relative deviation’ (MARD) of the objective is obtained from the robust optimization model and is developed to evaluate the accuracy quantitively. The MARD is defined as:

(20)

where is the objective in the deterministic model and TcostkR is the objective in the robust optimization model with the CLB k kW. As can be seen in Fig. 10, the ARDkTcost presents the trend of divergence along with the increasing width of the CLBs. When the width is Tcost less than 10 kW, the ARD10 is less than 0.5%, indicating that the robust optimization model is as accurate as the deterministic model. A second observation is that the solution time of the robust optimization model as shown in Fig. 11, across all cases (k = 1 kW, 2 kW, 3 kW, …, 100 kW), is always less compared to that resulting from the deterministic model, which required 568 s. Moreover, the solution time of the

Tcost D

MARDkTcost =

TcostkR+ 1 TcostkR k = 1 kW, 2 kW, 3 kW, TcostkR+ 1

(21) In total, 1300 calculations are performed by the robust optimization model and the results are summarized in Fig. 12 for the convenience of comparison. Across all the cases considered (1 simulation, 10 simulations, …), it can be observed that the MARDkTcost converges to about

10.0%

8.0%

ARD k

6.0%

4.0%

2.0%

0.0%

1

10

19

28

37

, 100 kW

46

55

64

73

82

91

100

width of cooling load bin (kW) Fig. 10. Absolute relative deviation of the objectives in the robust optimization model and the deterministic model. 400

solution time (s)

35.0

1600

30.0

1400 1200

25.0

1000

20.0

800

15.0

600

10.0

400

5.0

0.0

200 1

10

19

28

37

46

55

64

73

82

91

100

number of cooling load bins

Applied Energy 241 (2019) 390–403

J. Niu, et al.

0

width cooling load bin (kW) Solution time (s)

Number of cooling load bins

Fig. 11. Solution time and number of CLBs for different bin width.

0.7% when the width of CLBs decreases to 10 kW. Therefore, as long as the width of the CLBs is less than 10 kW, a stationary solution could be obtained using the proposed robust optimization model. The solution time is also summarized and drawn in Fig. 12. It can be seen that the efficiency of the robust optimization model is not dependent on the number of cooling load scenarios. Therefore, with the robust optimization model proposed in this paper, a large number of cooling load scenarios could be considered to design an optimal robust cooling source within a reasonable computational effort.

of the proposed model. Conclusions are summarized as follows: 1) The CLB method is developed and proven to reduce the time series cooling load scenarios effectively. Furthermore, the size of the robust optimal model for the design of building cooling sources is reduced correspondingly and can be solved by Gurobi. 2) An optimal cooling source is configured using the proposed robust model, which considers lots of cooling load scenarios. The effectiveness of the proposed robust optimal model is further discussed and verified by comparing the cooling source configuration with one obtained by using the traditional method. 3) In terms of the accuracy and efficiency of the proposed model, the results showed that when the width of the CLBs is less than 10 kW, the absolute relative deviation of the objective is less than 0.5% compared to the deterministic model that uses one cooling load scenario; the moving absolute relative deviation of the objective can also converge under various cooling load scenarios. Furthermore, it is interesting that the efficiency of the robust optimization model does not depend on the number of cooling load scenarios. Therefore,

5. Conclusions This paper proposed a robust optimization model for the design of building cooling sources under a large number of cooling load scenarios. A CLB method is proposed to reduce the numerous uncertain cooling load scenarios. Therefore, this proposed model balances the consideration of a large number of cooling load scenarios and solving the model with a reasonable computational effort. A case study is given as an example to demonstrate the effectiveness, efficiency and accuracy 6.0%

0.6

5.0%

30 25

0.4 0.3

3.0%

9

10

11

12

13

14

15

20

0.2

15

2.0%

10

1.0% 0.0%

solution time (s)

0.5

4.0%

MARDk Tcost

35

0.7

5

1

8

15

22

29

36

43

50

57

64

71

78

85

92

99

0

width of cooling load bin (kW) 1 simulation 200 simulations 600 simulations 1000 simulations

10 simulations 300 simulations 700 simulations solution time (s)

50 simulations 400 simulations 800 simulations

100 simulations 500 simulations 900 simulations

Fig. 12. The MARD of the objective and the solution time based on different numbers of cooling load simulations and different widths of the CLBs.

401

Applied Energy 241 (2019) 390–403

J. Niu, et al.

the uncertainty parameters can be estimated accurately by performing lots of MC simulations.

should also be introduced. Acknowledgements

Regarding future work, a possible model extension might be to introduce uncertain electrical demand, uncertain heating demand and their correlation. Additionally, the renewable energy uncertainty

This work is supported by the National Key Research and Development Program of China (Project No. 2017YFB0903404).

Appendix A. Deterministic mathematical model for the optimal cooling source design In this appendix, the formulation of the deterministic cooling source design is provided.

• Objective function The objective and investment expenditure in the deterministic model are the same as those in the robust model, as shown in Eq. (1) and Eq. (2). The annual average operating cost is calculated with the uncertain cooling load, which is presented using Y times of MC simulations with H hours in each simulation, expressed as follows: NT

Y

H

EleTj, y, ·Ep ·

Op =

T

Y

(A.1)

T = A, B, C j = 1 y = 1 = 1

• Thermal balance The output and demand of the cooling (Eq. (A.2)) should be balanced at any point in time. It is expressed as: NT

Y

H

T = A, B , C j = 1 y = 1 = 1

qTj, y,

cload y,

(A.2)

• Energy conversion and other constraints Additional constraints are added to impose operational constraints on the chiller considered. Eq. (A.3) is introduced to limit the maximum number of configured chillers.

NT

(A.3)

5T = A, B, C The number of operating chillers should be less than N T , as follows:

5

binTj, y,

NT

(A.4)

j=1

For the precise calculation of electricity consumption, a cubic function is formulated, as shown in Eq. (A.5). A binary variable is introduced to ensure that there is no electricity consumption when the chiller is off.

EleTj, y, = aT ·(qTj, y, )3 + bT ·(qTj, y, )2 + cT ·qTj, y, + dT · binTj, y,

(A.5)

The output of each chiller is constrained within its lower and upper bounds of capacity at any point in time, as follows:

0

qTj, y,

CaT

(A.6)

The outputs of same type of chillers are constrained to be equal via introducing a variable

qTj, y,

binTj, y,

q yT,

(1

q yT,

CaT

q yT,

and a large constant (big-M) M.

·M

(A.7)

binTj, y, )· CaT

qTj, y,

q yT,

(A.8) (A.9)

In the deterministic model, the SOS2 variables are also introduced to linearize Eq. (A.5), as follows:

seqsT, y, ·qsT

(1

seqsT, y, ·ElesT S

(1

binTj, y, )· CaT binTj, y, )· CaT

qiT, j

seqsT, y, · qsT

EleTj, y,

(A.10)

seqsT, y, · ElesT

(A.11)

seqsT, y, = 1

(A.12)

s=0

×Y×H+T×Y×H+T×S×Y×H+6 In Eqs. (A.1)–(A.12), the bolded italic symbols are variables. The deterministic model includes 3 × T × variables, which closely relates to the number of MC simulation. The MILP in this paper is exponential time solvable at least. Namely, the time complexity of the MILP in this paper is O(2n). Therefore, finding the right balance between considering a large number of MC simulations and solving the deterministic model with a reasonable computational effort might be difficult or even impossible. NT

402

Applied Energy 241 (2019) 390–403

J. Niu, et al.

References [22]

[1] Handbook A. ASHRAE Handbook–Fundamentals. Atlanta, GA; 2009. [2] Gang W, Wang S, Augenbroe G, Xiao F. Robust optimal design of district cooling systems and the impacts of uncertainty and reliability. Energy Build 2016;122:11–22. [3] Cheng Q, Wang S, Yan C, Xiao F. Probabilistic approach for uncertainty-based optimal design of chiller plants in buildings ☆. Appl Energy 2015. [4] Milan C, Stadler M, Cardoso G, Mashayekh S. Modeling of non-linear CHP efficiency curves in distributed energy systems. Appl Energy 2015;148:334–47. [5] Guo L, Liu W, Cai J, Hong B, Wang C. A two-stage optimal planning and design method for combined cooling, heat and power microgrid system. Energy Convers Manage 2013;74(10):433–45. [6] Marnay C, Stadler M, Siddiqui A, DeForest N, Donadee J, Bhattacharya P, et al. Applications of optimal building energy system selection and operation. Proc Instit Mech Eng Part A J Power Energy 2013;227(1):82–93. [7] Zhou Z, Zhang J, Liu P, Li Z, Georgiadis M, Pistikopoilos E. A two-stage stochastic programming model for the optimal design of distributed energy systems. Appl Energy 2013;103(1):135–44. [8] Mavromatidis G, Orehounig K, Carmeliet J. A review of uncertainty characterisation approaches for the optimal design of distributed energy systems. Renew Sustain Energy Rev 2018;88:258–77. [9] Mahadevan S. Uncertainty quantification for decision-making in engineered systems. Proceedings of the International Symposium on Engineering under Uncertainty: Safety Assessment and Management (ISEUSAM – 2012). 2013. p. 97–117. [10] Lu Y, Wang S, Yan C, Huang Z. Robust optimal design of renewable energy system in nearly/net zero energy buildings under uncertainties. Appl Energy 2017;187(2):62–71. [11] Gang W, Wang S, Xiao F, Gao D. Robust optimal design of building cooling systems considering cooling load uncertainty and equipment reliability. Appl Energy 2015:265–75. [12] Gang W, Augenbroe G, Wang S, Fan C, Xiao F. An uncertainty-based design optimization method for district cooling systems. Energy 2016:516–27. [13] Amusat OO, Shearing PR, Fraga ES. Optimal integrated energy systems design incorporating variable renewable energy sources. Comput Chem Eng 2016:21–37. [14] Gang W, Wang S, Shan K, Gao D. Impacts of cooling load calculation uncertainties on the design optimization of building cooling systems. Energy Build 2015;94:1–9. [15] Li H, Wang S. Probabilistic optimal design concerning uncertainties and on-site adaptive commissioning of air-conditioning water pump systems in buildings. Appl Energy 2017;202. [16] Cheng Q, Wang S, Yan C. Robust optimal design of chilled water systems in buildings with quantified uncertainty and reliability for minimized life-cycle cost. Energy Build 2016;126:159–69. [17] Tian W, Meng X, Yin B, Sun Y, Liu Y, Fu X. Design of Robust Green Buildings Using a Non-probabilistic Uncertainty Analysis Method. Procedia Eng 2017;205:1049–55. [18] Fu X, Sun H, Guo Q, Guo Q, Pan Z, Xiong W, et al. Uncertainty analysis of an integrated energy system based on information theory. Energy 2017:649–62. [19] Khosravi A, Nahavandi S, Creighton D. Construction of optimal prediction intervals for load forecasting problems. IEEE Trans Power Syst 2010;25(3):1496–503. [20] Khosravi A, Nahavandi S, Creighton D, Atiya AF. Lower upper bound estimation method for construction of neural network-based prediction intervals. IEEE Trans Neural Networks 2011;22(3):337–46. [21] Walker WE, Harremoës P, Rotmans J, Van Der Sluijs JP, Van Asselt MBA, Janssen P,

[23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43]

403

et al. Defining uncertainty: a conceptual basis for uncertainty management in model-based decision support. Integr Assessm 2003;4(1):5–17. Li H, Wang S, Cheng H. Sensitivity analysis of design parameters and optimal design for zero/low energy buildings in subtropical regions. Appl Energy 2018:1280–91. Li CZ, Shi YM, Huang XH. Sensitivity analysis of energy demands on performance of CCHP system. Energy Convers Manage 2008;49(12):3491–7. Lu Y, Wang S, Yan C, Shan K. Impacts of renewable energy system design inputs on the performance robustness of net zero energy buildings. Energy 2015;93(2):1595–606. Mavromatidis G, Orehounig K, Carmeliet J. Uncertainty and global sensitivity analysis for the optimal design of distributed energy systems. Appl Energy 2018;214:219–38. Saltelli A, Tarantola S, Campolongo F, Ratto M. Sensitivity Analysis in Practice. J Am Stat Assoc 1989;101(473):398–9. Lu Y, Wang S, Huang Z. A robust design method for renewable energy system of zero energy consumption buildings. Refrigeration 2016;35(4):54–61. Cheng Q, Yan C, Wang S. Robust optimal design of chiller plants based on cooling load distribution ☆. Energy Procedia 2015;75:1354–9. Mavromatidis G, Orehounig K, Carmeliet J. Design of distributed energy systems under uncertainty: a two-stage stochastic programming approach. Appl Energy 2018;222:932–50. Li F, Wu T, Badiru A, Hu M, Soni S. A single-loop deterministic method for reliability-based design optimization. Eng Optim 2013;45(4):435–58. Hu M, Cho H. A probability constrained multi-objective optimization model for CCHP system operation decision support. Appl Energy 2014;116(116):230–42. Li G, Sun W, Huang GH, Lv Y, Liu Z, An C. Planning of integrated energy-environment systems under dual interval uncertainties. Int J Electr Power Energy Syst 2018;100:287–98. Hussain A, Arif SM, Aslam M, Shah SDA. Optimal siting and sizing of tri-generation equipment for developing an autonomous community microgrid considering uncertainties. Sustain Cities Soc 2017;32:318–30. Akbari K, Jolai F, Ghaderi SF. Optimal design of distributed energy system in a neighborhood under uncertainty. Energy 2016;116:567–82. Wang C, Jiao B, Guo L, Tian Z, Niu J. Robust scheduling of building energy system under uncertainty. Appl Energy 2016;167:366–76. Wang S, Wang K, Teng F, Strbac G, Wu L. An affine arithmetic-based multi-objective optimization method for energy storage systems operating in active distribution networks with uncertainties. Appl Energy 2018;223:215–28. https://www.gams.com/latest/docs/S_GUROBI.html. https://www.energyplus.net/. Mavromatidis G, Orehounig K, Carmeliet J. Designing electrically self-sufficient distributed energy systems under energy demand and solar radiation uncertainty. Energy Procedia 2017;122:1027–32. Tian Z, Niu J, Lu Y, He S, Tian X. The improvement of a simulation model for a distributed CCHP system and its influence on optimal operation cost and strategy. Appl Energy 2016;165(5):430–44. Design Standard for Energy Efficiency of Public Buildings GB 50189-2015. Ministry of housing and urban-rural construction of the people’s republic of China. 2015 [Beijing (in Chinese)]. Heo Y. Bayesian calibration of building energy models for energy retrofit decisionmaking under uncertainty. Dissert Theses – Gradworks 2011. Huang P, Huang G, Wang Y. HVAC system design under peak load prediction uncertainty using multiple-criterion decision making technique. Energy Build 2015;91:26–36.