A two-stage stochastic programming model for the optimal design of distributed energy systems

Applied Energy 103 (2013) 135–144

Contents lists available at SciVerse ScienceDirect

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

A two-stage stochastic programming model for the optimal design of distributed energy systems Zhe Zhou a, Jianyun Zhang a, Pei Liu a,⇑, Zheng Li a, Michael C. Georgiadis b, Efstratios N. Pistikopoulos c a

State Key Laboratory of Power Systems, Department of Thermal Engineering, Tsinghua University, Beijing 100084, China Department of Chemical Engineering, Aristotle University of Thessaloniki, Thessaloniki 54124, Greece c Centre for Process Systems Engineering (CPSE), Department of Chemical Engineering, Imperial College London, London SW7 2AZ, UK b

h i g h l i g h t s " The optimal design of distributed energy systems under uncertainty is studied. " A stochastic model is developed using genetic algorithm and Monte Carlo method. " The proposed system possesses inherent robustness under uncertainty. " The inherent robustness is due to energy storage facilities and grid connection.

a r t i c l e

i n f o

Article history: Received 24 May 2012 Received in revised form 17 August 2012 Accepted 10 September 2012 Available online 9 October 2012 Keywords: Stochastic programming Distributed energy system Uncertainty Genetic algorithm The Monte Carlo method

a b s t r a c t A distributed energy system is a multi-input and multi-output energy system with substantial energy, economic and environmental benefits. The optimal design of such a complex system under energy demand and supply uncertainty poses significant challenges in terms of both modelling and corresponding solution strategies. This paper proposes a two-stage stochastic programming model for the optimal design of distributed energy systems. A two-stage decomposition based solution strategy is used to solve the optimization problem with genetic algorithm performing the search on the first stage variables and a Monte Carlo method dealing with uncertainty in the second stage. The model is applied to the planning of a distributed energy system in a hotel. Detailed computational results are presented and compared with those generated by a deterministic model. The impacts of demand and supply uncertainty on the optimal design of distributed energy systems are systematically investigated using proposed modelling framework and solution approach. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction A distributed energy system is typically based on the concept of ‘‘local production of energy for local consumption’’. It refers to an advanced energy supply system which consists of many small scale energy generation technologies and locates at or near end-users [1]. Recently, interest has been intensifying in the development of distributed energy systems owing to their high overall efficiency, excellent environmental performance and other benefits [2–4]. However, balancing demand and supply is a challenging task for a distributed energy system, because it always faces instantaneously varying loads and small number of equipments within the system provide limited operation flexibility to cope with the fluctuation in energy demands [4]. Therefore, distributed energy systems may not be able to produce the potential benefits due to ⇑ Corresponding author. Tel.: +86 10 62795734x333; fax: +86 10 62795736. E-mail address: [email protected] (P. Liu). 0306-2619/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.apenergy.2012.09.019

lack of appropriate system configurations and suitable operation strategies [1,5]. To address issues related to the optimal design and operation of distributed energy systems, several mathematical models have been proposed using different mathematical programming techniques such as mixed-integer programming (MIP) and multiobjective programming (MOP) [6,7]. Most models focus on the combined heating cooling and power (CCHP) system, which is a representative type of distributed energy system [8–19], whereas a few models investigate more complex systems where renewable energy resources are introduced [1,5,20–23]. A common feature of most of these approaches is that they are formulated as deterministic mathematical programming models, ignoring the uncertainties of parameters. If model uncertainties are not adequately identified and handled, the actual economic and feasible operation of the designed distributed energy system may deviate from the optimal one. Thus, the assessment of uncertainty in the modelling of distributed energy systems has been recently received a lot of

136

Z. Zhou et al. / Applied Energy 103 (2013) 135–144

Nomenclature Parameters emax the maximum amount of primary energy available, MW D the maximum energy load, MW Inv unit investment cost, yuan/MW OM unit operation and maintenance cost, yuan/MW P fuel price, yuan/MW h Y lifetime of different technologies, year I discount rate, 10% in this study Continuous variables e energy flow, MW D energy demand, MW e+ energy stored to storage devices, MW energy released from storage devices, MW eE the inventory in storage devices, MW h Cap capacity of different technologies, MW Ctotal total annual cost, yuan capital cost, yuan Ccapital CO&M operation and maintenance cost, yuan Cfuel fuel cost, yuan g efficiency or COP TE total annual energy consumption, MW h Binary variables jheating on and off status of energy conversion technologies operating in heating mode jcooling on and off status of energy conversion technologies operating in cooling mode

attention. Gamou et al. proposed an optimal unit sizing method for cogeneration systems taking into consideration uncertainty in energy demands [24]. Yoshida et al. performed a sensitivity analysis on uncertainties caused by volatile equipment performance, energy prices and possible decline in equipment costs [25]. Li et al. conducted a sensitivity analysis on fluctuating energy demands to evaluate their influence on the performance of CCHP systems [26]. Smith et al. analyzed a CCHP system model under different operating strategies with uncertainty in inputs and parameters, such as energy demand and process efficiency [27]. Li et al. proposed a mathematical programming model to optimize a CCHP system under uncertainty in energy demands [28]; Mavrotas et al. developed an energy planning framework that can be used in buildings taking into account uncertainties in economic parameters [29]. There are different types of uncertainty in the optimal design of distributed energy systems, for example, uncertainty in energy demand and energy supply, uncertainty in economic parameters such as unit investment cost and energy price and uncertainty in technological parameters such as efficiency. In most of the aforementioned studies considering uncertainty in the planning of distributed energy systems, energy demands uncertainty catches the most significant attention among different sources of uncertainties. Energy demands uncertainty arises from the simplified description of energy demands patterns using 24-h profiles of representative days thus reducing the dimension of the underlying optimization problem and improving the computational efficiency. However, there is always a trade-off between the computational efficiency and the accuracy of modelling results. Much information about energy demands is lost due to the adoption of representative days, leading to less accurate model predictions. Hawkes and Leach and Mavrotas et al. concluded in their studies independently that grouping of demand data will jeopardize the accuracy of the optimization results [29,30], which implies that considering the

