Optimal design of hybrid renewable energy systems using simulation optimization

Simulation Modelling Practice and Theory 52 (2015) 40–51

Contents lists available at ScienceDirect

Simulation Modelling Practice and Theory journal homepage: www.elsevier.com/locate/simpat

Optimal design of hybrid renewable energy systems using simulation optimization Kuo-Hao Chang a,⇑, Grace Lin b a b

Department of Industrial Engineering and Engineering Management, National Tsing Hua University, Hsinchu, Taiwan Advanced Research Institute, Institute for Information Industry, Taipei, Taiwan

a r t i c l e

i n f o

Article history: Received 15 September 2014 Received in revised form 5 December 2014 Accepted 6 December 2014

Keywords: Hybrid renewable energy system Energy management Monte Carlo simulation Simulation optimization

a b s t r a c t A hybrid renewable energy system (HRES) provides a viable solution for electrification of remote areas when grid extension is impossible or uneconomical. Such a system has several power stations and each of which includes photovoltaics (PV) and wind power generators, along with a diesel power generator as backup when renewable energy supply is insufficient. While the HRES is attractive due to the minimal environmental and health impact compared to fossil fuels, the design of HRES, specifically the determination of the size of PV, wind, and diesel power generators and the size of energy storage system in each power station, is very challenging. This is mainly due to a large number of factors involved in the problem, the profound uncertainty arising from the renewable resources and the demand load, and the complex interaction among factors. In this paper, we investigate the use of Monte Carlo simulation, along with simulation optimization techniques, for obtaining the optimal design of HRES in uncertain environments. The proposed model considers not only the equipment installation, including PV, wind and diesel power generators and the energy storage systems in each power station, but also the power generation, allocation and transmission within the HRES so as to achieve minimum expected total cost, while satisfying the power demand. An extensive computational study shows that the proposed model of realistic size can be solved efficiently, enabling quality decisions to be generated in practice. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction Renewable energy, such as photovoltaics (PV) and wind power, is viewed as green and clean because they are non-depletable, non-polluting and have minimal environmental and health impacts compared to fossil fuels. As the price of fossil fuel continues to increase, renewable energy has become more popular over the decades. However, renewable energy is highly volatile and unpredictable due to the unknown and variable sunshine hours and wind speeds from day to day [12]. This raises some concern about the use of renewable energy as it may lead to unstable power supplies and inability to satisfy power demand in some areas. Consequently, according to Renewables 2013 Global Status Report, renewable energy only accounts for 19% of the total world energy demand in 2013. To this end, many power system designers started to look for other alternatives. A hybrid renewable energy system (HRES), which combines several power stations and each of which includes photovoltaics (PV) and wind power generators, along with a diesel power generator as backup when renewable energy supply is insufficient, is one viable and attractive solution that brings into play the advantages, while at the same ⇑ Corresponding author. http://dx.doi.org/10.1016/j.simpat.2014.12.002 1569-190X/Ó 2014 Elsevier B.V. All rights reserved.

K.-H. Chang, G. Lin / Simulation Modelling Practice and Theory 52 (2015) 40–51

41

time eliminating the disadvantages of renewable energy. Specifically, the HRES first capitalizes on the renewable energy and, if the renewable energy supply is not sufficient, the diesel power generator is used to cover the shortage, thereby ensuring the stability of the power supply. In the literature, extensive studies have been conducted on the feasibility and performance of HRES in various forms, e.g., Nehrir et al. [22], Celik [4], Kolhe et al. [17], Papadopoulos and Karagiannidis [24], Rentizelas et al. [25] and Xydis [29]. In particular, Nehrir et al. [22] developed a computer-based approach for evaluating the general performance of stand-alone PV/wind generating systems. Celik [4] proposed a technique to analyze the performance of the autonomous small-scale PV and wind hybrid energy systems. Kolhe et al. [17] presented an analysis of the performance of a hybrid PV/wind energy system with hydrogen energy storage for long-term utilization. Papadopoulos and Karagiannidis [24] presented a study on determining the achievable penetration of renewable energy sources into an insular system for the purpose of electricity generation. Rentizelas et al. [25] investigated the effect of various scenarios for emission allowance price and provided direction as to which technologies are most probable to dominate the market in the future. Xydis [29] studied the wind potential of Kythira Island and provided a techno-economic analysis aimed at identifying the optimal solution for the proposed Wind Farms (WF) to be installed so that this isolated island could be interconnected to the mainland. Some research has been focused on determining the optimal size of a hybrid PV/wind energy system. In particular, Ai et al. [1] presented a method to calculate the optimal size of a hybrid PV/wind energy system in which performance was compared on an hourly basis. The optimization approaches used to design the hybrid PV/wind energy system in the most cost effective way include linear programming, e.g., Wene and Rydén [28], Kellogg et al. [16], Chedid and Saliba [8], Ferrer-Martí et al. [14], probabilistic analysis, e.g., Karaki et al. [15], Bagul et al. [3], iterative technique, e.g, Kellogg et al. [16], dynamic programming, e.g., Musgrove [20], and multi-objective method, e.g., Yokoyama et al. [30]. In particular, Chedid and Rahman [7] used a linear programming approach to design a hybrid PV/wind power system for autonomous or gridlinked applications. Elhadidy and Shaahid [13] and Shaahid and Elhadidy [26] detailed a variety of aspects of PV, wind, and battery-based hybrid systems including optimal sizing and operation. Karaki et al. [15] presented a probabilistic assessment for an autonomous PV–wind energy conversion system (SWECS) composed of several wind turbines (wind farm), several PV modules (PV park), and a battery storage feeding a load. Furthermore, Makkonen and Lahdelma [19] formulated the decision problem of a power company as a large mixed integer programming (MIP) model and developed a customized Branch-and-Bound algorithm for solving the sub-problems efficiently. Ferrer-Martí et al. [14] proposed a mathematical programming model to optimize the design of hybrid wind–PV systems. Gomez et al. [11] proposed a new model for characterizing the energetic behavior of grid connected PV inverters. Chang [5] proposed a decision support system that integrates a proposed model and an efficient solution method to support the planning and coordination of HRES. Maheri [18] developed a robust method for obtaining the optimal design of wind– PV–diesel hybrid systems. Cicirello and Langley [9] developed efficient algorithms for propagating parametric uncertainty using the hybrid Finite Element/Statistical Energy Analysis (FE/SEA) approach to the analysis of complex vibro-acoustic systems. Alsayed et al. [2] used a suitable procedure based on a Multi Criteria Decision Analysis (MCDA) optimization approach to achieve the optimal design of grid connected PV–WT PGSs. In this research, we investigate the use of Monte Carlo simulation, along with simulation optimization techniques, for obtaining the optimal design of HRES in uncertain environment. As shown in Fig. 1, we consider a region consisting of several demand areas and several power stations and all the power stations jointly supply power for all demand areas. Each power station includes PV, wind and diesel power generators as well as the energy storage systems, and there is an electricity grid that connects all power stations and demand areas. The model is to seek the optimal size of PV, wind, and diesel power generators and the optimal size of the energy storage systems so as to achieve minimum expected total cost, while satisfying the power demand of each area. We propose a simulation optimization model to characterize the problem. Compared to the model in Chang [5] where only the size of PV and wind power generators is considered, the proposed model is more elaborate. This is not only because it further takes into consideration the diesel power generator, but also it integrates the decisions on the size of PV, wind and diesel power generators as well as on the size of the energy storage system in each power station. Other special characteristics of the proposed model distinct from the literature are described as follows: first, the renewable energy supply and the power demand are treated as continuous random variables chosen based on data collected in the real world. This contrasts with the models in the literature where uncertainties are represented by a few scenarios associated with a given probability that reflects how likely is the scenario to happen. In practice, however, a limited number of scenarios may not be sufficient to fully represent the outcome of the uncertainty. Second, the number of decision variables and the number of uncertain parameters pertaining to the problem are large due to more practical considerations. When further taking into account the complex interactions between decision variables and uncertain parameters, the traditional approaches are not applicable and the only feasible way is to resort to simulation. Third, the maintenance and operation cost as well as the carbon emission cost are taken into consideration in the objective function used for determining the optimal design of HRES. Details about the proposed model and the solution method will be presented in later sections. The remainder of this article is organized as follows. In Section 2, we present the proposed model that characterizes the decisions on the equipment installation of the HRES in uncertain environments. In Section 3, we present the metamodelbased solution method that allows us to obtain the optimal design of HRES efficiently. In Section 4, we conduct an extensive computation study to understand of performance of the proposed solution method when solving problems of realistic size. We also compare the performance of the proposed solution method with that of other two existing algorithms. We conclude with future research in Section 5.

