WIXLIB

E. Shahryari, H. Shayeghi, B. Mohammadi-Ivatloo, M. Moradzadeh: Optimal Energy Management of Microgrid in Day-Ahead 88Ref. [13] has formulated the microgrid EMS as astochastic problem utilizing auto-regressive movingaverage (ARMA) to model existing uncertainties. Inaddition, authors of this paper have utilized a DRP thatmodifies load pattern by transferring demand from peakperiod to other time intervals. Application of both pricebased and incentive-based DR program in microgridEMS is studied in Ref. [14]. Moreover, the uncertaintiesrelated to price, load, wind speed and solar radiation ismodeled using probabilistic distribution function (PDF)of each parameter.pollution as the objective functions, the EM problem issolved using GSO algorithm.An information gap decision theory (IGDT) ispresented in Ref. [15] to obtain the bidding strategy ofmicrogrid as an stochastic problem. In Ref. [16], anadvanced dynamic programming method is presented tosolve microgrid EMS problem where probabilisticmodels for uncertain parameters are also considered. Ref.[17] has developed a linear problem for EMS, and MonteCarlo simulation has been applied to cover uncertaintiesof biomass generation of microgrid. In Ref. [18], a novelmulti-agent system (MAS) based model is presented tosolve the optimization problem. The contribution of Ref.[19] is to solve the multi-objective problem of EMS usingaugmented ɛ-constraint method. In Ref. [20], analternating direction method for multiplier (ADMM) hasbeen proposed to schedule the operation of microgridcomponents.The remainder of this paper is organized as follows: InSection 2 the proposed copula-based uncertaintymodeling technique is presented. Section 3 describes l formulation of objective functions andconstraints are presented in Section 4. The numericalresults and conclusions are expressed in Section 5 and 6,respectively.In this paper, the studied microgrid is grid-connected,and is composed of WT, PV, fuel cell (FC), micro-turbine(MT) and ESS. This paper presents a new hybriduncertainty modeling method to cover uncertaintiesassociated with wind and solar power. The presentedmethod is based on copula function and scenariogeneration. Firstly, the real and forecasted values ofuncertain parameters are utilized to generate a twodimensional conditional PDF using copulas. Then, aspecified number of scenarios are generated by scenariogeneration technique. Introducing the improvedincentive-based DRP is the second contribution of thispaper. Computing more realistic values for penalty andincentive payments as a function of peak intensity, takinginto account intra-day market which results in moreaccurate participation of consumers in DR, andconsidering the elasticity values as a function of costumertype and peak intensity are the properties of our proposedDR program. Considering the intra-day market inaddition to day-ahead market, in the formulation andsolving of the EM problem, is the third contribution ofthis paper. The latter allows microgrid components tomake more accurate participation in EM of microgrid.Considering the operational cost and environmentalIn summary, this paper has the following prominentfeatures: Proposing a copula-scenario based techniquefor uncertainty modeling Introducing an improved incentive-based DRP Considering both day-ahead and intra-daymarket in microgrid EM problemProposed copula-based uncertaintymodeling techniqueIt is clear that there is a stochastic dependency betweenvarious random variables. For example, the forecastedand real-time wind power generations have a significantlevel of correlation. As a result, more accurate predictionof wind power results in more precise planning of realtime wind power utilization. This explanation can beapplied for other stochastic variables such as solarradiation, electricity price and load. Various methodshave been employed in order to manage the uncertaintyof the aforementioned uncertain parameters in EM ofmicrogrid. The scenario generation method in Ref. [21]generates various scenarios for uncertain parameter usingthe PDF of prediction error. However, in this method, thePDF is made by prediction error values regardless ofpairing each error value to its real amount. The MonteCarlo method [22] is relied on repeated randomsamplings. The PDF function used in this method is thesame as PDF of scenario generation technique. Anotheruncertainty modeling method which is based on PDF ispoint estimate method [23]. Along with aforementionedstochastic methods which are based on PDF, there existsother techniques such as fuzzy method [24], robustoptimization [25], IGDT [26] and conditional value atrisk (CVaR) [27]. However, the proposed copula-scenariobased uncertainty modeling technique is a hybrid methodwhich has utilized the copula function to improve theaccuracy of the PDF. The generated PDF is more accuratesince it is made by pairing real value and its associatedprediction error. It is worth mentioning that the2.

