An efficient scenario-based stochastic programming framework for multi-objective optimal micro-grid operation

Applied Energy 99 (2012) 455–470

Contents lists available at SciVerse ScienceDirect

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

An efficient scenario-based stochastic programming framework for multi-objective optimal micro-grid operation Taher Niknam ⇑, Rasoul Azizipanah-Abarghooee, Mohammad Rasoul Narimani Department of Electrical and Electronics Engineering, Shiraz University of Technology, Shiraz, Iran

h i g h l i g h t s " Proposes a stochastic model for optimal energy management. " Consider uncertainties related to the forecasted values for load demand. " Consider uncertainties of forecasted values of output power of wind and photovoltaic units. " Consider uncertainties of forecasted values of market price. " Present an improved multi-objective teaching–learning-based optimization.

a r t i c l e

i n f o

Article history: Received 27 January 2012 Received in revised form 29 February 2012 Accepted 12 April 2012 Available online 2 July 2012 Keywords: Improved teaching–learning-based algorithm Micro grid Multi-objective stochastic optimization Renewable energy management Self adaptive probabilistic modification strategy Uncertainty

a b s t r a c t This paper proposes a stochastic model for optimal energy management with the goal of cost and emission minimization. In this model, the uncertainties related to the forecasted values for load demand, available output power of wind and photovoltaic units and market price are modeled by a scenario-based stochastic programming. In the presented method, scenarios are generated by a roulette wheel mechanism based on probability distribution functions of the input random variables. Through this method, the inherent stochastic nature of the proposed problem is released and the problem is decomposed into a deterministic problem. An improved multi-objective teaching–learning-based optimization is implemented to yield the best expected Pareto optimal front. In the proposed stochastic optimization method, a novel self adaptive probabilistic modification strategy is offered to improve the performance of the presented algorithm. Also, a set of non-dominated solutions are stored in a repository during the simulation process. Meanwhile, the size of the repository is controlled by usage of a fuzzy-based clustering technique. The best expected compromise solution stored in the repository is selected via the niching mechanism in a way that solutions are encouraged to seek the lesser explored regions. The proposed framework is applied in a typical grid-connected micro grid in order to verify its efficiency and feasibility. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Deregulating in power systems has motivated the electricity companies to utilize the Distributed Generators (DGs) near the energy consumers [1,2]. DGs include different types of power sources such as diesel engines, Micro Turbines (MTs), Fuel Cells (FCs), Photo Voltaic (PV), and small Wind Turbines (WTs) [3,4]. From the consumers’ point of view, using DGs causes the decrease in electricity cost, higher service reliability, higher power quality, increase of energy efficiency, and energy independence [5]. Furthermore, using Renewable Energy Sources (RESs) like wind or photovoltaic can satisfy the entire environment concerns [6–8]. Therefore, nowadays with progress in power electronic equipments ⇑ Corresponding author. Tel.: +98 7117264121; fax: +98 7117353502. E-mail address: [email protected] (T. Niknam). 0306-2619/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.apenergy.2012.04.017

and consequently appearance of the more advanced energy storage facilities, these devices are playing an important role in low voltage power networks. They can save energy at low-price hours and sell it at high-price hours which this helps the network to work more efficiently and economically [9]. Meanwhile the operation and control of DG units in accompany with storage facilities such as flywheels, energy capacitors and batteries are getting more complicated and challengeable [10]. Above claims and important role of the Micro Grids (MGs) in modern power systems have caused many researches to be carried out in order to scrutinize different aspects of the MG problem. Chakraborty et al. [11] utilized a linear programming algorithm in order to minimize MG operation cost and to optimize the states of battery charge. Sortomme and El-Sharkawi [12] used Particle Swarm Optimization (PSO) for decreasing the costs of MGs including controllable loads and battery storage by selling stored energy

456

T. Niknam et al. / Applied Energy 99 (2012) 455–470

(a) Ordinary MCS Random Value

1 0.8 0.6 0.4 0.2 0

0

20

40

60

80

100

120

140

160

180

140

160

180

Sample No.

(b) Rank-1 LMCS Random Value

1 0.8 0.6 0.4 0.2 0

0

20

40

60

80

100

120

Sample No. Fig. 1. Random numbers generated by (a) ordinary MCS and (b) rank-1 LMCS.

at high prices. Hernandez-Aramburo et al. [13] decreased the fuel consumption rate of the system while constrained this rate to satisfy the local electrical and thermal energy demands and provided a certain minimum reserve power, this study considered operation management of MGs as a single objective optimization problem. Also, storage facilities and selling or purchasing power to/from utility were not taken into account. Mohamed and Koivo considered the presence of storage systems in [14] but they did not formulate the problem on period with several time intervals and therefore did not consider the issue of charge and discharge cycles and other pricing issues in [15]. Logenthiran and Srinivasan [16] developed the classical unit commitment problem using a three steps method for minimizing the cost objective function. Morais et al. [17] proposed the classical unit commitment. Also, all previous researches have considered the operation management of MGs as a deterministic problem. Because of uncertainty in forecasted values of available output powers of WT and PV units, market prices and load demand considering this problem as deterministic one fails to obtain reliable solutions. Dukpa et al. [18] proposed a new optimal participation strategy for a wind power generator that employs an energy storage device for participating in a day-ahead unit commitment process considering stochastic power output. The authors have also proposed a novel smart energy management system which can coordinate power forecasting, energy storage and energy exchanging and then make a proper short-term scheduling to minimize the total operation cost. The important drawback of above study is that it does not consider all uncertainties of the problem. Also, nowadays restraining emissions of greenhouse gases becomes one of the most important issues of the electricity generation sector. The MG operators need powerful energy management tools which can dispatch the power between DG sources, substation and storage facilities in order to satisfy economic and environmental goals. In this regard, solving a Multi-objective Optimization Problem (MOP) [19,20] over 24 h considering the associated uncertainties is inevitable. The novelty of the proposed approach resides in considering of the optimal scheduling of generation units accompany with the inherent uncertainties which leads to getting close to actual conditions. Furthermore, the proposed approach considers Spinning Reserve Requirements (SSRs) in order to protect system from sudden fault, which can provide high reliability for MG systems. The two-stage

scenario-based stochastic programming framework is applied to solve the proposed problem in which at first, the discrete control variables (here-and-now) are generated for all generation units and after that the continuous dispatch variables (wait-and-see) are produced by scenario based technique [21]. These scenarios are generated using the Roulette Wheel Mechanism (RWM) and Lattice Monte Carlo Simulations (LMCS) method. Moreover, for avoiding the intractable computation, a scenario reduction technique should be applied. It should be noted that there are four techniques including scenario reduction, fine tuning of the solution algorithm, decomposition, and manipulation of the model for reducing the cumbersome computation [21]. To solve the proposed problem, two aforementioned first techniques are applied in this paper to decline the computation time. Also, the Teaching–Learning-Based Optimization (TLBO) algorithm as one of the fast evolutionary algorithms, is utilized to optimize the proposed problem [22]. It should be noted that the proposed Stochastic Multi-objective Optimal MG Operation (SMOMGO) is a mixed integer nonlinear and non-differentiable optimization problem. The TLBO algorithm is a new efficient optimization algorithm which has been inspired by a learning mechanism in a class [22]. To avoid trapping in local optima, the proposed optimization algorithm is equipped with a Self Adaptive Probabilistic Modification Strategy (SAPMS) which is called Improved TLBO (ITLBO) algorithm. Finding a set of solutions in accompany with the best compromise one can help to MG Central Controller (MGCC) to have more flexibility in their choices. In this respect, a fuzzy clustering method is used to prune the size of the repository without destroying its characteristics [23]. To incorporate Decision Makers (DMs) biased preference through the search process perfectly, a min– max approach is developed to decide on the best candidate solutions for the next optimization part of the proposed algorithm. In order to encourage individuals to seek lesser explored regions with higher probability of Pareto-Optimal Front (POF), a niching technique is employed to select the expected best compromise solution [23]. The main contributions of this work can be summarized as follows: (i) Formulate the SMOMGO by consideration of cost and emission objectives simultaneously, (ii) Implement a Pareto-based approach based on a novel ITLBO for solving the proposed problem, (iii) Model and formulate the Spinning Reserve Requirements (SRRs) to overcome the sudden fault in the system, (iv) Decrease the computation time with backward scenario reduction method (v) The performance and potential of the proposed approach is successfully validated with numerical simulations. 2. Stochastic model description To solve the proposed problem, a two-stage stochastic scenariobased method is implemented in this paper. At the first stage, all discrete control variables (here-and-now) are made for all generation units before the uncertainty of input random variables is realized. The uncertain variables are the load demands, the market prices, and powers output of PV and WT units in the study horizon. These uncertainties are realized using a scenario-based technique in which scenarios are generated by RWM and LMCS. In the next stage, continuous dispatch variables (wait-and-see) are generated. The scenario generation and reduction techniques are formulated as follows. 2.1. Scenario generation Random sampling is a corner stone of the Monte Carlo method, hence using a most efficient method for scenario generation is

457

T. Niknam et al. / Applied Energy 99 (2012) 455–470

Interval 7

Interval 6

Fig. 3. Accumulated normalized probabilities of the forecast error intervals. Interval 4 Interval 6

Interval 5 Interval 7

-3 σ

-2 σ

-σ

0

σ

2σ

3σ

Load Forecast Error ( Δ PD ,ld ,t ,s ) Fig. 2. Typical discretization of the PDF of the load level forecast error.

crucial. In this regard, a random number between [0, 1] is separately generated for each input random variable. In the state of the ordinary MCS, the random numbers are uniformly distributed in the intervals, but in the LMCS, an N-point lattice rule of rank-r in d-dimension is computed as follows [24]:

! r N X kj X v i mod1; nj i¼1 j¼1

kj ¼ 1; . . . ; N j ¼ 1; . . . ; r

ð1Þ