42

K.-H. Chang, G. Lin / Simulation Modelling Practice and Theory 52 (2015) 40–51

2. The proposed model In this section, we develop a model to characterize the decisions on equipment installation of the PV–wind–diesel hybrid systems and the energy storage systems in uncertain environments. Specifically, the model aims to find the optimal size of PV, wind and diesel power generators and the optimal size of the energy storage systems, respectively, in each power station so as to achieve the minimum expected total cost, while satisfying the power demand of each area. The total cost considered in the model consists of (1) the equipment installation cost of PV, wind and diesel power generators as well as the energy storage systems in each power station (henceforth called ‘‘equipment installation cost’’); (2) the maintenance and operation cost of PV, wind, and diesel power generators and energy storage systems (henceforth called ‘‘maintenance cost’’); (3) the penalty cost incurred when the total power demand is greater than the total power supply (henceforth called ‘‘power shortage cost’’); (4) the power generation cost for using the diesel power generators to reduce the difference between the power demand and supply (henceforth called ‘‘power generation cost’’); (5) the energy storage cost of storing excess power when the energy supply is greater than the total power demand (henceforth called ‘‘energy storage cost’’); and (6) the cost of carbon emission associated with three different types of power generators, including PV, wind and diesel power generators (henceforth called ‘‘carbon emission cost’’). Among these costs, the equipment installation cost is a one-time charge, whereas other costs are incurred at each time period and will be accumulated till the end of the planning horizon. Also, only the equipment installation cost is deterministic; other costs all involve uncertain parameters and thus are stochastic. To present the model, suppose a region has M power stations and N demand areas. Let the equipment installation cost of one unit PV, wind and diesel power generators be C s ; C w and C d respectively. Moreover, let C e represent the equipment installation cost of one unit energy storage systems. Regardless of PV, wind or diesel power generators, we assume the power generation will emit carbon dioxide. The unit cost of carbon emission of PV, wind and diesel power generators are Ds ; Dw and Dd , respectively. To be realistic, it is assumed there is power loss when transmitting power among stations or from power stations to demand areas. Let 0 < gij ; gik < 1 denote the power transmission loss between the power stations i; j and from the power station i to the demand area k, respectively. The decision variables are x ¼ ½xs1 ; xs2 ; . . . ; xsM ; xw1 ; xw2 ; . . . ; xwM ; xd1 ; xd2 ; . . . ; xdM  and s ¼ ½s1 ; s2 ; . . . ; sM , where xsi ; xwi and xdi correspond to the size of PV, wind and diesel power generators installed in each power station and si corresponds to the size of energy storage systems selected for each power station. Assume the planning horizon is L time periods. Let ms ; mw ; md be the maintenance cost of one unit of PV, wind, diesel power generators, respectively and ms be the maintenance cost of one unit of energy storage systems. Other notations associated with the mathematical model are given as follows (see Table 1). The unit for Eil ; P si;l ; P wi;l ; bil ; dkl can be, for example, MW h/month, while the unit for h; zij ; yik ; eI il ; fkl can be MW h. Due to the uncertainty arising from the renewable energy supply, the power demand and the power transmission loss, Psi;l ; Pwi;l ; dkl ; gij ; gik are modeled as random variables in the proposed model. With the introduction of the notations above, the proposed model can be presented as follows:

Fig. 1. An illustration of the energy transmission network Chang [5].

K.-H. Chang, G. Lin / Simulation Modelling Practice and Theory 52 (2015) 40–51

43

Table 1 Notations used in the proposed model. Eil Ci h zij r yik P si;l P wi;l bil dkl eI il

The amount of power generated by diesel power generators at power station i in period l The unit cost of diesel power generation at power station i the amount of power that one unit energy storage system can store at power station i The amount of power transmitted from power station i to power station j The energy storage cost per unit power per period The amount of energy transmitted from station i to demand area k The amount of power generated by PV generators at power station i in period l The amount of power generated by wind generators at power station i in period l The amount of power generated by the renewable energy sources at power station i in period l The power demand of area k in period l The amount of energy stored at power station i in period l

fkl

The energy shortage cost at area k in period l

Problem 1.

Minx;s

qðx; sÞ ¼ E½Q ðx; s; xÞ

xsi 2 ½0; xsi  # Z þ ; xwi 2 ½0; xwi  # Z þ ; xdi 2 ½0; xdi  # Z þ si 2 ½0; si  # Z þ

i ¼ 1; . . . ; M

i ¼ 1; . . . ; M

ð1Þ ð2Þ

where Q ðx; s; xÞ represents the total cost, which depends on the decision variables x; s and the underlying randomness x ¼ ½Psi;l ; Pwi;l ; dkl ; gij ; gik f16i;j6M;16k6N;16l6Lg . Eqs. (1) and (2) state that the size of the PV and wind power generators, and the

size of the energy storage systems of each power station have upper bounds and are all integer-valued. The total cost is a function of the random vector x and thus it is random as well. Observe that the expectation of the objective function is taken with respect to all the random variables existing in the model. Let lðx; sÞ; aðx; sÞ; bðx; s; xÞ; kðx; s; xÞ; jðx; s; xÞ, and cðx; s; xÞ represent the equipment installation cost, the maintenance cost, the power shortage cost, the power generation cost, the energy storage cost, and the carbon emission cost, respectively. The objective function can hence be expressed as

E½Q ðx; s; xÞ ¼ lðx; sÞ þ aðx; sÞ þ E½bðx; s; xÞ þ E½kðx; s; xÞ þ E½jðx; s; xÞ þ E½cðx; s; xÞ:

ð3Þ

In particular, the equipment installation cost

lðx; sÞ ¼ C s

M X

xsi þ C w

i¼1

M M M X X X xwi þ C d xdi þ C e si ; i¼1

i¼1

i¼1

which is deterministic and is only dependent upon the decision variables x and s. Similarly, the maintenance cost over the whole planning horizon L can be expressed as

aðx; sÞ ¼ ms L

M M M M X X X X xsi þ mw L xwi þ md L xdi þ me L si : i¼1

i¼1

i¼1

i¼1

The remaining costs, including bðx; s; xÞ; kðx; s; xÞ; jðx; s; xÞ, and cðx; s; xÞ, are dependent upon not only the decision variables x; s, but also the realization of x. Due to the complex, nonlinear and unknown interaction between x; s and x, these expectations cannot be analytically expressed but can be estimated through fast-running simulation. The determination of the size of PV and wind power generators is actually a non-trivial trade-off. Specifically, when the size of PV and wind generators is selected large, it requires a high equipment installation cost at first. However, because sufficient renewable power will be likely to be supplied, a less power generation cost is incurred later. A similar argument also applies to the choice of the size of diesel power generators and the size of energy storage systems. For example, the selection of a large size of diesel power generators would incur a larger equipment installation cost. However, the likelihood that the power shortage can occur is decreased and so is the power shortage cost. To obtain estimates of the power shortage cost, the energy storage cost, the power generation cost and the carbon emission cost, we present the following power transmission model: Power transmission model (for any period l):

Minbil ;zij ;yik s:t:

M X

M X

C i Eil

i¼1

gik yik P dkl ; for k ¼ 1; . . . ; N

i¼1

eI il1 þ bil þ

X j–i

gji zji þ Eil P

ð4Þ

N X

X

k¼1

j–i

gik yik þ

gij zij for i ¼ 1; . . . ; M

bil ¼ Psi;l xsi þ P wi;l xwi for i ¼ 1; . . . ; M yik P 0; zij P 0; Eil P 0 for i; j ¼ 1; . . . ; M and k ¼ 1; . . . ; N

ð5Þ ð6Þ

44

K.-H. Chang, G. Lin / Simulation Modelling Practice and Theory 52 (2015) 40–51

In the power transmission model, Eq. (4) ensures that the power supply should satisfy the power demand of each area, after taking the power transmission loss into consideration. Eqs. (5) and (6) require that the total energy supply of each power station i, including the power generated by the PV, wind and diesel power generators, the power obtained from other power stations, and the power that is stored previously, should be at least equal to the total transmitted to other power stations or demand areas. Let the maximal amount of power that one unit diesel power generator can generate be C a . Then, in the period l the amount of power that the power station i would generate from the diesel power generator is minfEil ; C a xdi g. Based on the proposed power transmission model, for a fixed pair of ðx; sÞ, one realization of the random variables x P   yields one optimal solution bil ; zij ; yik with the corresponding objective value M i¼1 CiEil . Consequently, the amount of power that can be stored in the end of period l in the power station i can be calculated by

( ( !)) N X X X       eI il ¼ min hsi ; max 0; b þ e gji zji þ I il1 þ min Eil ; C a xdi  yik þ zij : i j–i

k¼1

j–i

One observation for the power shortage cost, the energy storage cost, the power generation cost and the carbon emission cost at the power station i over the whole planning horizon L can be obtained as follows:

bðx; s; xÞ ¼

( ) L X M N X X fil max 0; dkl  eI il1  minfEil ; C a xi g l¼1 i¼1

kðx; s; xÞ ¼

ð7Þ

k¼1