Journal of Operation and Automation in Power Engineering, Vol. 8, No. 2, Aug. 2020equivalent weighted hourly value for each of uncertainparameters at the last step of the proposed method hasimproved the computation speed in comparison to theexisted stochastic methods.89cumulative distribution functions are FX and FYrespectively, the copula C can join their distributionfunction as below:FX Y (x , y ) C ( FX ( x ), FY ( y ))(3)2.1. Measurement of stochastic dependencyThe rank correlation (ρr), is utilized to compute theseverity of dependency between random variables [28].The rank correlation of random variables X and Y withcalmative distribution functions (CDFs) FX and FY is: r (X ,Y ) (FX (X ), FY (Y ))(1)Here, ρ(FX(X),FY(Y)) is called the product momentcorrelation of two CDFs which is computed as:(2)Cov (FX , FY ) (FX (X ), FY (Y )) (FX ) (FY )Where, Cov and σ denote the covariance and standarddeviation, respectively. It can be seen that the rankcorrelation is within interval [-1, 1] and is a factor thattransforms random variables to uniform randomvariables. If ρr 1, the dependency between randomvariables is excellent. In contrast ρr 0 shows theindependency of random variables.Assume FX(x) u and FY(y) v, where u and v arerealizations of the uniform random variables U and V,respectively. As a result, Eq. (3) can be written as:CVU(u ,v ) F (x , y ) F (FX 1 (u ), FY 1 (v ))Where,CV U(4)is the conditional distribution of V U andF-1 is the inverse of univariate distribution function.2.3. Our proposed methodIn this section the steps describing our proposeduncertainty modeling technique is expressed. Theproposed method utilizes predicted and real-time data asdependent random variables. The overall steps are asbelow:1) Collect the hourly historical data for the uncertainparameter for 1 year2) Predict hourly values of the parameter for one weekusing the collected historical data2.2. CopulaCopulas are mathematical functions that formulatemultivariate distribution functions by coupling multipleone-dimensional distribution functions [29]. This processis performed by transforming uniform distribution of allmarginal variables into a multivariate distribution [30].In mathematical terms, a copula is a function (C) with nvariables in the range of [0,1]n with the followingspecifications:1) The amount of C is in the unit interval [0,1];2) Considering u (u1, ,un) in [0,1] , C(u) is zero ifat least one coordinate is zero;n3) C(u) uk if all of coordinates equal to 1 except uk;According to Sklar’s theorem in Ref. [31] which is thebasis of the copula application, a copula function joinsmultiple one-dimensional distribution functions togetherto generate multivariate distribution functions. Sklar alsoproved that if the marginal distributions are continuous,there exists only one copula representation. There existstwo families of copulas; namely Ellipse family andArchimedean family. The commonly-used Gaussian andStudent t are included in Ellipse category while Frank,Gumbel and Clayton are contained in Archimedeanfamily [32].Considering X and Y as random variables whose3) Consider the predicted hourly values and the realtime amounts for the proportional time period as therandom variables4) Fit the data with copula function and generateconditional distribution function for real-time andpredicted values of uncertain parameter5) Generate joint conditional distribution estimationfor prediction error6) Generate specific number of scenarios using themethod presented in Ref. [33]7) Make weighted equivalent value for each hour ofparameter as below:nSEq . param (t ) Pr e . param (t ) Pr obS (t ) *ValueS (t )(5)S 1Here, Eq.param and Pre.param are the equivalentvalue and predicted value of uncertain parameter.Furthermore, ProbS and ValueS are the probability andvalue of parameter at Sth scenario.3. Improved incentive-based DR programIn this section, the mathematical formulation of theimproved incentive-based DR (IBDR) program which isan expansion of the proposed model in Ref. [34] ispresented. To begin with, it is crucial to define the