k

on and off status of distributed technologies

Subscripts f different types of fuel sec different types of secondary energy ter different types of tertiary energy d different types of demand m month hr hour egt energy generation technologies ect energy conversion technologies es energy storage technologies p electricity h heat c cooling hw hot water sh space heating gp grid power ng natural gas Superscripts g energy c energy s energy in energy out energy

generation section conversion section storage section input output

volatility of energy demands is of great importance towards the optimal design of CCHP system. In our previous work, a superstructure based mixed-integer linear programming (MILP) model has been proposed for the optimal design of distributed energy systems under purely deterministic conditions [31]. The salient feature of this model lies in its capability to take into account a large variety of potential energy resources and alternative technologies. This feature makes it possible to utilize distributed renewable energies that are utterly fluctuating. The method of representative days is also used for the description of energy demands as well as hourly resource availabilities of solar energy and wind energy. Therefore, lessons learned from previous works on the uncertainty of CCHP systems cannot be directly applied to more complex distributed energy systems proposed by our model. The external uncertainties should be addressed by considering fluctuations in both demand side and supply side. This work extends our previous contribution and presents a two-stage stochastic programming model for the optimal design of distributed energy systems with a stage decomposition based solution strategy. Both demand uncertainty and supply uncertainty are taken into consideration, making the uncertainty analysis in this study more comprehensive than previous studies which consider demand uncertainty only. In this model, genetic algorithm performs the search on the first stage variables while the Monte Carlo method is adopted in the second stage to deal with uncertainty in energy demand and supply. As the complex distributed energy system with multiple energy resources is a new concept proposed by us, the stochastic model, in which uncertainties in energy demand and supply are considered simultaneously, is novel in comparison to existing contributions. This paper is organized as follows: First a deterministic model for the optimal design of distributed energy systems is presented in Section 2. Then, a stage decomposition based solution strategy is briefly introduced

Z. Zhou et al. / Applied Energy 103 (2013) 135–144

in Section 3. The subsequent Section 4 provides a quantitative description of uncertainties of energy demands, solar energy availability and wind energy availability. Section 5 presents the results of the stochastic programming problem in comparison to those of the deterministic one. The impacts of external uncertainties on the optimal design of distributed energy systems are investigated based on the modelling results.

2. Mathematical model for distributed energy systems design under uncertainty 2.1. Deterministic model for the optimal design of distributed energy systems A superstructure based mixed-integer linear programming (MILP) model for the optimal design of distributed energy systems under purely deterministic conditions is presented in this section. The superstructure representation of a distributed energy system, i.e., its potential energy flow network, is shown in Fig. 1. The structure is comprised of an energy generation section, an energy conversion section and an energy storage section. The energy conversion processes take place in the following steps: Firstly, different primary energy resources are converted into heat and electricity by a variety of generation technologies in the energy generation section. Then, the secondary energy carriers, i.e., heat and electricity, are then converted to different forms of tertiary energy carriers, i.e., heat, cooling and electricity, by various conversion technologies in the energy conversion section. For example, water-source heat pump is a representative energy conversion technology that can produce heat with electricity. Next, energy storage technologies are installed to shift energy generated during low usage periods for consumption during peak hours aiming at a

137

more efficient use of energy. Finally, different forms of energy are supplied to end users to satisfy their energy demands. The framework includes twenty types of equipments, i.e., internal combustion engine, gas turbine, natural gas boiler, fuel cell, diesel engine, wind turbine, solar photovoltaic (PV), solar thermal collector, biomass boiler, biogas engine, biogas boiler, electric heater, absorption chiller/heat pump, heat exchanger, water-source chiller/heat pump, air-source chiller/heat pump, ground-source chiller/heat pump, battery, ice storage facility and thermal storage facility. These equipments are the most prevailing kinds of distributed energy technologies, while the model is capable for covering more technologies when they become commercially available. The inclusion of all possible distributed energy technologies does not imply that all of them will be adopted, but this provides various possible configurations for the model to choose from. The mathematic formulation of the deterministic programming problem as presented in our previous work is summarized in Appendix A [31]. The MILP problem can be recast in the following compact form:

min f ðb; d; xt Þ b;d;xt

s:t: udc ðb; dÞ ¼ 0 wdc ðb; dÞ  0

uoc ðb; d; xt Þ ¼ 0

ð1Þ

woc ðb; d; xt Þ  0 b 2 f0; 1gm ; d 2 R1 ; xt 2 Rn where b is a vector of binary design variables, representing the selection (or not) of types of equipments, d is a vector of continuous design variables, which represent continuous decisions to be made at the design stage, e.g., the capacities of equipments, xt denotes a vector of continuous operational variables, representing quantitative decisions to be made at any point in time t, f represents the

Fig. 1. The superstructure representation of a distributed energy system.

138

Z. Zhou et al. / Applied Energy 103 (2013) 135–144

objective function. In this work, the objective function is the total annual cost, udc and wdc are equality and inequality design constraints, which involve design variables b and d only, uoc and woc are equality and inequality operational constraints, which involve operational variables xt and/or design variables b and d. 2.2. Two-stage stochastic programming formulation The MILP design problem features simultaneous determination of design and operational variables. However, when uncertainty in energy demand and supply is considered, design variables and operational variables should be decided in different stages, because it cannot be anticipated what kinds of demand and supply profiles will realize and decisions on design variables must be able to adapt to various possible demand and supply profiles. As a result, binary design variables b and continuous design variables d are classified as the first stage variables that should be decided ‘‘here-and-now’’ prior to the resolution of uncertainty and operational variables xt that are ‘‘wait-and-see’’ variables in the second stage, which can be decided when all uncertain parameters have been observed. The optimization problem involves simultaneous determination of the optimal configuration of a distributed energy system and its optimal operating conditions. It is a typical two-stage recourse problem, which is usually addressed using the two-stage stochastic programming approach. By incorporating the uncertainty into the MILP problem, a twostage stochastic programming problem can be formulated as a master problem in combination with a sub-problem as follows:

ðMasterÞ : min f d ðb; dÞ þ E½fs ðb; d; nÞ b;d

s:t: udc ðb; dÞ ¼ 0 wdc ðb; dÞ  0 b 2 f0; 1gm ; d 2 R1 ðSubÞ : fs ðb; d; nÞ ¼ minfs ðb; d; xt ; nÞ

ð2Þ

xt

s:t: uoc ðb; d; xt ; nÞ ¼ 0

ered to be a suitable approach. Similarly, the master problem in this study is addressed by searching the solution space with genetic algorithm that uses techniques inspired by natural evolution, such as inheritance, crossover and mutation [8,44–46]. The evolutionary process guarantees offspring with good quality to survive to the next generation and eventually the algorithm has the ability to evolve to the optimum solution for the problem. During the search processes of the master problem, the recourse term is obtained through the second-stage sub-problem when the design variables are fixed. In this study, in contrast to scenario programming used in Till et al. [42], a sample average approximation method, i.e., the Monte Carlo (MC) method [47–49], is adopted to calculate the expected value of the stochastic term fs. Monte Carlo methods are a class of computational algorithms that rely on repeated random sampling to compute their results. These methods are most suited to calculation by a computer and tend to be used when it is infeasible to compute an exact result with a deterministic algorithm. In a Monte Carlo simulation, software samples a random value from each input distribution and runs the model using those values. After repeating the process a number of times, it estimates probability distributions for the uncertain outputs of the model from the random sample of output values. The larger the sample size, the more accurate the estimation of the output distributions. In a Monte Carlo based programming problem, uncertain parameters are given in the form of specific probability distributions rather than point values. The optimization process takes place iteratively with distinct values of uncertain parameters generated according to their probability distributions. After a large number of Monte Carlo runs, the results for the objective function and the main decision variables are obtained and recorded in the form of probability distributions. According to law of large numbers, the average of the results obtained from a large number of trials should be close to the expected value, and will tend to become closer as more trials are performed [50]. The stage decomposition based hybrid evolutionary approach for solving the two-stage stochastic programming problem is shown in Fig. 2. The algorithm is briefly described as follows:

woc ðb; d; xt ; nÞ  0 x t 2 Rn

Algorithm.

where, the objective function is split into a deterministic term fd representing decisions at the design stage and the expectation of a stochastic term fs which depends on the realization of uncertain parameters n at the operation stage.

Initialize the multi-set (population) bi and di, where the subscript i refers to the ith individual. The initial population bi and di are generated randomly. Evaluate the initial individual applying the fitness function F. For simplification purposes, the objective function, i.e., the total annual cost, is selected as the fitness function, given by

2.3. Solution strategy for the two-stage stochastic optimization problem

F ¼ fd ðb; dÞ þ E½fs ðb; d; nÞ Two-stage stochastic programming has been applied to various fields, such as supply chain planning [32–35], process design and operation [36–39] and infrastructure planning [40,41]. In this paper, a hybrid evolutionary algorithm, proposed by Till et al. [42], is adopted to solve the two-stage stochastic optimization problem. Their proposed method features that an evolutionary algorithm performs the search on the first-stage variables while the second-stage sub-problems are solved by a mixed-integer programming solver. Tometzki and Engell investigated the performance of the hybrid evolutionary algorithm and of the standard L-shape method with respect to the problem size [43]. For small scale problems, these two approaches exhibit similar performance in terms of convergence rate and solution quality. However, for larger problem sizes, the hybrid evolutionary algorithm is able to find optimal solutions within a short period. For long computation times, the hybrid evolutionary algorithm does not reach the solution quality of CPLEX. Due to the large size of the optimization problem in this work, the hybrid evolutionary algorithm is consid-

ð3Þ

While termination criteria not fulfilled in the kth generation, do Select best-ranking individuals to reproduce. The selection criterion is that the worst individuals, i.e., with large values of fitness function F, are replaced by the best individual. Breed new generation through crossover and/or mutation and give birth to offspring. The crossover operation is undertaken for each individual by

If Pc1 < Pcrossov er ; kþ1

¼ bj

kþ1

¼ Pc2 di þ ð1  Pc2 Þdj ;

kþ1

¼ ð1  Pc2 Þdi þ Pc2 dj

bi

di

dj

kþ1

k

k

¼ maxðbi ; bj Þ; k

k

k

ð4Þ

k

where Pc1 and Pc2 are both random numbers between 0 and 1, generated by a uniform random number generator; Pcrossover is the threshold of crossover operations.

139

Z. Zhou et al. / Applied Energy 103 (2013) 135–144

where N refers to the number of iterations in a Monte Carlo simulation. For each set of uncertainty parameters n, minxt fs;N ðb; d; xt ; nÞ is provided by a standard MILP solver. 3. Quantitative description of uncertain parameters In this study, three types of parameters, with significant volatility, are selected as uncertain parameters, i.e., energy demands, solar energy availability and wind energy availability. The probability distributions of these parameters are presented in this section. 3.1. Energy demands According to Gamou et al., it is assumed that probability distributions of energy demands at each sampling time roughly obey normal distribution in which 95% of the whole area is within the range of ±20% of the average energy demands [24]. Li et al. [26] and Li et al. [28] applied this assumption to their studies on CCHP systems. Therefore, it is also assumed that each type of energy demands at any point in time follows a normal distribution h  N(l, r2), where l represents the value of each point on the demand profile and 1.96r equals 20% of l because a 95% of confidence interval of normal distribution is corresponding to 1.96r. 3.2. Solar energy availability Methodologies for the estimation of hourly global solar radiation have been proposed by many researchers [51–54]. Kaplanis and Kaplani carried out a remarkable analysis about the stochastic characteristics of hourly global solar radiation with their observed data [55]. Implications drawn from their studies are as follow:  Fluctuations in hourly global solar radiation are strong in winter but weak in summer.  Fluctuations in hourly global solar radiation are strong for morning and late afternoon but weak around the solar noon.  Hourly global solar radiation follows a normal distribution. However, the eigenvalues vary at different point in time. In this study, the probability distribution of hourly global solar radiation is taken from Kaplanis and Kaplani [55] as listed in Table 1. Fig. 2. A flow chart of the stage decomposition based hybrid evolutionary approach.