L X M X C i minfEil ; C a xi g

ð8Þ

l¼1 i¼1

jðx; s; xÞ ¼

L X M X reI il

ð9Þ

l¼1 i¼1

cðx; s; xÞ ¼

L X M X Ds  Psi;l þ Dw  P wi;l þ Dc minfEil ; C a xdi g:

ð10Þ

l¼1 i¼1

In Section 3, we will present the solution method that utilizes a Monte Carlo simulation approach to estimate the objective function, followed by the application of a simulation optimization algorithm to solve the proposed model. 3. Solution method 3.1. Exploiting the structure of the proposed model We first exploit the structure of the proposed model and then develop a solution method to derive the optimal size of the PV, wind and diesel power generators and the optimal size of the energy storage systems in each power station so as to achieve the minimum expected total cost. A detailed examination of the proposed model indicates that the solution space is quite large when the number of power stations is large. For example, when M ¼ 10, there will be 40 decision variables. If each decision variable has 50 possible values, there will be a total of 4050 alternative configurations for the HRES design. This makes it impossible to exhaustively evaluate the total cost of each configuration and select the one with the minimum cost. One way to handle this problem is to employ efficient experimental designs and develop a metamodel (also called a local model) to characterize the behavior of the objective function in a local region. By exploring the local region and optimizing the metamodel in succession, the optimal solution can be located more efficiently. Chang [5] employs a metamodel-based method, called metamodel-based simulated annealing (MSA) to solve the problem. The main framework of MSA is like that of simulated annealing [27] but further includes effective modifications to improve the computational efficiency—a new solution is selected based on the constructed metamodels. In this research, we propose a solution method, called A-STRONG, which is based on the stochastic trust-region responsesurface method (STRONG) [6]—an improved framework based on response-surface-methodology (RSM)—but with appropriate modifications to enable the solving of the proposed model. As shown in Fig. 2, RSM is one of the most widely-used metamodelling-based methodology that sequentially approximates the objective function by employing efficient experimental designs and solving for the optimal solution based on the metamodels [21]. The advantages of STRONG are that, like RSM, it incorporates the well-established statistical tools, including design of experiment (DOE) and regression analysis, in its framework and hence has better computational efficiency than other approaches when solving higher-dimensional problems. The integration of the concept ‘‘trust region,’’ originally developed in deterministic nonlinear programming [10], further enables STRONG to be automated and possess the desirable convergence guarantee, i.e., the algorithm is guaranteed to converge to the true optimum with probability one (w.p.1). While STRONG is efficient for high-dimensional problems, some modifications are required before applying it to solve the proposed model. First, STRONG assumes a credible simulation model ready for output analysis. The estimate of objective value of each configuration is available whenever needed. Second, STRONG is focused on continuous, rather than integer-valued, decision variables. To overcome these two problems, we first propose a Monte Carlo simulation approach to obtain estimates of the objective value. Then, we utilize the STRONG algorithm, coupled with the Monte Carlo simulation approach and

K.-H. Chang, G. Lin / Simulation Modelling Practice and Theory 52 (2015) 40–51

45

Fig. 2. An illustration for RSM.

appropriate modifications, to enable the solving of the proposed model. The adapted algorithm is called A-STRONG. The main framework of A-STRONG is to treat the problem as a continuous one first and utilize STRONG to find a nearly-optimal solution. Then, the integer-valued constraints are taken into consideration by searching the ‘‘neighborhood’’ of the nearly-optimal solution. More details will be presented later. 3.2. Estimation of energy storage cost, diesel power generation cost and carbon emission cost Suppose the sample size is N. We propose to estimate the objective function of the proposed model as follows. As shown in Fig. 3, for one given solution x; s; N’s x is first generated. Specifically, each element of x is generated according to their respective distribution. Then, the power coordination model is solved N times. By using Eqs. (8)–(10), N observations can be obtained and the sample averages bðx; s; xÞ, kðx; s; xÞ; jðx; s; xÞ, and cðx; s; xÞ, given as Eqs. (11)–(14), serve as the estimate of E½bðx; s; xÞ, E½kðx; s; xÞ; E½jðx; s; xÞ, and E½cðx; s; xÞ, respectively.

PN

xi Þ N PN kðx; s; xi Þ kðx; s; xÞ ¼ i¼1 N PN i¼1 jðx; s; xi Þ jðx; s; xÞ ¼ N PN c ðx; s; xi Þ cðx; s; xÞ ¼ i¼1 N bðx; s; xÞ ¼

i¼1 bðx; s;

ð11Þ ð12Þ ð13Þ ð14Þ

The total cost is then estimated by

E½Q ðx; s; xÞ  lðx; sÞ þ aðx; sÞ þ bðx; s; xÞ þ kðx; s; xÞ þ jðx; s; xÞ þ cðx; s; xÞ: 3.3. An overview of STRONG STRONG is a metamodel-based framework. For each iteration, STRONG employs an experimental design and builds a metamodel to characterize the behavior of the objective function in a local region. As soon as the metamodel is built, the

Fig. 3. Estimating the objective function by Monte Carlo simulation.

46

K.-H. Chang, G. Lin / Simulation Modelling Practice and Theory 52 (2015) 40–51