E. Shahryari, H. Shayeghi, B. Mohammadi-Ivatloo, M. Moradzadeh: Optimal Energy Management of Microgrid in Day-Ahead 90concept of elasticity. This coefficient is defined as theratio of demand modification at t 1 in response to relativeprice changes during t 2 as below:E (t1, t 2 ) 0, DA (t 2 )D 0 (t1 ). D (t1 ) DA (t 2 )(6)Where, E(t1,t2), D(t1) and ρ (t2) are the elasticitycoefficient, load change at t t1 and price modification att t2, respectively. The time-dependent nature of elasticityrelates the price changes at t t2 to the load variation att t1. While t1 t2, increment of price at t t2 will causeincrement of load at t t1. When t1 t2, decrement of priceat t t2 will result in increment of load at t t1 and viceversa. This is expressed in Eq. (7) as below: E (t1, t 2 ) 0 if t1 t 2 E (t1, t 2 ) 0 if t1 t 2(7)The aim of the IBDR program is to maximize the totalbenefit function which consists of incentive payment,penalty values and revenue.thIncentive payments are paid to the i consumer, whenit changes its demand from D0,i(t) to Di(t) as a result ofprice variation. The amount of incentive payment iscomputed as:inc i (t ) t ( peak density )[ D 0,i (t ) D i (t )](8)Here, αt,i represents the incentive payment ( ) per kWhat tth hour for ith consumer. This coefficient is a functionof peak intensity of tth hour, and results in more incentivepayment for load reduction during peak hours.The penalty value of consumers for breach of contractis:Microgrid energy management problemformulation4.The mathematical formulation of microgrid energymanagement problem considering objective functionsand constraints are discussed in this section.4.1. Objective functionsIn this paper, the aim of the EMS problem formulation isto find the optimal hourly planning for the microgridcomponents including generators (WT, PV, MT, FC) andESS while satisfying various technical constraints. Thecost and revenue of microgrid as a result of powerexchange by utility grid is included in mathematicalmodelling. In addition, in the paper the EMS of microgridin both day-ahead and intra-day markets are considered.Two objective functions including total operational costof microgrid and environmental pollution are to beminimized simultaneously.Total operational cost of microgrid4.1.1.The first objective function is total operational cost ofmicrogrid expressed by:24F (Cost ) (CWind (t ) C PV (t ) C MT (t ) C FC (t )(12)t 1 C Buy (t ) C Sell (t ) C ESS (t ))Where,TCWind (t ) COp WT .PWT (t ) CCons WT(13)t 1TC PV (t ) COp PV .PPV (t ) CCons PV(14)t 1pen i (t ) t ( peak density )[RLi (t ) [D 0,i (t ) D i (t )]](9)Where, RLi(t) and [D0,i(t)-Di(t)] demonstrate the realand promised load reduction, respectively. βt representsthe penalty value ( ) per kWh of contract violation.Considering rev(Di(t)) as the revenue of ith consumer,the total benefit function for ith consumer is:beni (t ) rev ( D i (t )) D i (t ). (t ) inc i (t ) peni (t )(10)According to classical optimization rules, the benefitfunction is maximized when its deviation equals to zero.After mathematical simplifications [35], the optimalconsumption of customer during 24 h is as below:24 24D i (t1 ) D 0,i (t1 ).exp( E i (t1, t 2 ).t1 1 t 2(15)24 C.P (t )C FC (t ) Fuel FC C Op FC .PFC (t ) C M FC FC t 1 (16)The costs related to WT, PV, MT, FC are explained inEqs. (13)-(16). According to Eqs. (13) and (14),operational costs of WT and PV are composed of variableand fixed terms. The first term of Eq. (15) represents thepower generation cost of MT while the second and thirdterms are the operational and maintenance costs,respectively. Similarly, the first term of Eq. (16) indicatesthe fuel cost of FC while the second and third terms arethe operational and maintenance costs, respectively.T (t1) t ( peak density ) t ( peak density ) 0 (t1) )1(11)24 C.P (t )C MT (t ) Fuel MT C Op MT .PMT (t ) C M MT MT t 1 1 0 (t1 )C Buy (t ) (t ).PBuy (t )t 1(17)