The mutation operation is undertaken for each individual by

If Pm1 < Pmutation ; kþ1

bi

kþ1 di

k

k

¼ maxðbi ; bbest Þ; ¼

k di

þ 0:1ðP m2 

ð5Þ k 0:5Þdbest

where Pm1 and Pm2 are both random numbers between 0 and 1, generated by a uniform random number generator; Pmutation is the k k threshold of mutation operations; bbest and dbest represent the best individual in the kth generation.

3.3. Wind energy availability The probability density function (PDF) is a widely used approach to depict wind speed frequency distribution. A variety of PDFs have been proposed for the estimation of wind speed [56], among which the two parameter Weibull distribution is the most frequently used and widely accepted distribution for wind energy [57–60]. In this study, a three parameter generalized extreme value (GEV) distribution is adopted, which is a more general form of the two parameter Weibull distribution [61]. This PDF was fitted

Table 1 Probability distributions of hourly global solar radiations.

Evaluate the offspring applying the fitness function F. End while In the evaluation operations, the expected value of fs is given by

E½fs ðb; d; nÞ ¼

1X minfs;N ðb; d; xt ; nÞ N N xt

ð6Þ a

Time frame

Distributiona

Standard deviation, r

November–April, 9:00 am–15:00 pm November–April, the rest of the day May–October, 9:00 am–15:00 pm May–October, the rest of the day

h  N(l, r2) h  N(l, r2) h  N(l, r2) h  N(l, r2)

12%l 25%l 3%l 8%l

l represents the value of each point on the solar radiation profile.

140

Z. Zhou et al. / Applied Energy 103 (2013) 135–144

with Chinese wind resource data, which makes it suitable for our analysis. The cumulative distribution function is as follows:

 h x  li1=k  Fðx; l; r; nÞ ¼ exp  1 þ k

ð7Þ

r

where l is the location parameter, r the scale parameter and k the shape parameter; l equals 0.654; r equals 0.498; k equals 0.081. Many types of software have built-in library functions that can generate basic distributions, such as uniform distribution and normal distribution. However, specific techniques are desired to generate random numbers that follows an unconventional distribution. A commonly used technique is called the inverse transform technique. Suppose that a series of random numbers have a density function f. Let F denote the corresponding cumulative distribution function (CDF). If X = F1(U), then X is a random variable with CDF FX(x) = F, because:

F X ðaÞ ¼ PðX  aÞ ¼ PðF 1 ðUÞ  aÞ ¼ PðU  FðaÞÞ ¼ FðaÞ

ð8Þ

By this means, fluctuating wind energy can be simply simulated on the basis of a uniform random number generator. 4. Case study With the optimization model and its solution strategy described in Section 2 and probability distributions of energy demands, solar

energy availability and wind energy availability provided in Section 3, the deterministic and stochastic models are applied to the planning of a distributed energy system for a hotel in Beijing. The hourly demand data of this hotel in some representative days are shown in Fig. 3. Two cases are investigated for comparison purposes. Case 1 is a deterministic programming problem which takes energy demands, solar energy availability and wind energy availability of representative days as the model inputs. Case 2 is a two-stage stochastic programming problem for evaluation of the impact of multiple uncertainties. For the first stage optimization, population size for GA is 8, crossover rate is 0.7, and mutation rate is 0.1. When the number of generations exceeds 100, searching will stop automatically. For the second stage, the sample size is set 100 for each Monte Carlo simulation. Assumptions, data, and other details of the case studies can be obtained from our earlier publication [31]. All case studies have been solved using the GAMS Modelling and Optimization Platform [62]. The deterministic MILP problem (Case 1) involves 13,377 continuous variables, 6624 binary variables and 28,906 constraints. It takes 0.25 s of CPU time to solve this problem on an Intel(R) Core(TM) i3-2100 CPU of 3.10 GHz with a memory size of 2.99G. Therefore, it can be deduced that each Monte Carlo run in the stochastic MILP problem requires less than 1 s of CPU time. However, due to the large number of generations in the first stage and Monte Carlo simulations in the second stage, it takes about 14 h to find the optimal solution.

Fig. 3. The hourly energy loads of the hypothetic hotel in Beijing.

141

Z. Zhou et al. / Applied Energy 103 (2013) 135–144

5. Results and discussions The solution strategy for the two-stage stochastic optimization problem is a sample based algorithm, therefore it must employ a large number of samples in the Monte Carlo method and adequate number of generations in the genetic algorithm to guarantee the accuracy of the modelling results. The number of samples in each Monte Carlo simulation should be large enough so that the expected value of the stochastic term fs as shown in Eq. (3) does not change significantly with the sample size. Fig. 4 shows the average total annual cost under different sample sizes in one of the Monte Carlo simulations (Generation: 1st, Population: 1st). The average total annual cost converges quickly as the number of samples increases, therefore, a hundred samples for each Monte Carlo simulation are sufficient to ensure the accuracy of the modelling. The generation size is another crucial factor for provision of credible results. Fig. 5 shows the minimal total annual cost of each