optimal solution of the metamodel, subject to the ‘‘trust region,’’ is solved. The term ‘‘trust region’’ is analogous to that in trust region method in nonlinear programming [10]. It refers to a subregion in the solution space where the local model is believed to well represent the objective. The size of trust region is updated by the algorithm according to the quality of the new solution. If the new solution is satisfactory, it is accepted and the trust region either remains the same size or expands to allow for a larger moving step in the next iteration. Otherwise, the new solution is rejected and the trust region shrinks to restrict the moving step in the next iteration. The main framework of STRONG is illustrated in Fig. 4. In words, depending on the size of the trust region, STRONG has two stages (Stages I and II), which we call outer loop, and an inner loop. If the size of the trust region is large, it is suggested that the previous iteration(s) were successful, so a first-order model should be able to represent the objective function. Thus, STRONG constructs and optimizes the first-order model (Stage I). On the other hand, if the size of trust region is small, it is implied that the previous iteration(s) cannot find satisfactory solutions, and a more complicated model, i.e., a second-order model, is necessary. Thus, STRONG constructs and optimizes the second-order model (Stage II). If STRONG still cannot find a satisfactory solution in Stage II, the inner loop is initiated, where some effective mechanisms are employed to ensure that the quality of the local model is continuously improved and finally is able to result in a satisfactory solution. In what follows, we provide a schematic representation of Stages I, II and the inner loop, while leaving the technical details in Appendix. A full presentation about STRONG is given in Chang et al. [6]. 3.4. Stages I and II Suppose the problem dimension is m and the solution at iteration k is tk ¼ ðxk ; sk Þ and the radius of the trust region is Dk > 0, which is automatically determined by the algorithm. The trust region is defined as:

Bðtk ; Dk Þ ¼ ft 2 Rm : kt  tk k 6 Dk g;

ð15Þ

where k  k denotes the Euclidian norm.

Fig. 4. The main framework of STRONG.

K.-H. Chang, G. Lin / Simulation Modelling Practice and Theory 52 (2015) 40–51

47

For any iteration k, the schematic steps of Stages I and II are given as follows. Step 1. Construct Resolution III fractional factorial designs (Stage I) or central composite designs (Stage II) around the center point tk . Step 2. Build a first-order model (Stage I) or second-order model (Stage II) r k ðtÞ (Appendix A.1). Step 3. Solve tk 2 argminfrk ðtÞ : t 2 Bðtk ; Dk Þg (Appendix A.2). Step 4. Simulate observations at tk and tk , and estimate qðtk Þ. Conduct ratio-comparison (RC) and sufficient-reductions (SR) tests to examine the quality of xk and update the size of trust region based on the result of the RC test. If tk passes both the RC and SR tests, let tkþ1 ¼ tk . Otherwise, go back to the main framework if currently in Stage I or go to the inner loop if currently in Stage II (Appendix A.3). Note that the differences between Stages I and II is that Stage I employs fractional factorial designs and builds a first-order model. When the algorithm switches to Stage II, it is implied that Stage I cannot find a satisfactory solution based on firstorder models, Stage II therefore employs central composite designs and builds a more elaborate model, i.e., second-order model. Without causing confusion, we suppress x and let Q 1 ðtÞ; . . . ; Q n ðtÞ denote independent and identically distributed (i.i.d.) random variables of Q ðtÞ. Further let the average

Q ðt; nÞ ¼

n 1X Q ðtÞ: n i¼1 i

ð16Þ

q k ðtk Þ and b q k ðtk Þ, by Q ðtk ; N k Þ and STRONG estimates the objective value of the current solution and the new solution, b  Q ðtk ; N k Þ, respectively. The RC test, defined by

qk ¼

b q k ðtk Þ  b q k ðtk Þ ; r k ðtk Þ  r k ðtk Þ

ð17Þ

compares the ‘‘observed reduction’’ estimated based on the difference of sample averages of the two solutions and the ‘‘predicted reduction’’ estimated based on the difference of values of the built local model of the two solutions. If the ratio is large, i.e., qk P g1 , which implies that the new solution can yield significant reductions on function value, the new solution is accepted and the trust region is enlarged. On the other hand, if the ratio is moderate (g0 6 qk < g1 ), which implies that the ‘‘observed reduction’’ has fair agreements with the ‘‘predicted reduction,’’ the new solution is also accepted and the trust region should remain the same size. Lastly, if qk is close to 0 or negative (qk < g0 ), which implies that the new solution is poor, the new solution should be rejected and the trust region is shrunk. Here g0 and g1 satisfying 0 < g0 < g1 < 1 are determined by users; for example, we can let g0 ¼ 1=4 and g1 ¼ 3=4. Since the RC test contains sampling error, when a new solution passes the RC test, the SR test is further employed to ascertain that the new solution can yield ‘‘sufficient reductions’’ in the objective function. The SR test is defined as follows:

H0 : qðtk Þ  qðtk Þ 6 g20 fk

vs: H1 : qðtk Þ  qðtk Þ > g20 fk

where g20 fk is so-called ‘‘sufficient reductions,’’ which is defined in Lemma 1 of Chang et al. [6]. When H0 is rejected, it is concluded that the new solution can yield ‘‘sufficient reductions.’’ In this case, the solution is considered a satisfactory solution. In order to avoid that a poor solution is (incorrectly) accepted, the type-I error of the SR test, ak , which refers to the probability that a poor solution is (incorrectly) accepted, has to be controlled. For STRONG to achieve convergence, it is P required that 1 k¼1 ak < 1. 3.5. Inner loop When Stage II fails to find a satisfactory solution, the inner loop is initiated. For each loop, more design points are employed and accumulated to construct/improve the metamodel. The sample size of the current and new solutions is also increased to enhance the correctness of the result of the RC and SR tests. Moreover, the trust region is continued to shrink, restricting the optimization step to a smaller region. In what follows, we present the steps of the inner loop. Before doing so, we first introduce a few notations. Let ki be the ith inner loop of the kth iteration and Dki be the trust region in ki . Moreover, tki 2 Rm denotes the new solution found in ki . For any inner loop i, the steps taken in the inner loop are given as follows. Step 1. Construct/accumulate central composite designs (CCDs) around the center point tk . Increase the sample size of the current and new solutions. Step 2. Build/improve the second-order model r ki ðtÞ. Step 3. Solve tki 2 argminfrki ðtÞ : t 2 Bðtk ; Dki Þg. Step 4. Simulate observations at tk and tki and estimate qðtk Þ and qðtki Þ.

48

K.-H. Chang, G. Lin / Simulation Modelling Practice and Theory 52 (2015) 40–51