where v1, . . . , vN are vectors with dimension d obtained by the ordinary MCS. Dimension d indicates the number of random values in each scenario. One LMCS scenario includes a vector with dimension d of random numbers in the range of [0, 1] which is constructed by means of a set of values {kj, j = 1, . . . , r}. The difference among scattering of the scenarios related to MCS and LMCS in two-dimensional space is shown in Fig. 1 which reveals the superiority of the LMCS rank-1. All uncertainties are modeled based on their error determined by the Probability Distribution Function (PDF). As an example, a typical continuous PDF accompany with its discretization is shown in Fig. 2 for the load forecast error in time period t [25]. The probability distribution of a random variable is represented by a finite set of scenarios. This procedure is common in a stochastic programming. Also, each scenario corresponds with a single realization of the random variables throughout the study horizon. In addition, each scenario has an associated probability of occurrence. Thus load demand, market price, powers output of PV and WT levels for each scenario can be presented as follows:

PWT;t;s ¼ Pforecast þ DP WT;t;s ; WT;t

Pforecast PV;t

þ DPPV;t;s ;

ð2Þ

PV ¼ 1; . . . ; NPV ; t ¼ 1; . . . ; NT ; s

¼ 1; . . . ; Ns P forecast D;ld;t

error for unit WT, photovoltaic power forecasted error for unit PV, load demand forecasted error for load level ld and price forecasted error at time t in scenario s, respectively. NWT, NPV, ND, T, and Ns are the numbers of WT and PV units, load level, time intervals and scenarios, correspondingly. Different realizations of the random variables and their probabilities are obtained by usage of the scenario-generation method described in [24]. This scenario-generation process is based on LMCS and RWM [26] summarized as follows. At first, according to the desired preciseness, the distribution function centered on the zero mean (which represents the distribution mean) is divided into some class intervals (here, seven segments). Each class interval determines one standard deviation error (r) wide [27]. Moreover, each interval l is associated with a probability denoted by bl,t. Subsequently, according to different intervals and their probabilities obtained by the PDF, RWM [26] are applied to generate scenarios for each hour. In this regard, at first, the probabilities of different intervals are normalized in which their summation becomes equal to unity. In addition, each interval is associated with an accumulated normalized probability, as shown in Fig. 3. Therefore, each scenario comprises a vector of binary parameters identifying the load demand, market price, powers output of PV and WT in each period: n o WT PV PV Price Price Scenario s ¼ W L1;t;s ;...;W L7;t;s ;W WT 1;t;s ;...;W 7;t;s ;W 1;t;s ;...;W 7;t;s ;W 1;t;s ;...;W 7;t;s

ð3Þ

ð6Þ

W Ll;t;s ;

W WT uw;t;s ;

W Price pri;t;s

respectively. In other word, for each time interval and each random variable, a random number is produced between 0 and 1 which follows the LMCS strategy. The first interval with an accumulated norselected so that its associated binary parameter becomes equal to 1, the parameters of the non-selected intervals equal to 0, reversely. This trend is continued until the desired number of scenarios is generated. Finally, the normalized probability of each scenario is calculated according to the following equation: t¼1

ld ¼ 1; . . . ; ND; t ¼ 1; . . . ; NT ; s

ld¼1

l¼1

W Ll;t;s bl;t

ps ¼ PNs QN QND P7  T s¼1

¼ 1; . . . ; Ns

W PV pv ;t;s ,

where and are binary parameters indicating whether the lth load interval, the uwth wind power interval, the pvth photovoltaic power output, and the prith market price PV interval which are selected in scenario s (W Ll;t;s ; W WT uw;t;s ; W pv ;t;s , and Price L L WT PV Price W pri;t;s W l;s ¼ 1) or not (W l;t;s ; W uw;t;s ; W pv ;t;s ; and W pri;t;s W Ll;s ¼ 0),

QT QND P7 

þ DP D;ld;t;s ;

t¼1;...;T

malized probability less than or equal to this random number is

WT ¼ 1; . . . ; NWT ; t

¼ 1; . . . ; NT ; s ¼ 1; . . . ; Ns

PD;ld;t;s ¼

1

Interval 2

Interval 3

β 4,t , β 5,t β 6,t , β 7,t

PPV;t;s ¼

Interval 5

0 β 2,t , β 3,t

Interval 4

Interval 3

Interval 1

Interval 2

Interval 1

Probability Density

β 1,t

ð4Þ

t¼1

ld¼1

l¼1

P

W Ll;t;s bl;t

P  P   7 7 PV Price W WT uw;t;s buw;t pv ¼1 W pv ;t;s bpv ;t pri¼1 W pri;t;s bpri;t  P  P   s 7 7 7 WT PV Price uw¼1 W uw;t;s buw;t pv ¼1 W pv ;t;s bpv ;t pri¼1 W pri;t;s bpri;t

7 uw¼1



P

¼ 1;...;N s ð7Þ

Priceutility;t;s ¼

forecast Priceutility;t

þ DPriceutility;t;s ;

¼ 1; . . . ; Ns

t ¼ 1; . . . ; NT ; s ð5Þ

where PWT,t,s, PPV,t,s, PD,ld,t,s, Priceutility,t,s are the wind power output of unit WT (kW), photovoltaic power output of unit PV (kW), ldth load level (kW), and market price (€ct) at time t in scenario s, respecforecast

forecast forecast tively. P forecast PD;ld;t ; Priceutility;t are the forecasted wind WT;t ; P PV;t

power of unit WT (kW), forecasted photovoltaic power of unit PV (kW), the ldth forecasted load level (kW), and forecasted market price (€ct) at time t, respectively. DPWT,t,s (kW), DPPV,t,s (kW), DPD,ld,t,s (kW), and DPriceutility,t,s (€ct) are wind power forecasted

where ND is the total number of load levels in each hour. bl,t, buw,t, bpv,t, bpri,t are the probability of the lth load interval, uwth wind power interval, pvth photovoltaic power interval and prith market price interval, respectively. 2.2. Scenario reduction A higher number of generated scenarios results provide a better modeling of uncertainties but with the cost of higher cumbersome computation. Therefore, scenario-reduction techniques are essential to achieve computation tractability while keeping the stochas-

458

T. Niknam et al. / Applied Energy 99 (2012) 455–470

tic information embedded in the original scenario set as much as possible. The simultaneous backward method is implemented to trim down the number of deteriorated scenarios, as least as possible, this procedure depends on the accuracy of the approximation. Consider ns(s = 1, . . . , Ns) as different Ns scenarios, each with a probability of ps and DT s;s0 as distance of scenario pair (s, s0 ). For using the simultaneous backward method following steps should be applied [24]: Step1: Consider S as an initial set of scenarios, DS is the scenarios which should be deleted. The initial DS is null. Compute the distances of all scenario pairs: DT s;s0 ¼ DTðns ; ns0 Þ; s; s0 ¼ 1; . . . ; N s as ffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pd s s0 2 DT s;s0 ¼ i¼1 ðv i  v i Þ . Step2: For each scenario k, DT k;r ¼ min DT k;s0 ; s0 ; k 2 S and s0 – k, r is the scenario index which has the minimum distance with scenario k. Step3: Calculate PDk,r = pk  DTk,r, k e S. Select d in which PDd = min PDk, k e S. Step4: S = S  {d}, DS = DS + {d}, pr = pr + pd. Step5: Repeat steps 2–4 till the number is mitigated meet our favorite request.

Priceutility;s ¼ ½Priceutility;1;s ; Priceutility;2;s ; . . . ; Priceutility;T;s 1T Putility;s ¼ ½P utility;1;s ; Putility;2;s ; . . . ; Putility;T;s 1T PD;s ¼ ½PD;1;s ; PD;2;s ; . . . ; PD;T;s ð1ðNDTÞÞ ;

PD;t;s

¼ ½P D;1;t;s ; PD;2;t;s ; . . . ; PD;ND;t;s ð1NDÞ PWT;s ¼ ½PWT;1;s ; PWT;2;s ; . . . ; PWT;T;s ð1ðNWT TÞÞ ;

PWT;t;s

¼ ½P1;t;s ; P2;t;s ; . . . ; PNWT ;t;s ð1NWT Þ PBATT;s ¼ ½PBATT;1;s ; PBATT;2;s ; . . . ; PBATT;T;s ð1ðNBATT TÞÞ ;

PFC;s ¼ ½PFC;1;s ; PFC;2;s ; . . . ; PFC;T;s ð1ðNFC TÞÞ ;

PFC;t;s

¼ ½P1;t;s ; P2;t;s ; . . . ; PNFC ;t;s ð1NFC Þ PPV;s ¼ ½PPV;1;s ; PPV;2;s ; . . . ; PPV;T;s ð1ðNPV TÞÞ ;

PPV;t;s

¼ ½P1;t;s ; P 2;t;s ; . . . ; PNPV ;t;s ð1NPV Þ PMT;s ¼ ½PMT;1;s ; PMT;2;s ; . . . ; PMT;T;s ð1ðNMT TÞÞ ;

2.3. Stopping rule

PMT;t;s

¼ ½P1;t;s ; P 2;t;s ; . . . ; PNMT ;t;s ð1NMT Þ

In this paper, a criterion is utilized for choosing a number of samples in order to decide during the course of simulation, whether or not the estimation seems to be accurate enough. This method divides a simulation results into a number of batches, which each batch includes a specific number of scenarios. With this manner it is not essential to check a stopping rule after each scenario because the stopping rule is checked after each batch and if the results are not acceptable, another batch should be run. This criterion is called coefficient of variation cvf and is defined as follows:

cv f ¼

PBATT;t;s

¼ ½P1;t;s ; P2;t;s ; . . . ; PNBATT ;t;s ð1NBATT Þ

rf

pffiffiffiffiffiffi lf N s

ð8Þ

where rf and lf are the standard deviation and mean values of the output random variable f, respectively. If the value of cvf is less than a pre-specific tolerance, then the outcomes can be considered as acceptable results and the simulation can be terminated [28,29].

Us ¼ ½UG;s ; US;s 1ððNG þNS ÞTÞ UG;s ¼ ½UWT;s ;UFC;s ;UPV;s ;UMT;s ;UBATT;s ;Uutility;s 1ððNWT þNFC þNPV þNMT þNBATT ÞTÞ

UWT;s ¼ ½UWT;1;s ; UWT;2;s ; . . . ; UWT;T;s ð1ðNWT TÞÞ ;

UWT;t;s