Fig. 4. The average total annual cost under different sample sizes.

Fig. 5. The minimal total annual in each generation.

generation in the genetic algorithm. The historical minimal cost drops gradually and finally converges to 7.453 million RMB yuan, which is considered to be an optimal solution. The evolution of the total annual cost shown in Fig. 5 indicates that one hundred generations are adequate for accurate results. The obtained optimal design under uncertainty and the optimal design under deterministic conditions are shown together in Table 2 for comparison. Variables investigated include the objective function, i.e., total annual cost and key decision variables, i.e., capacities of various equipments that are selected by the optimization model, including internal combustion engine, wind turbine, solar heat collector, air source chiller/heat pump, water source chiller/heat pump, absorption chiller/heat pump, thermal storage facility and ice storage facility. It can be observed that the optimal design under uncertainty does not significantly deviate from the deterministic optimal design. Compared with the deterministic case, the minimum expected total annual cost of the stochastic case is brought up slightly and most devices exhibit minute differences in their capacities between the deterministic design and the stochastic design. The small difference between the stochastic design and the deterministic design mainly stems from the introduction of energy storage technologies and grid connection. Energy storage technologies can shift extra energy generated at one point in time for consumption at another point in time. Therefore, the system is capable to switch its operation schedule of energy storage facilities to adjust to different types of demand and supply profiles, leaving the core devices working in a nearly fixed mode. The grid power serves as a complementary power source. When there is a sudden increase in electricity consumption, the core power generation facilities can keep working in a fixed mode and the electricity shortage can be met by introducing more electricity from the power grid. To demonstrate the effects of installation of energy storage devices and grid connection, a separate case study was performed, in which installation of energy storage devices and grid connection are prohibited. The problem is solved with the same algorithm developed in this study. The results of the additional case study are presented in Table 3. In this case, the total annual cost and most equipment capacities of the stochastic case become much larger than those of the deterministic case. The different results between the cases with and without energy storage facilities and grid connection indicate that the introduction of energy storage and grid connection can reduce the negative impacts of energy demand and supply uncertainties. The small difference between the results of deterministic and stochastic models also indicates that a distributed energy system with energy storage facilities and grid connection possesses inherent robustness in dealing with demand and supply uncertainty, which makes such type of systems favored by applications with significant demand and supply uncertainties.

Table 2 Results comparison between deterministic and stochastic programming problems. Deterministic

Stochastic (expected value)

Objective function [cost, million RMB yuan] Total annual cost

7.359

7.453

Key design variables [capacity of equipment, MW] Internal combustion engine Wind turbine Solar heat collector Air source heat pump/chiller Water source heat pump/chiller Absorption heat pump/chiller Thermal storage facility Ice storage facility

1.425 0.038 0.234 0.682 2.476 2.584 2.565 4.541

1.425 0.038 0.247 1.037 2.490 2.584 3.811 6.125

Difference (%) 1.28 0.0 0.0 5.9 52.1 0.6 0.0 48.6 34.9

142

Z. Zhou et al. / Applied Energy 103 (2013) 135–144

Table 3 Results comparison between deterministic and stochastic programming problems when energy storage technologies and grid connection are prohibited. Deterministic

Stochastic (expected value)

Objective function [cost, million RMB yuan] Total annual cost

8.416

8.755

Key design variables [capacity of equipment, MW] Internal combustion engine Gas boiler Wind turbine Solar heat collector Air source heat pump/chiller Water source heat pump/chiller Absorption heat pump/chiller

2.014 0.623 0.040 0.361 1.394 1.217 4.379

2.152 1.001 0.060 0.459 0.905 2.664 5.320

For the optimal design of distributed energy systems, the system configuration and capacities of equipment are very similar for the deterministic and stochastic results. Therefore, redundancy design for coping with uncertainty is almost dispensable for distributed energy systems with energy storage facilities and grid connection. In light of the high computational efficiency of the deterministic model, it is preferred for the optimal design of distributed energy systems with fluctuating energy demand and supply. 6. Conclusion This work presents the development and application of a twostage stochastic model for the optimal design of distributed energy systems under energy demand and supply uncertainty. The potential impacts of uncertain parameters are accounted by reformulating a deterministic design problem into a two-stage stochastic programming problem. A solution strategy which combines genetic algorithm and the Monte Carlo method is proposed to solve the stochastic programming problem. A case study is used to illustrate that the optimal design under uncertainty does not significantly deviate from the deterministic one. The small difference between the stochastic design and the deterministic design mainly stems from the introduction of energy storage technologies and grid connection. The proposed distributed energy system, with energy storage facilities and grid connection, possesses inherent robustness under uncertainty. It makes this type of system suitable by applications with demand and supply uncertainty. The deterministic model is preferred for the optimal design of distributed energy systems with fluctuating energy demand and supply due to its computational efficiency. Acknowledgements The authors gratefully acknowledge the financial support from National Natural Science Foundation (Project No. 51106080), from the IRSES ESE Project of FP7 (Contract No: PIRSES-GA-2011294987), from BP company in the scope of the Phase II Collaboration between BP and Tsinghua University and from Mitsubishi Heavy Industries.

Difference (%) 4.03 6.9 60.6 51.0 27.2 35.1 119.0 21.5

A.1. Energy generation section In this section, various types of primary energy are converted to electricity and heat by different energy generation technologies. The primary energy consumptions are restrained by local resources availability at any point in time as follows: max eg;in f ;m;hr  ef ;m;hr

ðA:1Þ

where eg;in refers to the primary energy consumption; emax is the f f maximum amount of energy available for energy use; f stands for different types of primary energy resources, including natural gas, diesel, wind energy, solar energy and biomass; m stands for different months in a year and hr stands for different hours in a day. Each type of primary energy can be consumed by any technology that uses it as the fuel. The total consumption of each type of primary energy resource at any point in time is given by