Journal of Operation and Automation in Power Engineering, Vol. 8, No. 2, Aug. 2020TC Sell (t ) (t ).PSell (t )(18)t 1TC ES (t ) C Op ES .PEs (t ) CCons ESS(19)t 1The costs related to power exchange by the upstreamgrid are described in Eqs. (17) and (18). Eq. (17)represents the costs of buying power while Eq. (18) is thecost of selling power. In addition, operational andmaintenance costs of ESS are given in Eq. (19).4.1.2.Pollution cost of microgridThe total emission cost of microgrid and utility grid(UG) is calculated as below:24F (Emission ) (E MT (t ) E FC (t ) EUG (t ))(20)t 1Where,24E MT (t ) PMT (t ).EFMT(21)t 124E FC (t ) PFC (t ).EFFC(22)t 124EUG (t ) PBuy (t ).EFUG(23)t 1Here, Eqs. (21)-(23) are the emitted pollution by MT,FC and UG.4.2. Constraints24equal to the total load of that period as given by the powerbalance Eq. (24). The limitations related to powergeneration of FC and MT are given in Eq. (25) and (26),respectively. The constraints on the charging-dischargingpowers of ESS is specified by Eq. (27). In this equation,state of charging and discharging are separated usingpositive and negative values, respectively. The thermalcapacity of line which connects microgrid to main-grid islimited by Eq. (28).5.In this section, numerical studies are presented. First, thestudied microgrid have been introduced. Then, thesimulation results of copula-based uncertainty modelingtechnique are presented. At this step, two-dimensionalprobability distribution function (PDF) of uncertainparameters including wind and solar power have beenextracted. Subsequently, the updated load patterns as aresult of utilizing incentive-based DRP in intra-day andday-ahead market are depicted. The last part ofsimulation results is assigned to the hourly powermanagement values of the microgrid components and themain grid. This section has presented hourly EMS valuesusing GSO in day-ahead and intra-day markets. Adetailed introduction about GSO algorithm has beenpresented in our previous work [7]. It is worthmentioning that the simulations are performed usingMATLAB R2013a, running on a PC with a 1.5 GHzAMD Quad core A4 CPU and 4GB RAM.(24)t 1The studied microgrid has been adopted from Ref. [36]and shown in Fig. 1. It consists of a WT, two PVs, FC,MT and an ESS system. It is clear that the microgrid isconnected to the utility grid. The values of emissionfactors of MT, FC and UG are given in Table 2. Thevalues of power limitations, efficiency values along withoperation and maintenance costs for microgridcomponents are given in Table 3.PFC (t ) PESS (t ) PBuy (t ) PSell (t ))maxPFC (t ) PFCNumerical Results5.1. Studied microgridMicrogrid EMS contains multiple constraints such aspower balance, and constraints associated with ES,power generation by DGs etc. Meeting these constraintscan result in reaching a feasible solution for EMSproblem. In this paper, the following constraints areconsidered for the operation of microgrids:load (t ) (PWT (t ) PPV (t ) PMT (t ) 91Utility GridUtility Grid(25)20kV/400VPMT (t ) maxPMT20kV/400V(26) PESS (t ) PEmaxfor disch (PESS (t ) 0) disch arg e D .P (t ) P maxfor ch (PESS (t ) 0)E ch arg e C ESS(27)PBuy (t ) or PSell (t ) PLine(28)The summation of the produced power by all members ofmicrogrid (DGs ESS) in each time interval must bePV (20kW)LoadLoadPV (20kW)LoadMT (30kW)WT (50kW)LoadFC (30kW)ES (30kW)Fig. 1. The studied microgrid [36]

E. Shahryari, H. Shayeghi, B. Mohammadi-Ivatloo, M. Moradzadeh: Optimal Energy Management of Microgrid in Day-Ahead 92Table 2. Emission factors (kg/MWh) [37]Emission 0752.295Table 3. Power limitation, efficiency factors and costs of microgridcomponents [37]DGtypeMTFCWTPVESSPower 0.30.30.590.95Maintenance( )Operation( .0045.2. Uncertainty modeling based on copula function402000102030Predicted wind power (kW)40(b)Fig. 3. 3-D view of the joint conditional distribution estimation of theactual and forecast values of (a) wind power, (b) solar powerBased on the conditional distribution functions of Fig. 3,the joint conditional distribution function for predictionerror of each parameter is extracted and shown in Fig 4.To do this, the prediction error for each pair of real-timeand predicted values are computed, and the probabilityamount are taken from conditional distribution values ofFig. 3.After generating the joint distribution function ofprediction errors, it is time to generate appropriatescenarios to cover uncertainties. Firstly, random valuesare generated and then by comparing the generatedrandom amounts by prediction error values of Fig. 4, theoccurrence probability for each scenario is extracted. Allof the generated scenarios for four uncertain parametersof this study are shown in Fig. 5.0.25ProbabilityReal-time wind power (kW)This section focuses on the simulation results of theproposed uncertainty modeling technique. As mentionedbefore, the presented method requires the real-time andpredicted values of uncertain parameter for one year timeduration. The data related to actual and forecasted valuesof wind and solar power are extracted from [38].Utilizing these data, the probabilistic relationshipbetween each pair of real-time and predicted hourlyvalues of wind power and solar power are shown byscatter plot in Fig. 2. The marginal distribution functionfor each of the actual and estimated values are alsoplotted by its proportional axis. According to Ref. [39],the Gaussian copula is more appropriate in modelingmultivariate PDFs. Considering this fact, this paper usesGaussian copula in the study case. Figure 3 demonstratesthe joint conditional distribution function for actual andforecasted values of uncertain parameters. The proposedresults approve logical conditional distributions whichare close to real PDFs.0.20.150.10.0500(a)0.050.10.15Wind power prediction error(kW)0.2(a)Probability0.8(b)Fig. 2. Scatter plot of the joint distribution of the actual andforecast values of (a) wind power, (b) solar power0.60.40.2000.050.10.150.2Solar power prediction error (kW)0.25(b)Fig. 4. Joint conditional distribution estimation of the predictionerror for (a) wind power, (b) solar power

Wind power (kW)Journal of Operation and Automation in Power Engineering, Vol. 8, No. 2, Aug. 20204030201001471013Hour1619222416192224(a)Solar power (kW)504030201001471013Hour(b)Fig. 5. Generated scenarios for (a) wind power, (b) solar powerWind power (kW)402010471013Hour1619222416192224(a)Solar power (kW)4030201001time value and its proportional prediction value.However, in PDF generation of the existed methods, theutilized data are not pair and it is not important that whichprediction is proportional to which real-time value.Because, the PDF of copula-scenario based method isgenerated by pairs of prediction-actual amounts, its PDFhas more accuracy than existed methods. Totally, thismethod has tried to improve the accuracy of the PDFfunctions using copulas. Making a weighted equivalentvalue for each hour of uncertain parameters can beconsidered as another prominent feature of this techniquesince this point can improve the computing speed incomparison to the existed stochastic methods. Along withabovementioned advantages, the proposed uncertaintymodeling technique is not a certain method and is basedon the probabilities. This is the first disadvantage of theproposed technique. Furthermore, requirement of a largeamount of historical hourly data can be considered asanother drawback of this method.5.3. Load patterns as a result of incentive-based DRP300193471013Hour(b)Fig. 6. Hourly equivalent amounts after uncertainty modelingtechnique for (a) wind power, (b) solar powerTo simplify the computation and speed up theoptimization process, a weighted equivalent scenario isgenerated using Eq. 5. The resulted equivalent scenariofor each uncertain parameter is shown in Fig. 6. Thegenerated scenario will cover effects of all predictionerror values and their proportional occurrenceprobability.Combining copula as a mathematical function withscenario generation method proposed a new hybriduncertainty modeling technique which was not presentedup to now. Application of the copula function inmodeling uncertainties is another prominent feature ofthis method. The copula generates 3-D PDF based on thescatter plot of the predicted and real-time values. It canbe seen that in the scatter plot each point is a pair of real-The numerical studies related to efficiency of ourproposed incentive-based DRP is presented in thissection. The amounts of self and cross elasticity valuesare extracted using fitting process of hourly historicaldata which is explained more in detail in [35]. Table 4provides the elasticity coefficient values in day-aheadand intra-day market. The positive values of this tablemean an inverse relation between price and demand, i.e.price increment results in demand reduction and viceversa. The opposite statement is true for the negativevalues.The equ

fuel cost ( ) C ESS costs of electrical energy storage ( ) C Op-ESS operation cost of ES ( ) C Cons-ESS constant cost of ES ( ) C Op-MT C Op-FC operation cost of MT and FC ( ) C M-MT C M-FC maintenance cost of MT and FC ( ) D 0,i (t) initial load of ith customer (kW) D(t1) load change at t t1 (kW) D i (t) load of ith customer after DR (kW) E .