¼ ½U 1;t;s ; U 2;t;s ; . . . ; U NWT ;t;s ð1NWT Þ UPV;s ¼ ½UPV;1;s ; UPV;2;s ; . . . ; UPV;T;s ð1ðNPV TÞÞ ;

UPV;t;s

¼ ½U 1;t;s ; U 2;t;s ; . . . ; U NPV ;t;s ð1NPV Þ UFC;s ¼ ½UFC;1;s ; UFC;2;s ; . . . ; UFC;T;s ð1ðNFC TÞÞ ;

UFC;t;s

¼ ½U 1;t;s ; U 2;t;s ; . . . ; U NFC ;t;s ð1NFC Þ UMT;s ¼ ½UMT;1;s ; UMT;2;s ; . . . ; UMT;T;s ð1ðNMT TÞÞ ;

UMT;t;s

¼ ½U 1;t;s ; U 2;t;s ; . . . ; U NMT ;t;s ð1NMT Þ 3. Stochastic multi-objective optimal micro-grid operation (SMOMGO) SMOMGO considering hybrid RESs is a nonlinear optimization problem with random continuous and discrete parameters and variables. Penetration of load and generation uncertainties into the problem can affect the power outputs of all other deterministic DGs and storage devices in the MG systems. In the following, the objective functions and constraints are presented. Superscriptsˆ denote the expectation of output random variables. 3.1. Decision variables

X ¼ ½X1 ; X2 ; . . . ; XNs 1Ns nT Xs ¼ ½Ps ; Us ; fs 1nT ;

n ¼ 2NG þ 2NS þ ND þ 2

Ps ¼ ½PFC;s ; PMT;s ; PBATT;s ; Putility;s 1ððNFC þNMT þNBATT þ1ÞTÞ fs ¼ ½PWT;s ; PPV;s ; PD;s ; Priceutility;s 1ððNWT þNPV þNDþ1ÞTÞ

UBATT;s ¼ ½UBATT;1;s ; UBATT;2;s ; . . . ; UBATT;T;s ð1ðNBATT TÞÞ UBATT;t;s ¼ ½U 1;t;s ; U 2;t;s ; . . . ; U NBATT ;t;s ð1NBATT Þ Uutility;s ¼ ½U utility;1;s ; U utility;2;s ; . . . ; U utility;T;s ð1TÞÞ where X is the set of scenarios. NG, NS are the number of generating units and storage devices, respectively. Ps, Us, and fs are the vector of control variables, state variables and output random variables for scenario s, respectively. PFC,s, PMT,s, PBATT,s, Putility,s are the FC, MT, Battery (BATT), and utility production vectors in scenario s, respectively. PWT,s, PPV,s, PD,s, Priceutility,s are the WT production, PV production, load demand, and market price vectors in scenario s, correspondingly. PFC,t,s, PMT,t,s, PBATT,t,s, Putility,t,s are the FC, MT, BATT, and utility production vectors at time t in scenario s, correspondingly. PWT,t,s, PPV,t,s, PD,t,s, Priceutility,t,s are the WT production, PV production, load demand, and market price vectors at time t in scenario s, respectively. UG,s, US,s are the status vectors for generation and storage units in scenario s, respectively. UWT,s, UFC,s, UPV,s, UMT,s, UBATT,s, Uutility,s are the WT, FC, PV, MT, BATT, and utility status vectors in scenario s, correspondingly. UWT,t,s, UFC,t,s, UPV,t,s, UMT,t,s, and

T. Niknam et al. / Applied Energy 99 (2012) 455–470

Fig. 4. Flowchart of the proposed multi-objective ITLBO algorithm.

459

460

T. Niknam et al. / Applied Energy 99 (2012) 455–470

UBATT,t,s are the WT, FC, PV, MT, and BATT status vectors at time t in scenario s, correspondingly. UWT,t,s, UFC,t,s, UPV,t,s, UMT,t,s, UBATT,t,s, and Uutility,t,s are the states of the WT, FC, PV, MT, BATT, and utility at time t in scenario s, respectively.

NG NS ND X X X Pi;t;s U i;t;s þ Pj;t;s U j;t;s þ Putility;t;s U utility;t;s ¼ PD;ld;t;s ; i¼1

ð14Þ

Active power limits of units [5]

Pmin i;t

Minimize [5]

6 Pi;t;s 6 Pmax i;t ;

i ¼ 1; . . . ; NG; t ¼ 1; . . . ; T; s

¼ 1; . . . ; Ns

Ns T X X ps f1;s ðXs Þ ¼ ps Costs;t

ð15Þ

Ns X s¼1

t¼1

s¼1