eg;in f ;m;hr ¼

X

eg;in egt;f ;m;hr

ðA:2Þ

egt

where the subscript egt stands for different types of energy generation technologies. The production of secondary energy carriers, i.e., electricity and heat, by various energy generation technologies depends on the consumptions of primary energy resources, and the efficiency of the on-site generation technologies, which can be characterized as follows:

eg;out sec;m;hr ¼

X g;in eegt;m;hr  gegt;sec

ðA:3Þ

egt

where sec stands for secondary energy carriers, which could be either electricity or heat; sec e {p, h}, p and h refer to electricity g;out and heat; esec refers to the production of secondary energy; gegt,sec refers to the efficiency of converting a primary energy resource into a secondary energy carrier. Electricity can also be directly purchased from the power grid. Therefore, the total amount of electricity as input to the energy conversion section is the sum of electricity drawn from the grid and the electricity generated in the energy generation section. This can be denoted as

Appendix A. Mathematical formulation of the MILP model

g;out gp ec;in p;m;hr ¼ ep;m;hr þ ep;m;hr

A detailed MILP model for distributed energy systems design, based on the superstructure modelling approach as shown in Fig. 1, is presented here. Symbol e represents all energy flows and its superscripts indicate its position in the superstructure while the subscripts represent the types of energy flows. With respect to superscripts, the first superscript, i.e., g, c and s, stands for the section that this energy flow belongs to, while the second one, i.e., in and out, represents whether it is an input or an output.

The other type of energy carrier, heat, can only be generated by different energy generation technologies, as follows: g;out ec;in h;m;hr ¼ eh;m;hr

ðA:4Þ

ðA:5Þ

where ec;in sec is the amount of secondary energy for conversion in the energy conversion section; egp p denotes the amount of electricity drawn from the power grid.

143

Z. Zhou et al. / Applied Energy 103 (2013) 135–144

A.2. Energy conversion section In the energy conversion section, the consumption of each type of secondary energy is equal to the sum of energy consumed by energy conversion technologies which convert this type of secondary energy, as follows:

ec;in sec;m;hr

¼

X

ec;in ect;sec;m;hr

ðA:6Þ

Finally, each type of end use energy demand product must be satisfied by its corresponding tertiary energy carrier coming from the energy storage section, which is given by

es;out h;m;hr  Dsh;m;hr þ Dhw;m;hr es;out p;m;hr  Dp;m;hr es;out c;m;hr

ðA:13Þ

 Dc;m;hr

ect

where ect stands for different types of energy conversion technologies. Energy conversion technologies such as air-source chiller/heat pump usually have two working modes: heating mode and cooling mode. The coefficient of performance (COP) of an energy conversion technology differs when working in the heating mode and cooling mode. The productions of these tertiary energy carriers can be given by

ec;out ter;m;hr ¼

X

heating

ec;in ect;sec;m;hr  ðjm;hr

cooling

 gheating  gcooling ect;sec;ter þ jm;hr ect;sec;ter Þ

ðA:7Þ

ect

where ect stands for different types of energy conversion technologies; ter stands for tertiary energy carriers, i.e., electricity, heat and cooling; ter e {p, h, c}, c refers to cooling; ec;out refers to the producter tion of tertiary energy in the energy conversion section; gheating ect;sec;ter and gcooling ect;sec;ter are COPs of energy conversion technologies operating heating cooling in heating mode and cooling mode; jm;hr and jm;hr are binary variables. In this equation, a nonlinear term is introduced due to the multiplication of a binary variable and a continuous variable. For model simplification purposes, a model linearization procedure is performed to convert the nonlinear term into linear terms. As an energy conversion technologies cannot produce heat and heating cooling cooling simultaneously, jm;hr and jm;hr are restrained by heating

jm;hr

cooling

þ jm;hr

¼1

ðA:8Þ

A.3. Energy storage section

¼

ec;out ter;m;hr

ðA:9Þ

where es;in ter denotes the input of an energy storage equipment. Each end use energy product is equal to the product generated in the energy conversion section plus the net energy flow from the energy storage unit, namely energy stored in the storage unit less released energy, at any point in time, which is denoted as

es;out ter;m;hr

¼

es;in ter;m;hr

þ

eter;m;hr



eþter;m;hr

es;out ter

ðA:10Þ e ter

where refers to the output of an energy storage device; refers to energy released from the energy storage device; eþ ter refers to energy stored into the energy storage device. The amount of energy stored at the beginning of each time interval is equal to the energy stored at the beginning of the previous time interval plus the net energy flow, given by

Eter;m;hrþ1 ¼ Eter;m;hr þ eþter;m;hr  eter;m;hr =ges;ter  0

ðA:11Þ

where Eter represents the inventory of a storage device; ges,ter refers to the efficiency of an energy storage technology. Within a certain period of time, T, the total net energy flow from the energy storage units sums to zero as given by T X Eter;m;hr ¼ 0

A.4. Facility capacity The capacities of all energy generation, storage and conversion technologies must be large enough to meet the maximum amounts of energy generation, conversion and storage over the whole period of time, which are given by

Captech  max eout tech

ðA:14Þ

where tech refers to each technology appears in the energy generation section and energy conversion section, i.e., tech e {egt, ect}. The outputs of different technologies are constrained within their lower and upper bounds of capacity via introduction of binary variables, which can ensure that these technologies are running in reasonable ranges, as follows: Max out ktech;m;hr  CapMin tech  etech;m;hr  ktech;m;hr  Captech

ðA:15Þ

where k is binary variable; CapMin and CapMax are minimal and maximal capacities of different technologies. The capacity of energy storage technologies depends on the maximum inventory of a storage device, which can be denoted as:

Capes  max Eter;m;hr

ðA:16Þ

A.5. Objective function

The energy inputs of the energy storage section are given by

es;in ter;m;hr

where D refers to end user’s demands; sh and hw represent space heating and hot water.