Step 5. Conduct RC and SR tests to examine the quality of tki ; update the size of the trust region. If tki passes both the RC and SR tests, let tkþ1 ¼ tki and go back to the main framework. Otherwise, go back to Step 1 of the inner loop. It can be proved that as the number of inner loop iterations increases the estimated gradient of the metamodel converges to the true gradient of the underlying response surface and the ‘‘observed reduction’’ converges to the difference of true function values. These altogether ensure the correctness of the RC test and reduce the probability of (incorrectly) accepting a poor solution and (incorrectly) updating the size of the trust region. In Lemma 4 of Chang et al. [6], it is shown that if the algorithm is not at a stationary point, i.e., krqðtk Þk > 0, the algorithm can always find a new satisfactory solution, i.e., a solution that passes both the RC and SR tests. 3.6. Determining the optimal solution of the proposed model By using STRONG, the nearly-optimal solution of the proposed model can be found. In this section, we discuss how to adapt STRONG to determine the optimal size of the PV and wind power generators and the optimal size of the energy storage systems. The basic idea is to search for the optimal integer solution in the neighborhood  of the nearly-optimal solution obtained through STRONG. Specifically, suppose tk ¼ xks1 ; . . . ; xksM ; xkw1 ; . . . ; xkwM ; sk1 ; . . . ; skM is a nearly-optimal solution. Further suppose the final metamodel before the algorithm stops is r ki . The optimal solution of the proposed model is determined by solving the following problem: Problem 2.

min r ki ðx; sÞ   for i ¼ 1; . . . ; M: subject to xsi 2 dxksi e; bxksi c  k  k xwi 2 dxwi e; bxwi c for i ¼ 1; . . . ; M:   for i ¼ 1; . . . ; M: xdi 2 dxkdi e; bxkdi c   si 2 dski e; bski c : for i ¼ 1; . . . ; M: Problem 2 seeks the integer-valued optimal size of the PV and wind power generators and the optimal size of the energy storage systems that require the minimal cost within the neighborhood of the nearly-optimal solution tk .

4. Numerical experiments In order to understand the computational performance of the adapted STRONG algorithm (A-STRONG) relative to other existing algorithms, we create six scenarios based on the proposed model and use three algorithms to solve them. The algorithms under comparison are A-STRONG, simulated annealing (SAN) [27] and Nelder-Mead simplex method (NM) [23]. In particular, SAN is a random-search-based algorithm developed based on the observation of the cooling process of a liquid or solid in which molecules have higher mobility at high temperatures and the mobility is lost or declines as the temperature decreases. This concept has been successfully applied to solve some optimization problems, finding the lowest (or highest) value of the objective function. On the other hand, NM, originally developed in deterministic optimization, is an efficient direct search method that optimizes the response function merely by comparing function values. The six scenarios are constituted by three kinds of problem sizes and two kinds of uncertainty intensities. Specifically, the three kinds of problem sizes include small ðM ¼ 3; N ¼ 10Þ, medium ðM ¼ 6; N ¼ 30Þ and large number of power stations and demand areas ðM ¼ 20; N ¼ 100Þ. The two kinds of uncertainty intensities correspond to the settings of small and large variances for the power demand and the renewable energy supply: dkl  Normal ð100; 52 Þ; Psi;l ; P wi;l  Normal ð5; 12 Þ, and dkl  Normal ð100; 502 Þ; P si;l ; Pwi;l  Normal ð5; 32 Þ. In particular, the power transmission loss gij is set as a beta random variable, whose mean and variance depend on the distance between i and j, and 0 6 gij 6 1. Detailed settings associated with each scenario are summarized in Table 2.

Table 2 The problem setting associated with six scenarios. Problem settings Scenario 1

M ¼ 3; N ¼ 10; dkl  Normalð100; 52 Þ; P si;l and P wi;l  Normalð5; 12 Þ

Scenario 2

M ¼ 3; N ¼ 10; dkl  Normalð100; 502 Þ; P si;l and P wi;l  Normalð5; 32 Þ

Scenario 3

M ¼ 6; N ¼ 30; dkl  Normalð100; 52 Þ; P si;l and P wi;l  Normalð5; 12 Þ

Scenario 4

M ¼ 6; N ¼ 30; dkl  Normalð100; 502 Þ; P si;l and P wi;l  Normalð5; 32 Þ

Scenario 5

M ¼ 20; N ¼ 100; dkl  Normalð100; 52 Þ; P si;l and P wi;l  Normalð5; 12 Þ

Scenario 6

M ¼ 20; N ¼ 100; dkl  Normalð100; 502 Þ; P si;l and P wi;l  Normalð5; 32 Þ

49

K.-H. Chang, G. Lin / Simulation Modelling Practice and Theory 52 (2015) 40–51 Table 3 Parameter settings for all scenarios. Parameters

n0

nd

D0

~ D

g0

g1

c1

c2

ak

Settings

P3

P2

2

1.2

0.01

0.3

0.9

1.11

0:5  0:98k

Table 4 Performance of SAN, NM and STRONG in 6 scenarios.

Scenario Scenario Scenario Scenario Scenario Scenario

1 2 3 4 5 6

SAN

NM

A-STRONG

2192 (1057) 2957 (934) 10,270 (3250) 24,037 (12,229) 50,270 (12,574) 81,037 (29,305)

1823 (812) 2824 (1135) 8723 (3230) 19,831 (3212) 32,419 (11323) 63,416 (21345)

1391 (703) 2050 (1012) 6730 (2140) 14,921 (8731) 23,174 (9075) 21,042 (12035)

14000

A−STRONG NM SAN

12000

Total cost

10000 8000 6000 4000 2000 0

0

1000

2000

3000

4000

5000

6000

7000

8000

Number of observations Fig. 5. The convergence behaviors of A-STRONG, SAN and NM with Scenario 1.