( T NG X X Ns ¼ ums¼1 ps ½U i;t;s Pi;t;s Bi;t þ Start i maxð0; U i;t;s  U i;t1;s Þ t¼1

i¼1

j¼1

þStart j maxð0; U j;t;s  U j;t1;s Þ þ Shut j maxð0; U j;t1;s  U j;t;s Þ  þU utility;t;s Putility;t;s Priceutility;t;s ð9Þ where Bi,t and Bj,t are the bids of the DGs and storage units at time t, (Starti, Startj) or (Shuti, Shutj) represent the start-up or shut-down costs for the ith DG and the jth storage unit, respectively. Minimize [5] Ns X

Ns X

T X

s¼1

s¼1

t¼1

ps f2;s ðXs Þ ¼

ps

t¼1

Emissions;t

i¼1

þU utility;t;s Putility;t;s Eutility;t



ð10Þ

ð11Þ

where CO2,DG,i,t, SO2,DG,i,t and NOxDG,i,t are the amounts of CO2, SO2 and NOx emissions related to the ith DG unit at time t, respectively[5].

ð12Þ

where CO2,STOR,j,t, SO2,STOR,j,t and NOx,STOR,j,t are the amounts of CO2, SO2 and NOx emissions related to the jth storage unit at time t, respectively [5].

Eutility;t ¼ CO2;utility;t þ SO2;utility;t þ NOx;utility;t

ð16Þ t ¼ 1; . . . ; T; s ¼ 1; . . . ; Ns

ð17Þ

min min where P min i;t ; P j;t ; and P utility;t are the minimum active powers of the ith DG, the jth storage device and the utility at the time t, respecmax max tively. Similarly, P max i;t ; P j;t ; and P utility;t are the capacity powers of corresponding units at time t. Spinning Reserve Requirements (SRRs) [5]

NG NS X X max U i;t;s Pmax þ U j;t;s Pmax i;t j;t þ U utility;t;s P utility;t i¼1

j¼1

t ¼ 1; . . . ; T; s ¼ 1; . . . ; Ns

ð18Þ

where SRt,s is the SRRs at time t in scenario s. Charge and discharge rate limit related to storage device [5]

where Ei,t, Ej,t and Euility,t are described as the amount of pollutants emission in kg MWh1 for each generation unit, storage unit and utility at time t, respectively. These emission variables can be expressed as follows [5]:

Ej;t ¼ CO2;STOR;j;t þ SO2:STOR;j;t þ NOx;STOR;j;t

¼ 1; . . . ; Ns

ld¼1

j¼1

Ei;t ¼ CO2;DG;i;t þ SO2:DG;i;t þ NOx;DG;i;t

j ¼ 1; . . . ; NS; t ¼ 1; . . . ; T; s

ND X P PD;ld;t;s þ SRt;s ;

( Ns T NG NS X X X X ¼ ps ½U i;t;s Pi;t;s Ei;t  þ ½U j;t;s Pj;t;s Ej;t  s¼1

max Pmin j;t 6 P j;t;s 6 P j;t ;

max Pmin utility;t 6 P utility;t;s 6 P utility;t ;

NS X ½U j;t;s Pj;t;s Bj;t þShut i maxð0; U i;t1;s  U i;t;s Þ þ

^f 2 ðXÞ ¼

t

ld¼1

¼ 1; . . . ; T; s ¼ 1; . . . ; Ns

3.2. Objective functions

^f 1 ðXÞ ¼

j¼1

ð13Þ

where CO2,utility,t, SO2,utility,t and NOx,utility,t are the amounts of CO2, SO2 and NOx emissions related to utility at time t, respectively. The multi-objective SMOMGO problem consists of two competing objectives, i.e., expected total operation costs (9) and expected total operation emissions (10). The total operating cost of the MG in €ct (Euro cent) includes the fuel costs of DGs, start-up/shutdown costs and the costs of power exchange between the MG and the utility [5]. The emission of pollutant gases (10) is modeled in which three of the most important pollutants are involved: CO2 (carbon dioxide), SO2 (sulfur dioxide) and NOx (nitrogen oxides) [5]. 3.3. Constraints [5] The constraints of the proposed problem are listed as follows: Power balance [5]

wj;t;s ¼ wj;t1;s þ gcharge;j Pcharge;j;t;s Dt 

1

gdischarge;j

Pdischarge;j;t;s Dt;

j

¼ 1; . . . ; NS; t ¼ 1; . . . ; T; s ¼ 1; . . . ; Ns ; Dt ¼ 1h

ð19Þ

6 wj;t;s 6 wmax ; j ¼ 1; . . . ; NS; t ¼ 1; . . . ; T; s wmin j j ¼ 1; . . . ; Ns Pcharge;j;t;s 6 Pmax charge;j ;

ð20Þ j ¼ 1; . . . ; NS; t ¼ 1; . . . ; T; s

¼ 1; . . . ; Ns Pdischarge;j;t;s 6 Pmax discharge;j ;

ð21Þ j ¼ 1; . . . ; NS; t ¼ 1; . . . ; T; s

¼ 1; . . . ; Ns

ð22Þ

where, wj,t,s and wj,t1,s are the amounts of energy storage inside the battery storage j at time t and t  1 in scenario s, respectively. Pcharge,j,t,s(Pdischarge,j,t,s) is the jth battery permitted rate of charge (discharge) during a definite period of time (Dt = 1 h) at time t in scenario s. gcharge,j(gdischarge,j) is the efficiency of the jth battery during the charge(discharge) process. wmin and wmax are the lower and j j upper limits on amount of energy storage inside the jth battery. max Pmax charge;j and P discharge;j are the maximum rate of the jth battery charge(discharge) during each time interval Dt = 1h, respectively. 4. Solution methodology In this section, the ITLBO algorithm is presented followed by the multi-objective solution methodology including a fuzzy model for objective functions and Pareto optimization method. Thirdly, the application of the proposed algorithm is presented. Finally, some applicable tool usages of the SMOMGO problem are provided. The flowchart of the whole process is given in Fig. 4.

461

T. Niknam et al. / Applied Energy 99 (2012) 455–470

4.1. Improved teaching–learning-based optimization (ITLBO)

xyiter m;disc;jj;new1 ¼

TLBO is a new evolutionary algorithm which uses a population of learners to proceed to the global solution [22]. In this algorithm, the population is considered as a class of learners. In TLBO, different design variables are resembled different subjects offered to learners and the learners’ remarks are resembled the ‘fitness function’. The teacher is considered as the best solution obtained so far. Learners can improve their knowledge through teacher teaching and interaction between themselves. Therefore, the TLBO’s process is divided into two parts include teacher part and learners part. 4.1.1. Teacher part A good teacher can bring up the mean of the class to some extend based on their talents and abilities [22]. To clarify this process, the following parameters are defined: MEiter is a vector including the mean of each decision variable for all learners and Titer is the teacher at the iterth iteration of the evolution procedure. Based on the TLBO concept, Titer tries to move mean MEiter up to its own grade. The structure of each learner is formed by putting the discrete and continuous control variables beside each other as follows:

h i iter iter Xiter m ¼ Xm;disc ; Xm;cont h i iter iter iter ¼ xiter m;disc;1 ; . . . ; xm;dis;Ndisc ; xm;cont;1 ; . . . ; xm;cont;Ncont

ð23Þ

xxiter m;disc;jj;new1  xxmin xxmax  xxmin

;

m ¼ 1; . . . ; Nlearner ; jj

¼ 1; . . . ; Ndisc

ð32Þ

iter where xxiter m;disc;jj;new1 and xym;disc;jj;new1 are the pseudo-probability and the actual probability for the jjth discrete variable of new solution of the mth learner, respectively. xxmin and xxmax are the predetermined as the minimum and maximum value of xxiter m;disc;jj;new1 , correspondingly. It is notable that the initial value of xx0m;disc;jj parameters are equal to zero and xxmin and xxmax are equal to 50 and 50, respectively. Thus, the initial values for xy0m;disc;jj;new1 are 0.5 for m = 1, . . . , Nlearner and jj = 1, . . . , Ndisc. After that, a random number (rand) is generated and compared iter with the probability xyiter m;disc;jj;new1 . The value of xm;disc;jj;new1 will be

determined as follows [30]: ( 1 if rand 6 xyiter m;disc;jj;new1 ; iter m ¼ 1;...;N learner ; jj ¼ 1;...;N disc xm;disc;jj;new1 ¼ 0 else; ð33Þ

The final control vector is generated by putting the discrete and continues control variables beside each other as follows:

h i iter iter Xiter m;new1 ¼ Xm;disc;new1 ;Xm;cont;new1 h i iter iter iter ¼ xiter m;disc;1;new1 ;...;xm;dis;N disc ;new1 ;xm;cont;1;new1 ;...;xm;cont;N cont ;new1 ð34Þ

where Ndisc and Ncont are the number of discrete and continuous variables, respectively. The structure of the teacher and the mean value of the class are defined as follows:

ME

iter

¼

h

iter MEiter disc ; MEcont

i

h i iter iter iter ¼ meiter disc;1 ; . . . ; medisc;Ndisc ; mecont;1 ; . . . ; mecont;Ncont meiter disc;jj ¼

meiter cont;ii ¼

xiter 1;disc;jj

þ

After generating all vectors, they are compared with the old one according to the min–max method as follows [31,32]:

xiter 2;disc;jj

þ  þ

xiter Nlearner ;disc;jj

Nlearner iter iter xiter 1;cont;ii þ x2;cont;ii þ    þ xNlearner ;cont;ii

Nlearner

ð24Þ

ð25Þ

ð26Þ

Since the proposed problem has both discrete and continuous control variables so the generations of population involves two steps. The continuous solutions (learners) can be updated according to the following equation [22]:

  iter DMiter ¼ randðÞ Titer  T iter F M

ð27Þ

h i iter DMiter ¼ DMiter disc ; DMcont h i iter iter iter iter ¼ dmdisc;1 ; . . . ; dmdisc;Ndisc ; dmcont;1 ; . . . ; dmcont;Ncont

ð28Þ

¼ roundð1 þ randÞ T iter F

ð29Þ

T iter F

where is the teaching factor in the iterth iteration. rand is a random function generator in the range [0, 1]. Now, the new solution can be produced according to the following equation [22]: iter iter Xiter m;cont;new1 ¼ Xm;cont þ DMcont ;

m ¼ 1; . . . ; Nlearner

ð30Þ

The discrete variables (learners) can be updated as follows [30]: iter xxiter m;disc;jj;new1 ¼ xxm;disc;jj þ

¼ 1; . . . ; Ndisc

iter dmdisc;jj ;

m ¼ 1; . . . ; Nlearner ; jj ð31Þ

l^ fq ðXÞ ¼

8 > 1 > > > < ^max

^f q ðXÞ 6 ^f min ; q

f q ^f q ðXÞ ^f max ^f min

> q > > > :0

q

^f min < ^f q ðXÞ 6 ^f max ; q i

q ¼ 1; 2

ð35Þ

^f q ðXÞ > ^f max ; q

^f ¼ Maxððl ^ ^ ref ;1  lf1 Þ; ðlref ;2  lf2 ÞÞ

ð36Þ

where ^f min and ^f max are the minimum and maximum acceptable exq q ^ fq is the expected level of the qth objective function, respectively. l pected membership value for the qth objective function. lref,qq = 1, 2 is the reference membership value for the qth objective function. If the new learner has lower value of Eq. (36) then it substitutes the old one. It should be noted that the final results of this part are the inputs for the next part, i.e. learner part. 4.1.2. Learner part In this part, interaction between learners can improve their knowledge (situations). Each learner interacts with the other learners randomly by the use of group discussions, presentations, formal communications, etc. [22]. Thus, each learner can gain knowledge of something new if the other ones have more knowledge than him/her. This procedure is described in the following form [22]. For the mth learner in the class, an individual (n) is selected randomly in a way that m – n. Now a new solution is generated as follows [22]:

      iter iter iter Xiter if ^f Xiter < ^f Xiter m;new2 ¼ Xm þ randðÞ Xm  Xn m n   iter iter iter otherwise Xiter m;new2 ¼ Xm þ randðÞ Xn  Xm

ð37Þ

iter iter After carrying out this part, Xiter m is replaced by Xm;new2 if Xm;new2 has a better function value according to Eq. (36). Note that like the ‘teacher part’, the final results of the learner part are the inputs for the next part i.e. improved part.

462

T. Niknam et al. / Applied Energy 99 (2012) 455–470

4.1.3. Improved part In comparison with other evolutionary algorithms, TLBO has many major advantages that have boosted its usage in several approaches such as economic dispatch applications [33]. These vantages are its easy implementation, simple concept, minimal storage requirements and no need for tuning the algorithm parameters. It means this algorithm reaches the optimal solution without any user defined parameter and this feature is a major advantage of this algorithm. Despite these features, TLBO may experience unsuitable convergence, so this paper devises the following modification to cope with this problem and avoid trapping in local optima. Since the proposed problem has continuous and discrete control variables so like the generation of each vector, this section has different rules for generating new vectors for continuous and discrete variables. Concerning the continuous variables, two different modified rules are applied which each learner according to a probability model chooses one of these methods. 4.1.4. Modification method 1 This modification can be considered as globally correlated and therefore effectively enhances the TLBO’s global optimization capabilities as follows [34]: iter Xiter mm;cont;method1 ¼ Xmm;cont

  iter þ rand Titer cont  Xmm;cont ;

mm

¼ 1; . . . ; N1

ð38Þ

Xiter mm;cont;method1

where is the new solution which is generated from modification method1. Titer cont is the continues variables part of teacher Titer. N1 is the number of learners which select the modification method 1. 4.1.5. Modification method 2 According to [34], the trigonometric modification operation of Eq. (39) biases the new trial solution heavily in the direction of the best one of three individuals chosen for modification process. It is worthwhile to note that it is a local search operator. In this way, the proposed modification operator fulfills a wider search in the solution space and has significant effect in many generations along the ITLBO algorithm. This method is formulated as follows [34]:

Xiter mm;cont;method2 ¼

! iter iter Xiter r1 þ Xr2 þ Xr3 þ ðR1 3     iter iter  R2 Þ Xiter þ ðR3  R2 Þ Xiter r1  Xr2 r2  Xr3   iter þ ðR1  R3 Þ Xiter ; mm r3  Xr1

¼ 1; . . . ; N2 j^f ðXiter r1 Þj

j^f ðXiter Þj R2 ¼ Rr20 ,

ð39Þ j^f ðXiter r3 Þj

     R ¼ ^f Xiter þ r1 0

where, R1 ¼ R0 , R3 ¼ R0 ,       ^ iter  ^ iter  iter f Xr2  þ f Xr3 . Xmm;cont;method2 is the new solution which is generated from modification method2. N2 is the number of learners which select the modification method 2. In each iteration (iter) and for each solution mm, three vectors, r1, r2, r3, are selected from the existing population such that r1 – r2 – r3 – mm. The probabilistic solution strategy is defined as follows: at first the probability of both modification methods are considered probmethod = 0.5 (method = 1, 2) and a parameter called accumulator is allocated to each of them in which acummethod = 0 (method = 1, 2). In each iteration, a weight factor is allocated to each learner after sorting the population according to Eq. (40). It is clear that the best learner gets the larger weight factor. After that the related accumulator of each method is brought up to date as Eq. (41):

wwm ¼

logðNlearner  m þ 1Þ ; logð1Þ þ    þ logðNlernear Þ

acummethod ¼ acummethod þ

wwmm ; Nmethod

m ¼ 1; . . . ; Nlearner

ð40Þ

mm ¼ 1; . . . ; Nmethod

ð41Þ

where Nmethod is a number of learners which select methodth modification method and wwmm (mm = 1, . . . , Nmethod) are the weight factors corresponding to them. The excitation probability is calculated as:

probmethod ¼ ð1  aÞ  probmethod þ a 

acummethod ðmethod Iter max

¼ 1; 2Þ

ð42Þ

where a is a learning rate to control the learning speed in the ITLBO algorithm and it is considered equal to a = 0.142 in this paper. Finally the RWM is applied to choose the methodth modification method for each learner based on normalized probability values as follows:

probmethod ¼

probmethod ðmethod ¼ 1; 2Þ prob1 þ prob2

ð43Þ

After the above procedure the new solution is generated for each learner m as Xiter m;cont;new3 . The proposed modification for all discrete variables has simple mechanism in which for each learner the discrete value is updated as follows [30]:

if rand < 0:08( xiter m;disc;jj;new3 ¼

1 if

xiter m;disc;jj ¼ 0

0

xiter m;disc;jj ¼ 1

if

ð44Þ

After generating discrete and continuous modification vector, they are put beside each other for generating a control variable vector. Eq. (36) is computed for each modification vector. ^f is considered as a criterion for substituting the modification vector with the relevant vector which had been generated by the learner section. 4.2. Multi-objective optimization problem (MOP) 4.2.1. Pareto-based approach Multi-objective optimization methods are specific tools to explore the search space of optimization problem for optimal solution/solutions in a reasonable time. The proposed Pareto method obtains several non-inferior solutions of the optimization problem in a single run and preserves the diversity of Pareto optimal solutions besides determines the best one. The Pareto method is based on the non-dominated solutions concept. In a minimization problem, a solution X1 dominates X2 if and only if the two following conditions are satisfied simultaneously [23,33,35]:

8q 2 f1; 2g : ^f q ðX1 Þ 6 ^f q ðX2 Þ 9q 2 f1; 2g : ^f q ðX1 Þ < ^f q ðX 2 Þ

ð45Þ

X1 is a non-dominated solution within the set {X1, X2}, if X1 dominates the solution X2. The solutions that are non-dominated within the entire search space are denoted as Pareto-optimal solutions. These non-dominated solutions are stored in the repository during the algorithm search process. 4.2.2. Fuzzy-based clustering to reduce the repository size The POF for most of the problems is enormously large and may even include an infinite number of individuals. The large number of Pareto-optimal solutions leads to more computation time [23]. Moreover, there are always memory constraints. Therefore, decreasing the size of the repository without destroying the characteristic of the POF is necessary. Fuzzy clustering method is applied for fixing the size of repository, it is worthwhile to note

T. Niknam et al. / Applied Energy 99 (2012) 455–470

that the applied method dose not destroy the characteristics of the POF [23]. In this regard, when the number of solutions in the repository reaches its maximum limit, pruning all non-dominated solutions (members of the repository and non-dominated solutions wishing to enter the repository) is done by the proposed fuzzy clustering method. In order to generate a compromise value for each member of the repository, the MGCC is asked to specify the desired weight factor for each of the objective functions which represent the importance of each objective named as wq, q = 1, 2. Then the following normalized weighted membership approach can be used to extract better solutions in the premiere region with iteratively search process which it is expected to find the closest solutions with respect to the MGCC’s requirements [23]:

P2 ^ q¼1 wq lfq ðn1 Þ Nl ðn1 Þ ¼ PNrep P2 ^ q¼1 wq lfq ðn1 Þ n1 ¼1

ð46Þ

In the proposed clustering approach, at first each non-dominated solution builds a separate cluster. After that the solutions which have small distance with this cluster adjoined to this cluster until the number of clusters becomes equal to the repository size. At this time, an individual with higher normalized weighted membership value from each cluster is selected to be stored in the repository. 4.2.3. Niching technique for diversity preservation Niching technique [23], is utilized in this paper for avoiding the learners accumulating towards densely populated area and consequently obtaining a uniformly distributed POF. This technique implements the fitness sharing which is based on the fundamental idea that the resources in a particular niche have to be shared between the present individuals [23]. Therefore, the higher number of the neighboring individuals leads to more inferiority of the fitness. In the proposed approach, updating rule of the teacher part is changed by incorporating a niching-based procedure that guides the learners towards sparse regions. In the proposed method, the position of the teacher is selected from the repository by usage of the niching method. The Pseudocode of this process is shown at the end this section. In this regard, sharing distance dn1 n2 between repository member n1 and n2 should be computed. The Shðdn1 n Þ which is the corner stone of the niching method computed as follow [23]:

Shðdn1 n2 Þ ¼

8 <

d 2 n n 1 r1 2 share : 0

if dn1 n2 6 rshare

ð47Þ

otherwise

where rshare is a niche radius which is set to 2 in this paper. After computing the Sh(dij), the niche count M n1 can be computed and after that fitness of each individual which is inversely proportional to its niche count will be calculated. The global best solution in the repository is selected by the RWM [36]. The niching approach can be expressed as follows [23]: For n1 ¼ 1 to N rep For n2 ¼ 1 to N rep Calculate sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ^f q ðXn Þ^f q ðXn Þ 2 P2 1 2 and Shðdn1 n2 Þ according to (47). dn1 n2 ¼ q¼1 ^max ^min fq

M n1 ¼ M n1 þ Shðdn1 n2 Þ

f q

EndFor (it refers to set n2)

fit n1 ¼ M1n EndFor (it refers to 1

set n1)

Normalize the fitness function and use RWM to select the Titer. The initial values for Mn1 ; n1 ¼ 1; . . . ; N rep are zero. Nrep is the number of repository members. fit n1 ; Mn1 are the fitness and niche count of the n1th repository member, respectively. 4.3. Application of the proposed algorithm Step 1: Input all required data for the system

463

Step 2: Scenario generation: Generate 1000 scenarios and reduce them with the backward scenario reduction to reach the desired batch including Ns scenarios. This batch should be satisfied Eq. (8). For each scenario, the hourly load demand, market price, powers output of PV and WT units are deterministic. Step 3: Represent and initialize the learners for the population: A learner is a vector representing the generation levels and statuses of each unit (i.e. MT, Phosphoric Acid Fuel Cell (PAFC), Battery, and utility) at each hour of the day. At the beginning, the learners of the population are generated randomly in the range of [0, 1] and located between the maximum and minimum operating limits of the unit. For each learner in the population, total operation costs f1,s s = 1, . . . , Ns and total operation emissions f2,s s = 1, . . . , Ns are calculated according to Eqs. (9) and (10), respectively. It should be noted that all the constraints should be satisfied for each scenario contributed to each learner. A penalty term is added to f1,s and f2,s for each violation of the constraints. Aggregated total operation costs and emissions for each learner are calculated using Eqs. (9) and (10), respectively. After that, compute Eq. (36). Step 4: Constitute the repository: Save the non-dominated solutions of the current population in the repository. Step 5: Reduce the repository size: Use the fuzzy clustering mechanism based on Section 4.2.2 to achieve the desired repository size. Step 6: Update: Titer: The niching mechanism is applied on the set of non-dominated solutions in the repository and the best compromise solution is selected as Titer. Step 7: Incorporate teacher part as described in Section 4.1.1. Step 8: Incorporate learner part: The learners would try to improve themselves by the use of knowledge interaction as described in Section 4.1.2. Step 9: Use SAPMS by incorporating modification part as described in Section 4.1.3. Step 10: Update the repository: Check for non-domination. Update the repository. Step 11: Check the convergence criteria: Go to step 5 for the next iteration. This loop can be terminated after a predetermined number of iterations and the teacher is selected as the expected best compromise solution. 4.4. Tool usage If the uncertainty of load demands, market prices, and available powers output of PV and WT units are neglected by the MGCC and the deterministic MOMGO is solved instead of the stochastic MOMGO problem, the MGCC will almost certainly fail to meet the challenges of supplying the load demands and satisfying the problem constraints. By the time that the dispatch results of MOMGO problem must be taken place at the beginning of the next day. Therefore, the sources of the system uncertainties should be considered in order to avoid inconsistency between the real conditions and ante-scheduling patterns of MOMGO problem. Additionally, in each time period, the forecast should be updated and a new SMOMGO should be run while the powers output of all DG sources, Nickel-Metal-Hydride-Battery (NiMH-Battery) storage and purchase power from utility in the previous hour are taken into account. Consequently, to handle the system uncertainties as well as the constraints of the proposed framework, implementing the SMOMGO to cope with the stochastic nature of MOMGO problem is a vital importance. In the conventional vertically integrated utilities, MGCCs are in charge of load supply and generation scheduling, simultaneously. Therefore, to cope with the load and market price uncertainties be-

464

T. Niknam et al. / Applied Energy 99 (2012) 455–470

Fig. 5. Typical low-voltage MG model.

Table 1 Results obtained by different methods for DOMGO1 (S1). Solution technique

Total operation costs (€ct) Best value

Mean value

Worst value

GA [37] PSO [37] FSAPSO [37] CPSO-T [37] CPSO-L [37] AMPSO-T [37] AMPSO-L [37] TLBO TLBO-method1 TLBO-method2

277.7444 277.3237 276.7867 275.0455 274.7438 274.4317 274.5507 272.7610 271.1361 270.9328

290.4321 288.8761 280.6844 277.4045 276.3327 274.5643 274.9821 273.1998 271.7890 271.2512

ITLBO

269.7600

269.7600

Table 2 Results obtained by different methods for DOMGO1 (S2). CPU time (min)

Solution technique

Total operation costs (€ct) Best value

Mean value

Worst value

304.5889 303.3791 291.7562 286.5409 281.1187 274.7318 275.0905 273.8411 272.0866 271.6319

NA NA NA NA NA NA NA 0.005 0.015 0.013

GA PSO TLBO TLBO-method1 TLBO-method2

298.1642 297.7435 291.6320 288.8204 288.3759

310.8519 309.2959 292.8431 289.4275 288.7199

325.0587 323.7989 293.5299 290.0067 289.0355

0.323 0.296 0.006 0.020 0.017

ITLBO

287.4798

287.4798

287.4798

0.014

269.7600

0.012

NA: not available in the referred literature.

side PV and WT power generation variations, the proposed SMOMGO framework is suitable for the DMs of the system. Also, in the deregulated power systems, both of MGCCs and upstream networks can use this framework. The MGCCs can use the proposed stochastic framework to handle their uncertainties and also, the upstream networks can implement the proposed framework for bidding strategy to manage their uncertainties regarding load, market price, PV and WT forecasted errors. Therefore, utilities can adopt appropriate strategy to participate in power markets while they consider the probable scenarios.

5. Simulation and numerical results In order to verify the effectiveness and superiority of the proposed ITLBO algorithm for the SMOMGO problem, it is applied on a typical low voltage MG. This system contains various DG sources

CPU time (min)

Table 3 Results obtained by different methods for DOMGO1 (S3). Solution technique

Total operation costs (€ct) Best value

Mean value

Worst value

CPU time (min)

GA PSO TLBO TLBO-method 1 TLBO-method 2

334.3679 327.9226 313.8403 306.1463 304.1127

336.3074 331.0284 315.7245 307.4896 305.5970

345.3569 340.2094 316.9180 308.5608 306.3305

1.294 1.040 0.021 0.065 0.061

ITLBO

302.2767

302.2783

302.2812

0.056

such as MT, PAFC, PV, WT, and a NiMH-Battery as portrayed in Fig. 5. The system data is extracted from [5], where a complete data set can be found. The SRRs are set to 5% of the load demand in each hour. It should be noted that a time period of 1 day with hourly time step is considered as the study horizon. All the DG units produce active power at unity power factor. It means that these units require or produce no reactive power. Also, the thermal load is not

465

T. Niknam et al. / Applied Energy 99 (2012) 455–470

Fig. 6. Generation levels of different units for DOMGO1 (S1) in kW.

Fig. 8. Generation levels of different units for DOMGO1 (S3) in kW.

650 600

Deterministic POF Stochastic POF

550 500

Emission (Kg)

considered in the proposed MG system. Furthermore, there is a power exchange link between the utility and the MG during the time step in the study period based on the decisions made by the MGCC. In this paper, it is assumed that the MGCC buys the available powers of the PV and WT units at each hour. In order to better exposition of our proposed framework, three different strategies are considered. In the first Strategy (S1), all the DG units are in service state during the time horizon and also the initial charge of the NiMH-Battery is infinite. In the second Strategy (S2), the NiMH-Battery starts with the 100 kW initial charge and all the units are in service. In the third Strategy (S3), the NiMH-Battery starts with no charge at the beginning of the first hour of the day and all units are allowed to start up or shut down for flexible operation of the MG. For comparison purposes, the single objective versions of SMOMGO have been solved, named the SOMGO1, referred to cost, and the SOMGO2, referred to emission. Similarly, the deterministic version of SMOMGO, referred to as DMOMGO, has also been solved. The simulations are carried out on a Pentium P4, Core 2 Duo 2.4 GHz personal computer with 1 GB RAM memory. The software is developed using MATLAB 7.8. The parameter setting of the optimization algorithm is as follows: (i) the number of the population is equal to 40. (ii) Maximum number of the iterations is 50. It should be pointed out that all evolutionary methods require tuning of different algorithm parameters for their proper searching. A small change in these parameters may result in a large change in the algorithm performance. TLBO overcomes such difficulties, because it does not require any parameter for tuning. It means this algorithm reaches the optimal solution without adjusting any parameter and this feature is a significant advantage of this algorithm. For a better illustration of the performance of the proposed method in the stochastic optimization, six batches of scenarios are considered. They include 6  30, 5  30, 4  30, 3  30, 2  30, and 30 samples of scenarios, respectively. As mentioned before, a scenario is constructed by samples of the market price, load demand, and available output powers of the PV and WT units for all hours in the time period. Thus, each scenario has 82 parts. Obviously, the computational effort of the scenario-based approach is directly proportional to the number of scenarios in each batch so the backward scenario reduction method would reduce the total number of scenarios considering a tradeoff between CPU time and solution accuracy. It should be mentioned that due to alternative nature of the evolutionary algorithm, 30 independent trial runs are carried out to extract the statistical information such as best, worst, and mean value of the obtained solution.

Fig. 7. Generation levels of different units for DOMGO1 (S2) in kW.

450 400

Expected best compromise value

350 300 250 200 150 200

400

600

800

1000

1200

1400

1600

Cost (Ect) Fig. 9. Deterministic and stochastic POF achieved by ITLBO algorithm for S1 (80 non-dominated solutions).

5.1. Deterministic single objective and multi-objective solutions of the DMOMGO This case is a base case which the input random variables are considered in the deterministic quantities. Firstly, the single

466

T. Niknam et al. / Applied Energy 99 (2012) 455–470 750 Deterministic POF Stochastic POF

700

Emission (Kg)

650 600 550 500 450 400 350 200

400

600

800

1000

1200

1400

1600

Cost (Ect) Fig. 10. Deterministic and stochastic POF achieved by ITLBO algorithm for S2 (80 non-dominated solutions).

Fig. 12. Generation levels of different units for the best compromise solution (S1) in kW.

750 Deterministic POF Stochastic POF

700

Emission (Kg)

650 600 550 500 450 400 200

400

600

800

1000

1200

1400

Cost (Ect) Fig. 11. Deterministic and stochastic achieved by ITLBO algorithm POF for S3 (80 non-dominated solutions).

objective versions of DMOMGO have been solved, namely DOMGO1, referred to cost, and DOMGO2, referred to emission. This procedure should be done to extract the extreme points of the tradeoff surface and judge the diversity characteristics of the POF. The quality of the solutions found by ITLBO and TLBO has been assessed through the comparison with the results reported in the literature by available methods [37], which did not model SRRs. Since this paper considers the uncertainty in load demands, the market prices, and powers output of PV and WT units so, for reliability reasons, reserve consideration is crucial. As can be inferred from Tables 1–3, the superiority of ITLBO over the other approaches is substantiated by the achievement of better values for cost objective function. It is obvious that the SAPMS is a useful tool which not only accelerates the conversions of the original TLBO, but also enables it to alleviate the stagnation and escape from local optima. Above tables confirm that the proposed ITLBO provides high quality solutions and faster convergence in a confidence manner in comparison with the other cited methods in the area. To show the effect of the proposed modification methods 1 and 2, the results of the optimization of each separate modification operator in different strategy cases (S1, S2 and S3) have been added to Tables 1–3, correspondingly. It is necessary to note that the basic idea behind this approach is to simultaneously select adaptively multiple effective operators from the candidate modification pool on the

Fig. 13. Generation levels of different units for the best compromise solution (S2) in kW.

Fig. 14. Generation levels of different units for the best compromise solution (S2) in kW.

basis of their previous experiences in the generated promising solutions and applied to perform the modification operation. Besides, the ITLBO with SAPMS can better manage transition from

467

T. Niknam et al. / Applied Energy 99 (2012) 455–470 Table 4 Comparison of the results obtained by different scenarios for SOMGO1 (S1). No. of scenarios

Total operation costs (€ct) Mean

SD

0.95% confidence intervals

Expected

VSS

1  30 2  30 3  30 4  30 5  30 6  30

272.9092 270.2148 270.2008 270.2118 270.2092 270.3120

5.9814 6.5252 6.4397 7.03878 7.5787 8.1009

2.1404 1.6511 1.3305 1.2594 1.2128 1.1835

274.6370 274.5769 274.5701 274.5676 274.5686 274.5686

4.8769 4.8169 4.8101 4.8076 4.8085 4.8085

Relative error (%)

Coefficient of variation (%)

CPU time (min)

0.7794 0.6013 0.4846 0.4587 0.4417 0.4310

0.4002 0.3118 0.2512 0.2378 0.2290 0.2234

0.352 0.684 1.014 1.332 1.695 1.988

Relative error (%)

Coefficient of variation (%)

CPU time (min)

0.7499 0.5670 0.4685 0.4388 0.4188 0.4099

0.3852 0.2936 0.2426 0.2271 0.2168 0.2121

0.390 0.785 1.212 1.583 1.942 2.396

Table 5 Comparison of the results obtained by different scenarios for SOMGO1 (S2). No. of scenarios

1  30 2  30 3  30 4  30 5  30 6  30

Total operation costs (€ct) Mean

SD

0.95% confidence intervals

Expected

VSS

291.5444 289.0636 288.9976 289.0253 289.0339 289.1054

6.1512 6.5743 6.6506 7.1904 7.6747 8.2284

2.2012 1.6635 1.3740 1.2865 1.2282 1.2021

293.5038 293.3664 293.2640 293.1851 293.2628 293.2633

6.0239 5.8866 5.7842 5.7052 5.7829 5.7834

Table 6 Comparison of the results obtained by different scenarios for SOMGO1 (S3). No. of scenarios

1  30 2  30 3  30 4  30 5  30 6  30

Total operation costs (€ct) Mean

SD

0.95% confidence intervals

Expected

VSS

Risk

306.5437 304.4576 303.9253 303.9592 304.0785 304.2799

6.7752 7.2178 7.1983 8.0202 8.9822 9.7806

2.4245 1.8263 1.4872 1.4350 1.4374 1.4288

308.5888 308.5241 308.4972 308.4946 308.4983 308.5006

6.3120 6.2474 6.2204 6.2179 6.2216 6.2239

2.0451 4.0665 4.5719 4.5354 4.4198 4.2207

each generation to the next one in comparison with each separate modification operator. The optimal generation levels of the different units for the case of DOMGO1 are shown in Figs. 6–8 for S1, S2, and S3, correspondingly. In the first two strategies (S1 and S2), all units are in service, hence the MGCC has to buy at least minimum power output of all units (because they are in service and they have to generate electricity power) even though this procedure burdens extra cost. Since the MT units are expensive so, they are restricted to their minimum value during most hours in S1. But the procedure is different in S2 and S3 scenarios, since the discharge action of the battery is restricted so, the MGCC should buy from the MT unit. It is clear that in all strategies economical choice such as PAFC and utility are dispatched as much as possible to satisfy the system load demand. Also, as shown in these figures, battery saves low-price energy of the utility at light-load hours and releases this energy during the peak-load hours which prepares good efficiency for the MGCC. The MGCC purchases power from the utility at the first hours in order to fulfill load demand, meantime the extra power is saved in battery. Meanwhile, the discharge procedure of battery is started form hour 9 (peak load) to fulfill the load demand as well as selling of the power to upstream network at high price. This routine is altered in the 17th–19th hours. Also, the battery is completely discharged in hour 20. In state of DOMGO2 which refers to the emission objective function, above debate can also be concluded. After determining the extreme points, the ITLBO is applied to optimize operation costs and emissions simultaneously in order to achieve POF. Figs. 9–11 show POF achieved in one single run for S1, S2, and S3, respectively. These POF are extracted in 0.018, 0.023, and 0.068 min, correspondingly. These figures demonstrate

Relative error (%)

Coefficient of variation (%)

CPU time (min)

0.7857 0.5919 0.4821 0.4652 0.4659 0.4631

0.4035 0.3061 0.2496 0.2409 0.2412 0.2396

1.330 2.891 4.840 6.422 8.122 9.790

that the POFs achieved by the proposed method are well-distributed and have satisfactory diversity which covers the entire POF effectively. The optimal generation levels of the different units for the best compromise solution are shown in Figs. 12–14 for S1, S2, and S3, correspondingly. 5.2. Stochastic single objective and multi-objective solutions of the SMOMGO From the point of MG operation, the stochastic approach by consideration of the intermittent nature of system components and loads can provide a more accurate solution for determining the allowance and optimum costs and emissions. The stochastic single objective SOMGO1&2 are firstly solved by usage of the scenario generation scheme described in Section 2.1. 1000 24-h scenarios [38] are generated and subsequently trimmed to six batches. The number of scenarios is usually selected in a way that the coefficient of variation becomes small i.e. between 0.1% and 1% [39]. Tables 4– 6 show the coefficient of variation for the aforementioned batches for S1, S2, and S3, respectively. As can be seen from these Tables, the total operation costs are decreased with increasing the number of the scenarios. The indexes such as the value of expected, mean, Standard Deviation (SD), 95% confidence intervals [24], relative error [24], coefficient of variation [28], and the Value of the Stochastic Solution (VSS) [40] for each batches are also shown in these Tables. The values of all indexes except SD are decreased slightly when the number of samples is increased. In order to quantify the relative impact of the uncertainty sources on the decision making of the MGCC, the VSS can be implemented. This metric measures how much the MGCC would be willing to spend cost and emission to know the

468

T. Niknam et al. / Applied Energy 99 (2012) 455–470 12

11 10 9

10

8 7

8

6 5

6

4 3 2

4

1 0 580

585

590

595

600

605

610

615

Fig. 15. PDF of total operation costs for the expected best compromise solution of S1.

future realizations of the under consideration stochastic problem. The VSS is the difference between the optimal solutions of the single objectives DMOMGO and SMOMGO. The additional cost and emission caused by considering of the uncertainties. In the other word, modeling the system uncertainties increases the operation cost and emission values as well as CPU time because of considering different most probable scenarios in the stochastic framework instead of single scenario in the deterministic scheme. Also, it should be noted that lower the value of confidence interval shows the precise degree of the expected value. When the interval is smaller, we expect the real time solution to be closer to expected value [24]. In addition, the relative error is computed as ðð95% confidence intervalÞ=expected valueÞ  100%. The smaller the relative error, the more efficient the proposed method is [24]. It should be noted that the main disadvantage of the batches with enormous scenarios is the great CPU computational time requirement. With the help of fast and efficient computing tools i.e. scenario reduction technique and the ITLBO algorithm, a huge stochastic model can be easily solved. From observing Tables 4–6, it can be found that the batch which includes 30 scenarios is computationally significantly faster the other batches while its indexes are also near the another batches. The values of the above tables show that the major difficulty with scenario reduction from 180 to 30 leads to a little increase in the total operation costs. Although using the 30 scenario increase 0.025%, 0.082%, and 0.029% operational costs of SOMGO1 and the CPU time is decreased 82%, 84%, and 86% for S1, S2, and S3, respectively. The tradeoff between the operation costs and the CPU execution time reveals the superiority of the proposed backward scenario reduction. Thus, the 30 scenarios are considered for the SMOMGO problem. In this paper, the stochastic approach tries to find expected POF solutions. The stochastic solutions may not be global optima solutions to the individual scenarios but they are a robust and also located near global solutions which this provides possible realizations of the uncertainties. To this end, the Risk Index (RI) is calculated and added to Table 6 for S3. The RI is the average value of the difference between the objective function value of the proposed approach and the deterministic approach which is specified as an index to compare the performance of the stochastic approaches. It is based on the difference between the objective function value of deterministic solution of the proposed problem corresponding to each scenario fq,s([Ps, Us, fs]) and the objective function value of the proposed problem of each scenario which is obtained by unit scheduling from applying the proposed approach fq,s([Ps, U, fs]) and is defined as follows [29]:

2

0 675

680

685

690

695

700

705

710

715

Fig. 16. PDF of total operation costs for the expected best compromise solution of S2.

11 10 9 8 7 6 5 4 3 2 1 0 690

695

700

705

710

715

720

725

730

735

Fig. 17. PDF of total operation costs for the expected best compromise solution of S3.

RIq ¼

Ns 1X ðfq;s ð½Ps ; U; fs Þ  fq;s ð½Ps ; Us ; fs ÞÞ; Ns s¼1

q ¼ 1; 2

ð48Þ

It is worthwhile to note that the smaller this criterion gets, the results of the proposed stochastic method are more close to the deterministic results. In the optimization process, each scenario considers different values for load, WT and PV powers output and market price. Also, each scenario has different decision variables but the same state variables. The expected operation costs and emissions of these 30 scenarios are the output random variables which should be optimized simultaneously in the algorithm procedure. This procedure is solved by the efficient ITLBO method and so the concluded stochastic POFs are extracted and shown in Figs. 9–11 for S1, S2, and S3, respectively. The total CPU times for these stochastic problems are about 0.512, 0.673, and 1.977 min for S1, S2, and S3, respectively. The PDF of the operation costs of the expected best compromise solutions for S1, S2, and S3 through scenarios are presented in Figs. 15 and 16, correspondingly. As can be seen from these figures, the output random variables e.g. operation costs tend to have the same PDF as the input random variables, which in this study is a normal distribution function. Therefore, a normal distribution is fitted to the scenario results using its mean and SD. However, it

469

T. Niknam et al. / Applied Energy 99 (2012) 455–470 Table 7 Results of different combinations of weight factors for SMOMGO (all strategies). Combination

wq

S1

S2

S3

0

w1 = 1.00

^f ¼ 274:6370 €ct 1 ^f ¼ 613:3658 kg 2

^f ¼ 293:5038 €ct 1 ^f ¼ 706:9867 kg 2

^f ¼ 308:5888 €ct 1 ^f ¼ 735:3178 kg 2

^f ¼ 361:2429 €ct 1 ^f ¼ 432:9583 kg 2

^f ¼ 388:4731 €ct 1 ^f ¼ 615:1857 kg 2

^f ¼ 407:2744 €ct 1 ^f ¼ 652:6287 kg 2

^f ¼ 597:8253 €ct 1 ^f ¼ 309:4458 kg 2

^f ¼ 697:1258 €ct 1 ^f ¼ 507:8045 kg 2

^f ¼ 713:2213 €ct 1 ^f ¼ 564:6447 kg 2

^f ¼ 983:0139 €ct 1 ^f ¼ 242:4756 kg 2

^f ¼ 1045:9978 €ct 1 ^f ¼ 438:4151 kg 2

^f ¼ 1072:3867 €ct 1 ^f ¼ 490:0273 kg 2

^f ¼ 1472:4627 €ct 1 ^f ¼ 192:1161 kg 2

^f ¼ $1410:0035 €ct 1 ^f ¼ 400:4604 kg 2

^f ¼ 1319:9998 €ct 1 ^f ¼ 440:1177 kg 2

w2 = 0.00 1

w1 = 0.80 w2 = 0.20

2

w1 = 0.50 w2 = 0.50

3

w1 = 0.20 w2 = 0.80

4

w1 = 0.00 w2 = 1.00

is well depicted from these figures that the discrete and nonsmooth behavior of the SMOMGO problem skews the PDF of the output random variables. This behavior is difficult to be predicted in advance, because it depends on the all of the input random variables uncertainties (see Fig. 17). The main goal of extracting the non-dominated solutions in the proposed SMOMGO problem during the search process is obtaining of the best compromise solution based on the need of the system operator. The MGCC tends to supply load with the lower cost as well as lower emission; therefore the MGCC should implement the proposed approach under uncertain environment to obtain the reliable solutions. The MGCC can select different combinations P of weight factors in such a way that 2q¼1 wq ¼ 1. The MGCC terminates the combination steps if the desired solution is found; otherwise, it selects a new wq. This procedure continues until the MGCC is satisfied with the solution according to the system operating conditions based on the objective functions. To inform the compromise solution values associated with different combinations of weight factors, Table 7 is offered. In this table, w1 and w2 are weight factors corresponding to the operation costs and emissions, respectively. In addition, the decision makers can bias preference through the search process perfectly according to the min–max approach which is described in (36) and presented in flowchart. This procedure is developed to decide on the best candidate solutions for the next optimization part of the proposed algorithm. In this work, the lref,qq = 1, 2 is set to one for both objective functions to assign the same importance to each of them. It can be inferred from the results that modeling of the system uncertainties will increase the operation cost and emission values because stochastic procedure considers different most probable scenarios instead of one scenario (as the deterministic scheme). Indeed, approaching to the real conditions of the power system in the ante-scheduling studies will cost some expenses which are expectable. In other words, SMOMGO will concurrently consider the most probable scenarios. Besides, using the proposed stochastic framework, all 30 accepted scenarios according to their probability values contribute into the output random variable results, whereas the deterministic method relies on only one scenario. The 30 accepted scenarios capture more of the uncertainty spectrum of the power system, which is approximately four times more than that of the deterministic framework. So, the MOMGO results of the stochastic framework are more realistic than the deterministic framework results. From the result of the case studies and strategies, it can be deduced that the proposed SMOMGO have benefit from the both aspects, i.e. presentation a flexible method to compromise the conflicting objectives, total operation cost and emissions minimization, with several constraints specially SRRs as well as more efficient handling of system uncertainties. The number of variables

Table 8 Optimization statistics for all strategies. Problem

No. of variables

No. of constraints

S1 DOMGO1&2 SOMGO1&2 DMOMGO SMOMGO

96 2880 96 2880

192 5760 192 5760

S2 DOMGO1&2 SOMGO1&2 DMOMGO SMOMGO

96 2880 96 2880

216 6480 216 6480

S3 DOMGO1&2 SOMGO1&2 DMOMGO SMOMGO

192 2976 192 2976

216 6480 216 6480

and constraints for the whole case studies and strategies are listed in Table 8. The advantages of the proposed framework have been investigated in the entire previous sections and concluded in Section 6. Beside, one of the major drawbacks of the proposed ITLBO method in comparison with the classical optimization method such as linear programming [11] and mixed-integer linear programming [17], is its execution burden that increases with the number of the DG and storage units which are used in the micro-grid network. However, it is noteworthy that these methods lead to non-optimal solution with great economic and emission loss. Although it is usually motivating to increase the computational performance of the method, but since the operation cost and emission are very important in optimal micro-grid operation, a reasonable increase in computation time would cause no severe problem. In order to overcome the drawback of mathematical methods in solution of optimal operation, the proposed ITLBO algorithm has been employed to solve the problem with practical modeling of operation constraints. It should be pointed out that another disadvantage of the meta-heuristic evolutionary algorithms is their less robustness in achieving optimal solution. The simulation results demonstrate the robustness and consistent of the proposed ITLBO approach 6. Conclusion In this paper, an ITLBO method has been developed and applied to solve the SMOMGO with clean sources such as WT, PV, and PAFC for economic and emission operation problem by considering uncertainties including WT and PV units powers output, load demand, and market price over the 24 h study horizon. Also for get-

470

T. Niknam et al. / Applied Energy 99 (2012) 455–470

ting closer to real condition as well as the reliability reasons, a reserve constraint is taken into account. The best advantage of the proposed algorithm is quick transfer of the information between agents which this gives more ability to the proposed algorithm in finding the global optima irrespective of the complexity of the problem. Besides, a SAPMS is complemented to the original TLBO to cope with the drawback of premature convergence. This modification includes two powerful knowledge interaction strategies. Each individual according to a probability model chooses one of these methods to improve its knowledge. Since some RESs such as WT, PV have intermittent characteristic, approaches to analyze MGs would be stochastic rather than deterministic. The proposed approach shows how the proposed formulation works in comparison with an unreal-case-based deterministic technique. To take the uncertainties into account, a two-stage scenario-based method according to the Monte Carlo technique is implemented. Using the LMCS rank-1, possible scenarios of power system operating states are generated and a probability is assigned to each scenario. For a tradeoff between computation time and accuracy, a backward scenario reduction technique is utilized. The main goal of the MOP is focusing on the minimization of the operation cost and emission level in MGs without sacrificing the realization of the system. Since the proposed algorithm profits the niching and fuzzy clustering methods so it can obtain a true and well-distributed set of Pareto-optimal solutions which gives the MGCC several chances to select an appropriate power dispatch plan according to environmental or economical considerations. The capability and performance of the proposed approach is tested on a typical MG including various DGs under three strategies. The PDF of the output random variables are achieved for each strategy. Also, a comparison between the deterministic and the stochastic approach in the term of the VSS and RI is carried out. As a result, this study provides a method for the MGCC with the blueprint for control and management of the challenging SMOMGO problem emerged from proliferation of RESs in MG systems

References [1] Llaria A, Curea O, Jiménez J, Camblong H. Survey on microgrids: unplanned islanding and related inverter control techniques. Renew Energy 2011;36:2052–61. [2] Hawkes AD, Leach MA. Modelling high level system design and unit commitment for a microgrid. Appl Energy 2009;86:1253–65. [3] Lasseter RH. MicroGrids. IEEE Power Eng Soc Winter Meet 2002;1:305–8. [4] Giannoulis ED, Haralambopoulos DA. Distributed Generation in an isolated grid: methodology of case study for Lesvos–Greece. Appl Energy 2011;88:2530–40. [5] Moghaddam AA, Seifi A, Niknam T. Multi-operation management of typical micro-grids using particle swarm optimization: a comparative study. Renew Sust Energy Rev 2012;16:1268–81. [6] Ren H, Zhou W, Nakagami K, 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. [7] Niknam T, Taheri SI, Aghaei J, Tabatabaei S, Nayeripour M. A modified honey bee mating optimization algorithm for multiobjective placement of renewable energy resources. Appl Energy 2011;88:4817–30. [8] Soroudi A, Ehsan M, Zareipour H. A practical eco-environmental distribution network planning model including fuel cells and non-renewable distributed energy resources. Renew Energy 2011;36:179–88. [9] Hessami MA, Bowly DR. Economic feasibility and optimisation of an energy storage system for Portland Wind Farm (Victoria, Australia). Appl Energy 2011;88:2755–63. [10] Kazempour SJ, Parsa Moghaddam M, Haghifam MR, Yousefi GR. Electric energy storage systems in a market-based economy: comparison of emerging and traditional technologies. Renew Energy 2009;34:2630–9.

[11] Chakraborty S, Weiss MD, Simoes MG. Distributed intelligent energy management system for a single-phase high frequency AC microgrid. IEEE Trans Ind Electron 2007;54:97–109. [12] Sortomme E, El-Sharkawi MA. Optimal power flow for a system of microgrids with controllable loads and battery storage. IEEE/PES Power Syst Conf Expos 2009:1–5. [13] Hernandez-Aramburo CA, Green TC, Mugniot N. Fuel consumption minimization of a microgrid. IEEE Trans Ind Appl 2005;41:673–81. [14] Mohamed FA, Koivo HN. System modelling and online optimal management of microgrid with battery storage. In: Proc 6th international conf renew energies power, quality; 2007. [15] Mohamed FA, Koivo HN. System modelling and online optimal management of microgrid using mesh adaptive direct search. Int J Electr Power Energy Syst 2010;32:398–407. [16] Logenthiran T, Srinivasan D. Short term generation scheduling of a microgrid. In: Tec IEEE region 10th conf; 2009. p. 1–6. [17] Morais H, Kàdàr P, Faria P, Vale ZA, Khodr HM. Optimal scheduling of a renewable micro-grid in an isolated load area using mixed-integer linear programming. Renew Energy 2010;35:151–6. [18] Dukpa A, Dugga I, Venkatesh B, Chang L. Optimal participation and risk mitigation of wind generators in an electricity market. IET Renew Power Gener 2010;4:165–75. [19] Akbari R, Hedayatzadeh R, Ziarati K, Hassanizadeh B. A multi-objective artificial bee colony algorithm. Swarm Evol Comput 2012;2:39–52. [20] Akbari R, Ziarati K. Multiobjective bee swarm optimization. Int J Innov Comput I (IJICIC) 2012;8:715–26. [21] Pablo Ruiz A, Russ Philbrick C, Peter Sauer W. Modeling approaches for computational cost reduction in stochastic unit commitment formulations. IEEE Trans Power Syst 2010;25:588–9. [22] Rao RV, Savsani VJ, Vakharia DP. Teaching–learning-based optimization: A novel method for constrained mechanical design optimization problems. Comput Aid Des 2011;43:303–15. [23] Agrawal S, Panigrahi BK, Tiwari MK. Multiobjective particle swarm algorithm with fuzzy clustering for electrical power dispatch. IEEE Trans Evol Comput 2008;12:529–41. [24] Wu L, Shahidehpour M, Li T. Stochastic security-constrained unit commitment. IEEE Trans Power Syst 2007;22:800–11. [25] Billinton R, Allan RN. Reliability evaluation of power systems. 2nd ed. New York: Plenum; 1996. [26] Michalewicz Z. Genetic algorithms + data structures = evolution programs. 3rd ed. New York: Springer; 1996. [27] Wu L, Shahidehpour M, Li T. Cost of reliability analysis based on stochastic unit commitment. IEEE Trans Power Syst 2008;23:1364–74. [28] Billinton R, Li W. Reliability assessment of electric power systems using Monte Carlo methods. New York: Plenum; 1994. [29] Siahkali H, Vakilian M. Stochastic unit commitment of wind farms integrated in power system. Electr Power Syst Res 2010;80:1006–17. [30] LanLan Z, Ling W, Xiuting W, Ziyuan H. A novel PSO-inspired probability-based binary optimization algorithm. Int Symp Inform Science Eng 2008:248–51. [31] Sakawa M, Yano H. An interactive fuzzy satisfying method using augmented minimax problems and its application to environmental systems. IEEE Trans Syst Man Cybern 1985;15:720–9. [32] Malekpour AR, Tabatabaei S, Niknam T. Probabilistic approach to multiobjective volt/var control of distribution system considering hybrid fuel cell and wind energy sources using improved shuffled frog leaping algorithm. Renew Energy 2012;39:228–40. [33] Azizipanah-Abarghooee R, Niknam T, Roosta A, Malekpour AR, Zare M. Probabilistic multiobjective wind-thermal economic emission dispatch based on point estimated method. Energy 2012;37:322–35. [34] Wong KP, Dong ZY. Differential evolution, an alternative approach to evolutionary algorithm. In: Proc 13th int conf intel syst appl power syst; 2005. p. 73–83. [35] Niknam T, Narimani MR, Jabbari M, Malekpour AR. A modified shuffle frog leaping algorithm for multi-objective optimal power flow. Energy 2011;36:6420–32. [36] Goldberg DE. Genetic algorithms for search. Optimization, and machine learning, reading. MA: Addison-Wesley; 1989. [37] Moghaddam AA, Seifi A, Niknam T, Pahlavani MR. Multi-objective operation management of a renewable MG (micro-grid) with back-up micro-turbine/fuel cell/battery hybrid power source. Energy 2011;36:6490–507. [38] Morales JM, Minguez R, Conejo AJ. A methodology to generate statistically dependent wind speed scenarios. Appl Energy 2010;87:843–55. [39] Breipohl A, Lee FN, Huang J, Feng Q. Sample size reduction in stochastic production simulation. IEEE Trans Power Syst 1990;5:984–92. [40] Morales JM. Impact on system economics and security of a high penetration of wind power. PhD dissertation, Dept Electr Eng, Univ Castilla-La Mancha, Ciudad Real, Spain; 2010.