ðA:12Þ

hr

T is set as a day, which means the energy storage technologies can only deal with short-term, daily fluctuations.

The total annual cost is selected as the objective function, consisting of the capital cost, the operation & maintenance (O&M) cost and the fuel cost. The optimization model obtains the most suitable configuration of distributed energy system by minimizing the total annual cost, as follows:

Min C total ¼ C capital þ C O&M þ C fuel

ðA:17Þ

All the associated costs are given by

X Inv tech  Captech 

C capital ¼ C O&M

m

tech

C fuel

I

1  ð1 þ IÞY tech tech X XX g X ¼ ðOMtech  etech;m;hr Þ þ OMes  Capes hr

m

ðA:19Þ

es

XXX g;in ¼ ðef ;m;hr  Pf ;hr þ egp m;hr  P gp;hr Þ f

ðA:18Þ

ðA:20Þ

hr

where Invtech refers to the unit investment cost; Ytech is the lifetime; I denotes the discount rate, specified at 10% in this study; OMtech and OMes represent the unit operation costs; Pf,hr and Pgp,hr refer to the hourly fuel cost and electricity price. References [1] Ren H, Zhou W, Nakagami Ki, Gao W, Wu Q. Multi-objective optimization for the operation of distributed energy systems considering economic and environmental aspects. Appl Energy 2010;87:3642–51. [2] Akorede MF, Hizam H, Pouresmaeil E. Distributed energy resources and benefits to the environment. Renew Sust Energy Rev 2010;14:724–34. [3] Angel AB-R. Future development of the electricity systems with distributed generation. Energy 2009;34:377–83.

144

Z. Zhou et al. / Applied Energy 103 (2013) 135–144

[4] Pepermans G, Driesen J, Haeseldonckx D, Belmans R, D’haeseleer W. Distributed generation: definition, benefits and issues. Energy Policy 2005;33:787–98. [5] Ren H, Gao W. A MILP model for integrated plan and evaluation of distributed energy systems. Appl Energy 2010;87:1001–14. [6] Liu P, Pistikopoulos EN, Li Z. Energy systems engineering: methodologies and applications. Front Energy Power Eng China 2010;4:131–42. [7] Liu P, Georgiadis MC, Pistikopoulos EN. Advances in energy systems engineering. Ind Eng Chem Res 2011;50:4915–26. [8] Wang J, Jing Y, Zhang C. Optimization of capacity and operation for CCHP system by genetic algorithm. Appl Energy 2010;87:1325–35. [9] Kong XQ, Wang RZ, Huang XH. Energy optimization model for a CCHP system with available gas turbines. Appl Therm Eng 2005;25:377–91. [10] Abdollahi G, Meratizaman M. Multi-objective approach in thermoenvironomic optimization of a small-scale distributed CCHP system with risk analysis. Energy Build 2011;43:3144–53. [11] Mago PJ, Chamra LM. Analysis and optimization of CCHP systems based on energy, economical, and environmental considerations. Energy Build 2009;41:1099–106. [12] Roque Díaz P, Benito YR, Parise JAR. Thermoeconomic assessment of a multiengine, multi-heat-pump CCHP (combined cooling, heating and power generation) system – a case study. Energy 2010;35:3540–50. [13] Li H, Nalim R, Haldi PA. Thermal-economic optimization of a distributed multigeneration energy system—a case study of Beijing. Appl Therm Eng 2006;26: 709–19. [14] Carvalho M, Serra LM, Lozano MA. Optimal synthesis of trigeneration systems subject to environmental constraints. Energy 2011;36:3779–90. [15] Carvalho M, Lozano MA, Serra LM. Multicriteria synthesis of trigeneration systems considering economic and environmental aspects. Appl Energy 2012;91:245–54. [16] Rong A, Lahdelma R. An efficient linear programming model and optimization algorithm for trigeneration. Appl Energy 2005;82:40–63. [17] Arosio S, Guilizzoni M, Pravettoni F. A model for micro-trigeneration systems based on linear optimization and the Italian tariff policy. Appl Therm Eng 2011;31:2292–300. [18] Arcuri P, Florio G, Fragiacomo P. A mixed integer programming model for optimal design of trigeneration in a hospital complex. Energy 2007;32: 1430–47. [19] Ortiga J, Bruno JC, Coronas A. Selection of typical days for the characterisation of energy demand in cogeneration and trigeneration optimisation models for buildings. Energy Convers Manage 2011;52:1934–42. [20] Ren H, Zhou W, Gao W. Optimal option of distributed energy systems for building complexes in different climate zones in China. Appl Energy 2012;91: 156–65. [21] Ren H, Gao W, Zhou W, Nakagami Ki. Multi-criteria evaluation for the optimal adoption of distributed residential energy systems in Japan. Energy Policy 2009;37:5484–93. [22] Liu P, Pistikopoulos EN, Li Z. An energy systems engineering approach to the optimal design of energy systems in commercial buildings. Energy Policy 2010;38:4224–31. [23] NREL. HOMER – the optimization model for distributed power. [accessed 2012]. [24] Gamou S, Yokoyama R, Ito K. Optimal unit sizing of cogeneration systems in consideration of uncertain energy demands as continuous random variables. Energy Convers Manage 2002;43:1349–61. [25] Yoshida S, Ito K, Yokoyama R. Sensitivity analysis in structure optimization of energy supply systems for a hospital. Energy Convers Manage 2007;48: 2836–43. [26] Li C, Shi Y, Huang X. Sensitivity analysis of energy demands on performance of CCHP system. Energy Convers Manage 2008;49:3491–7. [27] Smith A, Luck R, Mago PJ. Analysis of a combined cooling, heating, and power system model under different operating strategies with input and model data uncertainty. Energy Build 2010;42:2231–40. [28] Li C, Shi Y, Liu S, Zheng Z, Liu Y. Uncertain programming of building cooling heating and power (BCHP) system based on Monte-Carlo method. Energy Build 2010;42:1369–75. [29] Mavrotas G, Florios K, Vlachou D. Energy planning of a hospital using mathematical programming and Monte Carlo simulation for dealing with uncertainty in the economic parameters. Energy Convers Manage 2010;51: 722–31. [30] Hawkes A, Leach M. Impacts of temporal precision in optimisation modelling of micro-combined heat and power. Energy 2005;30:1759–79. [31] Zhou Z, Liu P, Li Z, Ni W. An engineering approach to the optimal design of distributed energy systems in China. Appl Therm Eng 2012.