Other parameter settings are introduced as follows: xsi ; xwi ¼ 30; xdi ¼ 10; si ¼ 15; C i ¼ 1 for all i; C s ¼ 5; C w ¼ 10; C d ¼ 10; C e ¼ 5, ms ¼ 1; mw ¼ md ¼ 3; me ¼ 5; C a ¼ 1, Ds ¼ Dw ¼ 1 and Dc ¼ 10. The parameter setting of the algorithm is given in Table 3. The final optimal solution is returned with the corresponding objective value whenever the algorithm has consumed 5000, 10,000 and 100,000 observations for the small, medium and large problems. All algorithms are run for 20 times for each scenario. The average total costs over the 20 macroreplications and the associated standard deviation (shown in the parenthesis) are given in Table 4. In Scenario 1, for example, it takes around 165 s to solve the problem on the computer with Intel(R) Core(TM)2 Duo CPU (2.26 GHz, 2.27 GHz) and 4 GB of RAM. Based on Table 4, it can be found that, given the same computational budget, A-STRONG can find the optimal solution with much better quality (i.e., the lower total cost) than that of SAN and NM for six scenarios. Specifically, the average total cost of the optimal solution found by A-STRONG is around 30–50% lower than that of SAN and NM depending on the scenario. Moreover, the computational gain of A-STRONG increases with the problem size and the uncertainty intensity. This is because both SAN and NM are random-search-based algorithms, which only explore one point at each iteration. By contrast, A-STRONG explores a ‘‘region’’ at each iteration and hence is much more efficient. To enhance the understanding of the convergence behavior of A-STRONG, SAN and NM, we report one sample path of the three algorithms with Scenario 1 in Fig. 5. It can be seen that A-STRONG can converge much faster than that of SAN and NM. The convergence behavior with other scenarios also exhibit the same pattern and hence are not shown here. 5. Conclusions In this paper, we propose a simulation optimization model to characterize the decision on the design of HRES in uncertain environments. The proposed model aims to seek the optimal size of PV, wind and diesel power generators as well as the

50

K.-H. Chang, G. Lin / Simulation Modelling Practice and Theory 52 (2015) 40–51

optimal size of energy storage systems, in each power station so as to achieve minimum expected total cost, while satisfying the power demand. A metamodel-based algorithm, called A-STRONG, is developed to solve the proposed model. An extensive numerical study shows that A-STRONG can solve the proposed model of realistic size efficiently. Moreover, because the search strategy of A-STRONG is to explore a region, instead of a point, at one time, its computational efficiency is superior to that of the existing algorithms, namely, SAN and NM. Two directions can be explored in future research. First, when there is limited historical data for estimating the distribution of the renewable energy supply and the power demand, the proposed model should account for the distribution ambiguities. A solution method that can produce a robust solution regardless of the distribution ambiguity should be developed. Second, the computational efficiency of the proposed algorithm may be further enhanced by applying useful variance reduction techniques, such as common random numbers. Acknowledgments The authors would like to thank two anonymous referees for their insightful comments and suggestions that have significantly improved this paper. This research is conducted under the ‘‘III Innovative and Prospective Technologies Project’’ subsidized by the Ministry of Economy Affairs of Taiwan. Appendix A A.1. Metamodel construction For the outer loop, suppose x is the center point, and l is number of design points. Let X l denote the design matrix with e l the design matrix with the main effect, quadratic and interaction columns. We use X l to only the main effect columns, and X e l to construct second-order models. construct first-order models, while using X For each iteration, STRONG constructs a first- or second-order model with the selected factors,

b qT ðtk Þðt  tk Þ and q k ðtk Þ þ r r k ðtÞ ¼ b k 1 b b k ðtk Þðt  tk Þ; b r k ðtÞ ¼ q k ðtk Þ þ r qTk ðtk Þðt  tk Þ þ ðt  tk ÞT H 2 respectively. Let

Y m ¼ ðy1 ; . . . ; ym ÞT  Qðt; nÞ be the centralized response vector, where yj is the output of the simulation experiment. The gradient rqðtÞ and Hessian matrix HðtÞ are estimated as follows:



b qðtÞ ¼ X T X l r l

"

b qðtÞ r b HðtÞ

#

1

X Tl Y l

ðestimating only the gradientÞ;

 1 eT X el eT Yl ¼ X X l l

ðestimating the gradient and Hessian matrixÞ;

ð18Þ ð19Þ

where the gradient and Hessian are estimated by applying the ordinary least squares (OLS) method. For the inner loop, STRONG accumulates the design points to estimate the gradient and Hessian of second-order models. Specifically, for the inner loop ki , suppose U lk is the newly-generated main-effect orthogonal design within the current trust i h iT includes all the design points up to the inner loop ki . region Dki , and the design matrix X lk ¼ X lk ; U lk i

i1

i

A.2. Solving the subproblem The subproblem is defined as ~tk 2 argminfr k ðtÞ : t 2 Bð~tk ; Dk Þg and can be solved by the following steps: n o b qT ðtk Þd : kdk 6 Dk . q k ðtk Þ þ r Step 1. Find the steepest descent direction dk ¼ arg min b k Step 2. Choose a step size sk ¼ arg min fr k ðsdk Þ : s > 0; ksdk k 6 Dk g. Step 3. Let tk ¼ tk þ sk dk . The similar idea applies to inner loop except k is replaced by ki . A.3. Updating the size of the trust region STRONG-X automatically updates the size of the trust region based on the results of the RC and SR tests. Let 0 < c1 < 1 < c2 , for example, c1 ¼ 1=2 and c2 ¼ 2. In Stage I, if tk fails in either the RC or SR test, the center point remains

K.-H. Chang, G. Lin / Simulation Modelling Practice and Theory 52 (2015) 40–51

51