[32] Gupta A, Maranas CD. A two-stage modeling and solution framework for multisite midterm planning under demand uncertainty. Ind Eng Chem Res 2000;39:3799–813. [33] Clay R, Grossmann I. A disaggregation algorithm for the optimization of stochastic planning models. Comput Chem Eng 1997;21:751–74. [34] You F, Wassick JM, Grossmann IE. Risk management for a global supply chain planning under uncertainty: models and algorithms. AIChE J 2009;55:931–46. [35] Mitra K, Gudi RD, Patwardhan SC, Sardar G. Midterm supply chain planning under uncertainty: a multiobjective chance constrained programming framework . Ind Eng Chem Res 2008;47:5501–11. [36] Ierapetritou MG, Pistikopoulos EN. Batch plant design and operations under uncertainty. Ind Eng Chem Res 1996;35:772–87. [37] Liu P, Pistikopoulos EN, Li Z. Decomposition based stochastic programming approach for polygeneration energy systems design under uncertainty. Ind Eng Chem Res 2010;49:3295–305. [38] Guillén G, Espuña A, Puigjaner L. Addressing the scheduling of chemical supply chains under demand uncertainty. AIChE J 2006;52:3864–81. [39] Terrazas-Moreno S, Flores-Tlacuahuac A, Grossmann IE. Simultaneous design, scheduling, and optimal control of a methyl-methacrylate continuous polymerization reactor. AIChE J 2008;54:3160–70. [40] Han JH, Lee IB. A two-stage stochastic programming model for planning CO2 utilization and disposal infrastructure considering the uncertainty in the CO2 emission. Ind Eng Chem Res 2011. [41] Tarhan B, Grossmann IE, Goel V. Stochastic programming approach for the planning of offshore oil or gas field infrastructure under decision-dependent uncertainty. Ind Eng Chem Res 2009;48:3078–97. [42] Till J, Sand G, Urselmann M, Engell S. A hybrid evolutionary algorithm for solving two-stage stochastic integer programs in chemical batch scheduling. Comput Chem Eng 2007;31:630–47. [43] Tometzki T, Engell S. Hybrid evolutionary optimization of two-stage stochastic integer programming problems: an empirical investigation. Evol Comput 2009;17:511–26. [44] Renner G, Ekárt A. Genetic algorithms in computer aided design. Comput Aid Des 2003;35:709–26. [45] Goldberg DE. Genetic algorithms in search, optimization, and machine learning. Addison-wesley; 1989. [46] Shopova EG, Vaklieva-Bancheva NG. BASIC – a genetic algorithm for engineering problems solution. Comput Chem Eng 2006;30:1293–309. [47] Robert CP, Casella G. Monte Carlo statistical methods. Springer Verlag; 2004. [48] van Peborgh Gooch JR, Hounslow MJ. Monte Carlo simulation of sizeenlargement mechanisms in crystallization. AIChE J 1996;42:1864–74. [49] Diwekar UM, Kalagnanam JR. Efficient sampling technique for optimization under uncertainty. AIChE J 1997;43:440–7. [50] Durrett R. Probability: theory and examples. Cambridge Univ Pr; 2010. [51] Kaplanis S. New methodologies to estimate the hourly global solar radiation; comparisons with existing models. Renew Energy 2006;31:781–90. [52] Cao J, Lin X. Study of hourly and daily solar irradiation forecast using diagonal recurrent wavelet neural networks. Energy Convers Manage 2008;49: 1396–406. [53] Janjai S, Pankaew P, Laksanaboonsong J. A model for calculating hourly global solar radiation from satellite data in the tropics. Appl Energy 2009;86:1450–7. [54] El-Sebaii A, Al-Hazmi F, Al-Ghamdi A, Yaghmour SJ. Global, direct and diffuse solar radiation on horizontal and tilted surfaces in Jeddah, Saudi Arabia. Appl Energy 2010;87:568–76. [55] Kaplanis S, Kaplani E. A model to predict expected mean and stochastic hourly global solar radiation I (h; nj) values. Renew Energy 2007;32:1414–25. [56] Carta J, Ramírez P, Velázquez S. A review of wind speed probability distributions used in wind energy analysis: case studies in the Canary Islands. Renew Sust Energy Rev 2009;13:933–55. [57] Keyhani A, Ghasemi-Varnamkhasti M, Khanali M, Abbaszadeh R. An assessment of wind energy potential as a power generation source in the capital of Iran, Tehran. Energy 2010;35:188–201. [58] Ucar A, Balo F. Evaluation of wind energy potential and electricity generation at six locations in Turkey. Appl Energy 2009;86:1864–72. [59] Carta J, Ramirez P. Analysis of two-component mixture Weibull statistics for estimation of wind speed distributions. Renew Energy 2007;32:518–31. [60] Islam M, Saidur R, Rahim N. Assessment of wind energy potentiality at Kudat and Labuan, Malaysia using Weibull distribution function. Energy 2011. [61] Chen Z, Ni W, Li Z. Preliminary study on wind power characteristics. Acta Energiae Solaris Sinica 2011;32:210–5 [in Chinese]. [62] GAMS Development Corporation. GAMS – A user’s guide. Washington, DC; 2008. [accessed 2012].