unchanged and the trust region shrinks, i.e., tkþ1 ¼ tk and Dkþ1 ¼ c1 Dk ; if qk P g1 and tk passes the SR test, the center point then moves to the new solution and the trust region enlarges, i.e., tkþ1 ¼ tk and Dkþ1 ¼ c2 Dk . If g0 6 qk < g1 and tk passes the SR test, then the center point will move to the new solution and the trust region remains the same size; i.e., tkþ1 ¼ tk and Dkþ1 ¼ Dk . References [1] B. Ai, H. Yang, H. Shen, X. Liao, Computer-aided design of PV/wind hybrid system, Renew. Energy 28 (10) (2003) 1491–1512. [2] M. Alsayed, M. Cacciato, G. Scarcella, G. Scelba, Design of hybrid power generation systems based on multi criteria decision analysis, Sol. Energy 105 (2014) 548–560. [3] A.D. Bagul, Z.M. Salameh, B. Borowy, Sizing of stand-alone hybrid wind–photovoltaic system using a three-event probability density approximation, Sol. Energy 56 (4) (1996) 323–335. [4] A.N. Celik, Optimization and techno-economic analysis of autonomous photovoltaic–wind hybrid energy systems in comparison to single photovoltaic and wind systems, Energy Convers. Manage. 43 (18) (2002) 2453–2468. [5] K.-H. Chang, A decision support system for planning and coordination of hybrid renewable energy systems, Decis. Support Syst. 64 (2014) 4–13. [6] K.-H. Chang, L.J. Hong, H. Wan, Stochastic trust-region response-surface method (STRONG)-a new response-surface framework for simulation optimization, INFORMS J. Comput. 25 (2) (2013) 230–243. [7] R. Chedid, S. Rahman, Unit sizing and control of hybrid wind–solar power systems, IEEE Trans. Energy Convers. 12 (1) (1997) 79–85. [8] R. Chedid, Y. Saliba, Optimization and control of autonomous renewable energy systems, Int. J. Energy Res. 20 (1996) 609–624. [9] A. Cicirello, R.S. Langley, Efficient parametric uncertainty analysis within the hybrid finite element/statistical energy analysis method, J. Sound Vib. 333 (2014) 1698–1717. [10] A.R. Conn, N.L.M. Gould, P.L. Toint, Trust-Region Methods, SIAM, 2000. [11] L. Davila-Gomez, A. Colmenar-Santos, M. Tawfik, M. Castro-Gil, An accurate model for simulating energetic behavior of PV grid connected inverters, Simul. Modell. Pract. Theory 49 (2014) 57–72. [12] M.K. Deshmukh, S.S. Deshmukh, Modeling of hybrid renewable energy systems, Renew. Sustain. Energy Rev. 12 (1) (2008) 235–249. [13] M.A. Elhadidy, S.M. Shaahid, Optimal sizing of battery storage for hybrid (wind + diesel) power systems, Renew. Energy 18 (1) (1999) 77–86. [14] L. Ferrer-Martí, B. Domenech, A. García-Villoria, R. Pastor, A MILP model to design hybrid wind–photovoltaic isolated rural electrification projects in developing countries, Eur. J. Oper. Res. 226 (2013) 293–300. [15] S.H. Karaki, R.B. Chedid, R. Ramadan, Probabilistic performance assessment of autonomous solar–wind energy conversion systems, IEEE Trans. Energy Convers. 14 (3) (1999) 766–772. [16] W.D. Kellogg, M.H. Nehrir, G. Venkataramanan, V. Gerez, Generation unit sizing and cost analysis for stand-alone wind photovoltaic and hybrid wind/ PV systems, IEEE Trans. Energy Convers. 13 (1) (1998) 70–75. [17] M. Kolhe, K. Agbossou, J. Hamelin, T.K. Bose, Analytical model for predicting the performance of photovoltaic array coupled with a wind turbine in a stand-alone renewable energy system based on hydrogen, Renew. Energy 28 (5) (2003) 727–742. [18] A. Maheri, Multi-objective design optimisation of standalone hybrid wind–PV–diesel systems under uncertainties, Renew. Energy 66 (2014) 650–661. [19] S. Makkonen, R. Lahdelma, Non-convex power plant modelling in energy optimisation, Eur. J. Oper. Res. 171 (3) (2006) 1113–1126. [20] A.R.D. Musgrove, The optimization of hybrid energy conversion system using the dynamic programming model – Rapsody, Int. J. Energy Res. 12 (1988) 447–457. [21] R.H. Myers, D.C. Montgomery, C.M. Anderson-Cook, Response Surface Methodology-Process and Product Optimization Using Designed Experiments, John Wiley & Sons, New York, 2009. [22] M.H. Nehrir, B.J. LaMeres, G. Venkataramanan, V. Gerez, L.A. Alvarado, An approach to evaluate the general performance of stand-alone wind/ photovoltaic generating systems, IEEE Trans. Energy Convers. 15 (4) (2000) 433–439. [23] J.A. Nelder, R. Mead, A simplex method for function minimization, Comput. J. 7 (1965) 308–313. [24] A. Papadopoulos, A. Karagiannidis, Application of the multi-criteria analysis method Electre III for the optimisation of decentralised energy systems, Omega-Int. J. Manage. Sci. 36 (2008) 766–776. [25] A.A. Rentizelas, A.I. Tolis, L.P. Tatsiopoulos, Investment planning in electricity production under CO2 price uncertainty, Int. J. Prod. Econ. 140 (2) (2012) 622–629. [26] S.M. Shaahid, M.A. Elhadidy, Prospects of autonomous/stand-alone hybrid (photo-voltaic + diesel + battery) power systems in commercial applications in hot regions, Renew. Energy 29 (2) (2003) 165–177. [27] J.C. Spall, Introduction to Stochastic Search and Optimization: Estimation, Simulation, and Control, John Wiley & Sons, New York, 2003. [28] C.-O. Wene, B. Rydén, A comprehensive energy model in the municipal energy planning process, Eur. J. Oper. Res. 33 (2) (1988) 212–222. [29] G. Xydis, A techno-economic and spatial analysis for the optimal planning of wind energy in Kythira island, Greece, Int. J. Prod. Econ. 146 (2013) 440– 452. [30] R. Yokoyama, K. Ito, Y. Yuasa, Multi-objective optimal unit sizing of hybrid power generation systems utilizing photovoltaic and wind energy, J. Sol. Energy Eng. 116 (4) (1994) 167–173.