# Probabilistic Operation of the Microgrids Including Active Loads and Distributed Generation Constrained to Flexibility and Reliability Indices

Authors: Hosein Hasan Shahi; Mehdi Nafar; Mohsen Simab
Journal: International Transactions on Electrical Energy Systems
 (2022)

---

## Abstract

This paper presents a flexible-reliable operation strategy for clean microgrids (MGs) consisting of power sources, energy storage systems (ESSs), and responsive loads. The proposed strategy attempts to provide the minimum expected operating cost for the MG, power sources, predicted pollutant emission, and expected energy not-supplied (EENS) because of anN − 1 event considering different applications. The objective function of the problem is subject to several constraints including AC optimal power flow equations, flexibility and reliability boundaries, and model of power sources and active loads. The Pareto optimization based on the weighted sum method is then incorporated to obtain a single-objective model. Probabilistic programming is also used to model the uncertainties of load, renewable energy, electrical energy price, and MG equipment accessibility. The ant-lion and crow search algorithms are merged to solve the problem and find a reliable optimal solution. A standard MG is employed to test the proposed strategy and indicate its capability in enhancing technical and economic indicators of the MG while ensuring that the MG is clean.


---

## Body

## 1. Introduction

With the developments in the application of green energy produced by distributed generation (DG), renewable distributed generation (RDG), energy storage systems (ESSs), and demand response programs (DRPs), local resources are a promising solution to supply most of the demand locally [1–5]. Such green resources will also help enhance technical and economic indicators by proper adoption of energy management systems (EMSs) [6–10] and considering the coordination between these power sources and distribution system operator (DSO) [11]. EMS can optimally manage the demand-side DGs and ALs, thus reducing the power loss and significant voltage variations on network buses, and resulting in enhanced operation of the network [12]. To reach higher-level flexibility for the network, flexibility resources like nonrenewable RDGs (NRDGs) and ALs can be utilized besides RDGs [13]. The adoption of local resources will also help reduce the rate of outages for customers in the case of internal faults in the equipment or external faults occurring due to natural disasters. Hence, proper EMS using DGs and ALs enhances the reliability of the network during an N − 1 event [14, 15].Energy management and optimal operation of distribution networks and MGs have widely been discussed in the literature. The stochastic programming and operation of a reconfigurable distribution network containing NRDGs was modeled and a network resilient against natural disasters was obtained [16]. The model included an objective function to minimize the resiliency cost (including repair and outage cost) and the planning cost (including capital cost of the network equipment and network operating cost). The aim was to enhance efficiency and profitability of the network. The same approach of Ref. [6] but using a stochastic-robust hybrid programming model was adopted in Ref. [17]. The results reported by Refs. [16, 17] show that the outstanding efficiency of 100% can be achieved by optimal allocation of sources and elements besides benefiting from an EMS for optimal reconfiguration of the distribution network. EMS implementation in a distribution network helps to provide flexibility by employing RGS such as batteries and electric vehicle (EV) parking lot [18] as these sources enjoy negligible time constant. They also provide a flexibility of 100%. DGs’ and EVs’ ability to manage electrical energy for the sake of enhancing the operation of a hybrid network was presented in Ref. [19], where DGs and EVs are used as energy hubs (EHs). By meeting the local demand, these sources help reduce power loss of and provide smoother profiles for voltage, pressure, and temperature [20]. Energy management is commonly modeled as nonlinear programming (NLP) or mixed-integer NLP (MINLP). A linear programming (LP) or mixed-integer LP (MILP) model was derived from the linear approximation method (LAM) [6–10] to find the optimal solution within the shortest possible time. However, as per Refs. [16–20], the computational error associated with calculating voltage, power, and loss variables is roughly 0.5%, 2.5%, and 10%, respectively. As a result, the LAM leads to high computational error, despite providing a rapid response. Operation and flexibility indicators in an MG consisting of RDGs, batteries, and electric springs were modeled in a power management problem [21]. The battery together with the system and electric spring managed to control reactive power and adjust the voltage of sensitive load buses while providing flexibility and helping in energy management.An efficient power management strategy was applied to a distribution network containing EVs parking lot [22] so that operation and power quality indicators were enhanced. Mathematical techniques like CONOPT and IPOPT [23] have also been used to solve the nonlinear energy/power management problem. The optimal solutions obtained in each solver are not the same because of the nonlinear and nonconvex nature of AC power flow (AC-PF) equations [21, 22]. An NLP model using nonhybrid evolutionary algorithms (NHEA) such as Antlion Optimizer (ALO), Differential Equations (DE), and Genetic Algorithm (GA) was discussed in Refs. [24–26]. The mentioned algorithms require many convergence iterations to find the optimal solution, so the computational time is considerable. Also, referring to statistical studies, such algorithms suffer from high standard deviation in the response; hence, the difference in the final solution of iterations is not negligible. Table 1 lists a summary of research works carried out so far. The authors in Ref. [27] present a model to preprocess simulation of large-scale converters in a DC microgrid. A current and voltage conditioning method based on distributed adaptive dynamic programming is used in Ref. [28] for a PV-wind system. Ref. [29] discusses a generative adversarial network based on trinetworks form (tnGAN). The network helps address the issues with leak detection caused by incomplete data of the sensors. Another leak detection method based on data is suggested in Ref. [30], where parameters of the pipelines, real flow variables, and previous values of pressure are adopted to present a pipeline model. To find the place of leak, a dynamic programming that depends on action and is subject to limitations of pressure and distance is incorporated. A robust optimization model of MGs is described in Ref. [31], where uncertain parameters are also considered to address the economical and robust operation of the MG. The model is structured in two stages and aims at balancing the economic and robust operation. The optimization helps find robust settable variables of the operation. The authors in Ref. [32] present a multi-stage stochastic programming model of MG operation, while the MG contains RDGs, storage, and thermal units. The model has two sub-divisions: short-term and long-term. The former adopts the predicted data during every six hours, and the latter deals with the energy storage value of the prediction period. Unit commitment of DGs using the PSO has also been addressed in the literature [33]. In [34], a stochastic optimization models various uncertain parameters. This model in its upper-level subproblem takes MGs as the leaders constrained by uncertain accessibility of PV and wind units and electrical energy price. Chance constraints help model the accessibility of RESs and thus evaluate the risk of excessive use of MG. Several scenarios are utilized to provide a model of electrical energy price. The lower-level subproblem deals with electrical energy dispatching for individual demand scenario. Still, there are some gaps in the research on energy management of distribution networks (or MGs) that should be addressed:(i)
Because of utilizing RDGs in power systems and considering the uncertainty associated with the output power of RDGs, the system flexibility becomes low. One common definition of flexibility is “Changing pathways of generation and consumption to react to an external price or activation signal, aiming to provide a service in the system” [21]. To realize flexibility, flexibility solutions including ALs and NRDGs can be incorporated in the system. This has been presented in many research pieces. However, the mathematical modeling of flexibility has not gained extensive attention [18, 21]. An index needs its related numerical results to be assessed accurately, which can be realized using mathematical modeling.(ii)
ESM of distribution networks or MGs, focusing on few indicators, has widely been discussed in the literature. Yet, improving mere one indicator cannot ensure improving other indicator(s). For instance, resiliency enhancement imposes high cost [16, 17]. So, several different indicators should be addressed by ESM at the same time so that the operation of the whole system is improved.(iii)
As mentioned, ESM problem is structured in the form of NLP or MINLP. Some research works like Refs. [16–20] adopt LAM to find the optimal solution of the problem, which demands low computational burden and time while leading to drastic errors. Mathematical approaches like CONOPT and IPOPT have also been incorporated [21, 22] to solve the EMS problem, even though the solutions provided are different. NHEA methods like GA have also been used to deal with the EMS, but the required time and the standard deviations in the responses are the challenges. Hybrid evolutionary algorithms (HEAs) help to reach the optimal solution with few iterations as several processes are performed in parallel while optimizing the given decision variables, leading to a low value of standard deviation.Table 1
A summary of the literature.
Ref.IndicesSolverFlexibility modelEconomicOperationEnvironmental emissionReliabilityFlexibility[16, 17]YesYesNoNoNoLAMNo[18]YesYesNoNoYesLAMYes[19]NoYesNoNoNoLAMNo[20]YesYesNoNoNoLAMNo[21]NoYesNoNoYesMAYes[22]NoYesNoNoNoMANo[24–26]NoYesNoNoNoNHEANo[27–30]NoYesNoNoNoMANo[31]YesYesNoNoNoLAMNo[32]NoYesNoNoNoMANo[33]YesYesNoNoNoNHEANo[34]NoYesNoNoNoMANoPMYesYesYesYesYesHEAYesPM: Proposed model.The present study attempts to cover and give solutions to the abovementioned research gaps by presenting EMS for MGs influenced by DGs and ALs. To this end, the study employs a flexible-reliable operation strategy (FROS) and takes into account technical indicators including operation status, flexibility, reliability, and economic and emission indicators. The proposed design establishes a three-objective problem to find the minimum expected operating cost of MGs and DGs, the estimated emission amount, and EENS for anN-1 event. By doing this, economic, environmental, and reliability indicators are modeled appropriately. The constraints are AC-PF equations of the MG, operating limits, flexibility, and reliability of the MG, and operational model of DGs and ALs, namely, DRP and battery storage. The weighted sum of functions is combined by the Pareto optimization to formulate the suggested design. Moreover, uncertainties associated with load, power production of RDGs, electricity price, and accessibility of MG devices are taken into account during internal fault occurrence. The uncertainties are modeled using the probabilistic programming. The roulette wheel mechanism (RWM) produces scenarios to be used in the design. Probability distribution function (PDF) relevant to the variables is established using the Gaussian method. An MINLP model is adopted for the problem and combined with ALO-CSA algorithm find the optimal solution with a low standard deviation. The innovations of the design include:(i)
Energy management of an MG penetrated with DGs and ALs is modeled to provide clean, reliable, and flexible energy for the network;(ii)
Optimal operation, flexibility, reliability, economic and environmental indicators are assessed in the MG at the same time;(iii)
The combined ALO-CSA is used to find the optimal solution with low standard deviation in the response; and(iv)
Probabilistic model of uncertainties is provided and Gaussian method is used to derive the exact PDF of variables.Also, the objectives of this paper are as follows:(i)
Simultaneous access to the optimal state of technical (operation, reliability, and flexibility), environmental, and economic indicators in the microgrid,(ii)
Deriving optimal power scheduling for sources, storage, and responsive loads according to the objectives of multi-criteria microgrids,(iii)
Extracting the almost unique optimal solution (with low dispersion in the final response) in a short computing time,(iv)
Determining the probability distribution function of problem variables.The layout of the paper is as follows. Problem formulation is given in Section2. The solution process is addressed in Section 3. The results are provided in detail in Section 4. In the end, Section 5 presents the conclusions.


## 2. Formulation of the Proposed Scheme

This section describes the model of energy management problem for an MG with DGs and ALs to improve economic, environmental pollution, reliability, operation, and flexibility conditions of the system based on the Pareto optimization technique. It is an optimization problem. The optimization model includes objective function and constraints [35–39]. The objective function minimizes the sum of the weighted functions of the expected operating cost of MGs and NRDGs, the expected pollution level, and the EENS resulting from the N -1 contingency. The scheme is also constrained to the AC optimal power flow (AC-OPF) equations, MG flexibility and reliability constraints, and the operation model of DGs and ALs. Thus, the problem will be formulated as follows:(1)minF=ωC∑w∈Wρw∑t∈Tλt,wPs,t,wU+∑n∈N∑l∈Lan,l+bn,lPn,l,t,wDG+cn,lPn,l,t,wDG2ืCost+ωE∑w∈Wρw∑t∈TγsPs,t,wU+∑n∈N∑l∈Lγn,lPn,l,t,wDGืEM+ωR∑w∈Wρw∑t∈T∑n∈NLn,t,wNSืEENS.Subject to:(2)ωC+ωE+ωR=1,(3)Ln,t,wNS+Pn,t,wU+Pn,t,wBdis−Pn,t,wBch+Pn,t,wDR+∑l∈LPn,l,t,wDG+∑j∈NIn,jPn,j,t,wF=Pn,t,wRL+Pn,t,wCL+Pn,t,wIL∀n,t,w,(4)Qn,t,wU+∑l∈LQn,l,t,wDG+∑j∈NIn,jQn,j,t,wF=Qn,t,wRL+Qn,t,wCL+Qn,t,wIL∀n,t,w,(5)Pn,j,t,wF=gn,jVn,t,w2−Vn,t,wVj,t,wgn,jcosφn,t,w−φj,t,w+bn,jsinφn,t,w−φj,t,w∀n,j,t,w,(6)Qn,j,t,wF=−bn,jVn,t,w2+Vn,t,wVj,t,wbn,jcosφn,t,w−φj,t,w−gn,jsinφn,t,w−φj,t,w∀n,j,t,w,(7)φs,t,w=0∀t,w,(8)Pn,j,t,wF2+Qn,j,t,wF2≤Sn,jFmax∀n,j,t,w,(9)Ps,t,wU2+Qs,t,wU2≤SsUmax∀t,w,(10)Vnmin≤Vn,t,w≤Vnmax∀n,t,w,(11)−εF≤Ps,t,wU−Ps,t,1U≤εF∀t,w,(12)0≤Ln,t,wNS≤Pn,t,wRL+Pn,t,wCL+Pn,t,wIL∀n,t,w,(13)Pn,l,t,wDG2+Qn,l,t,wDG2≤Sn,lDGmax∀n,l,t,w,(14)Pn,l,t,wDG=Pn,lDGmax∀n,l∈RES,t,w,(15)−ξnRLPn,t,wRL−ξnCLPn,t,wCL−ξnILPn,t,wIL≤Pn,t,wDR≤ξnRLPn,t,wRL+ξnCLPn,t,wCL+ξnILPn,t,wIL∀n,t,w,(16)∑t∈TPn,t,wDR=0∀n,w,(17)0≤Pn,t,wBdis≤DRnxn,tB∀n,t,w,(18)0≤Pn,t,wBch≤CRn1−xn,tB∀n,t,w,(19)Enmin≤Enini+∑τ=1tηnchPn,τ,wBch−1ηndisPn,τ,wBdis≤Enmax∀n,t,w.The objective function of the proposed problem is given in equation (1), which is based on the Pareto optimization technique based on the weighted functions sum method [40]. The first term in this equation expresses the minimization of the total expected operating cost of the MG and NRDGs, where the operating cost of the MG is equal to the cost of energy purchasing from the upstream network based on the energy market price [16]. Therefore, if the variable PU has a positive (negative) value, the first term of the Cost function in equation (1) represents the energy purchase cost (sales revenue) of the MG. In addition, since PU depends on the charging and discharging capacity of ALs, the operating cost of ALs is also estimated in the Cost function [41]. In the second tern of equation (1), the minimization of the expected pollution level resulting from receiving energy from the upstream network and NRDGs is presented. In this regard, the pollution coefficient, γ, is equal to the sum of the pollution coefficients resulting from NOX, SO2, and CO2 pollutants [42]. Next, the minimization of the EENS resulting from the N − 1 contingency during an internal failure of MG equipment is appropriate to the third term of equation (1) [43]. It is noteworthy that the mentioned objective functions are modeled as a weighted sum in an objective function F based on equation (1), in which the sum of weight coefficients, ωC, ωE, and ωR, should be 1 based on the constraint (2) [40].The AC-OPF constraints for the MG are expressed in equations (3)–(10) [21, 22]. Equations (3)–(7) represent the AC-PF equations, which model the active and reactive power balance in each bus, equations (3)–(4), active and reactive power flow through the distribution line, and the voltage phase angle of the slack bus [21]. Constraints (8)–(10) also introduce operation constraints of the MG, including the capacity limits of distribution lines and substations, equations (8) –(9), and the constraint of the bus voltage, Equation (10) [22]. In these equations, it is assumed that the MG is connected to the upstream network via a distribution substation connected to the slack bus, so the PU and QU variables for other buses are zero. The flexibility limit of the MG is further stated in constraint (11). Note that due to low flexibility in MG including RDG arising from forecasting error of RDG power generation, the resulting of the day-ahead and real-time scheduling of network is not the same [18]. Therefore, in a flexible MG, the distance between the active power of the distribution substation in scenario w should be small compared to the scenario corresponding to the deterministic programming (in this case, it is considered equal to the first scenario) [18]. As a result, the flexibility tolerance (εF) for the flexible MG is very low, where for 100% flexibility, its value will be zero. In the end, the limit of MG reliability is considered in equation (12), which refers to the limit of consumption load interrupted in the MG [16].The operation model of DGs and ALs is presented in constraints (13)–(19). Equations (13) –(14) concerns the performance of DGs in MG, where equation (13) refers to the capacity limit of these sources. Equation (14) also indicates the amount of power generation of RDGs, where PDGmax depends on climatic conditions such as solar irradiance and wind speed. Also, modeling of two types of ALs such as DRP and battery is presented in this problem. Constraints (15) and (16) are commensurate with the performance of the incentive-based DRP [43]. In this type of DRP, it is assumed that residential, commercial, and industrial consumers shift part of their energy consumption from peak load hours (corresponding to high energy prices) to nonpeak load hours (corresponding to low energy prices). Hence, it is expected that consumers participating in the proposed DRP will be able to reduce their energy costs, which is a good incentive for them. Thus, constraint (15) indicates the power variation limitation of the DRP, and equation (16) ensures that the reduced energy consumption of consumers during peak hours at off-peak hours is provided by the MG and various sources [43]. Finally, the formulation of the battery is mentioned in equations (17)–(19) [14], which indicates the limitations of discharge rate, charge rate, and energy storage in the battery, respectively. In these equations, the binary variable xB represents the battery charge and discharge performance, which prevents simultaneous battery charge/discharge operation [14]. Finally, it should be noted that the formulation presented in this section is suitable for AC microgrids of different sizes.In this case, parameters such as the active and reactive load of residential,PRL and QRL; commercial, PCL and QCL; and industrial consumers, PIL and QIL; energy price, λ; power generation of RDG sources, PDGmax; and availability of MG equipment such as lines and distribution substations during earthquakes are uncertain. Therefore, in this paper, probabilistic programming is used to model these parameters. To this end, the RWM produces a certain number of scenarios. In each scenario, the probabilities of load values and energy prices are calculated from the normal PDF, and the probability of the PDGmax value for the wind (photovoltaic) system is obtained from (Beta) Weibull PDF [1, 44, 45]. Finally, the probability uncertainty of the availability of MG equipment is based on Bernoulli PDF [43]. The Bernoulli PDF is based on the forced outage rate (FOR) of equipment against internal failures, and the normal, Weibull, and Beta PDFs are used for an uncertainty parameter based on its mean and standard deviation. Then, after solving the problem described by equations (1)–(19), the Gaussian method (GM) is used to obtain the standard PDF of the problem variables. In this method, for accurate PDF estimation, the probability function of a variable is equal to the sum of weighted normal PDFs [46].


## 3. Solution Process

### 3.1. Fuzzy Decision-Making Method

The problem presented in the previous section follows the Pareto optimization technique based on the sum of the weighted functions. In other words, for different values of the weight coefficients ،ωC, ωE, and ωR, different values are obtained for the Cost, EM, and EENS functions. Depicting the obtained points in 3D coordinates represents the Pareto front of the proposed scheme. Hence, there is a need to achieve an optimal point in this situation, which according to Ref. [47] is known as the best compromise solution. Therefore, in this paper, the fuzzy decision-making method is used to reach this goal. The implementation steps of this method are provided as follows [47]:Step 1. Determining the membership function: For each of the Cost, EM, and EENS functions, the value of the linear membership function (f^) is determined for different values of coefficients ωC, ωE, and ωR with the following equation:(20)f^m=1fm≤fmminfm−fmmaxfmmin−fmmaxfmmin≤fm≤fmmax0fmmax≤fmm=Cost,EM,EENS.In equation (20), the terms fmin and fmax denote the minimum and maximum values of an objective function, respectively. To calculate these parameters, the problem described by equations (1)–(19) is solved for the three cases ωC = 1, ωE = 1, and ωR = 1 [47]. More details will be discussed in Subsection A in Section 5.1Step 2. Determining the minimum membership function: In this step, the minimum value of the set f^Cost,f^EM,f^EENS is calculated for different values of ωC, ωE, and ωR [47]. It is worth mentioning that the value obtained in this section is considered as the term fi, where i represents the step of changes in the mentioned weight coefficients.Step 3. Determining the best compromise solutioI: The best compromise solution will be between the Cost, EM, and EENS functions corresponding to the maximum value of f [47].


### 3.2. Hybrid Meta-Heuristic Algorithm

The proposed problem by equations (1)–(19) is an MINLP problem. Therefore, in this paper, to compensate for the third research gap given in Section 1, the hybrid ALO [24] and CSA [48] (ALO + CSA) solver is used to achieve the optimal solution. Since decision variables in this algorithm are updated in two general processes, namely, ALO and CSA phases, it is expected that it can be a reliable optimal solution with a low standard deviation in the response. Further details on the capabilities of the algorithm are discussed in Section 4.2.A.To solve the proposed problem by evolutionary algorithms, even using the mentioned algorithm, the problem variables are divided into two general categories. The first category refers to the decision variables, which include variables of thePDG for the set of L - RES, QDG, PDR, xB, PBdis, PBch, and LNS, the values of which are determined using the hybrid ALO + CSA algorithm proportional to equations (21)–(27). The other category has dependent variables that include PDG for the set of RES, PU, QU, PF, QF, V, and φ. The value of PDG for the RES set is calculated based on the constraint (14), and the other dependent variables are determined by the AC-PF equations (3)–(7). In this paper, the backward-forward power flow method is used to solve the AC-PF problem [49].(21)Pn,l,t,wDG∈0,Sn,lDGmax∀n,l∈L−RES,t,w,(22)Qn,l,t,wDG∈0,Sn,lDGmax∀n,l,t,w,(23)Pn,t,wDR∈∀n,t,w,(24)xn,tB∈0,1∀n,t,(25)Pn,t,wBdis∈∀n,t,w,(26)Pn,t,wBch∈∀n,t,w,(27)Ln,t,wNS∈∀n,t,w.The penalty function method is employed in this study to estimate the limits of operation, flexibility, and reliability of the MG, (8)–(12); capacity limit of the DG, (13); DRP constraint, (16); and the limit of energy stored in the battery, (19). In this method, the penalty function for the constraints a ≤ b and a = b is μ.max (0, a − b) and κ.(a − b), which are added to the main objective function, i.e., to equation (1), for the mentioned constraints as given in equation (28) [25]. μ ≥ 0 and κ ∈ (−∞, +∞) represent Lagrange multipliers. Finally, the function modeled in equation (28) is known as the fitness function [50].(28)minF=ωC×Cost+ωE×EM+ωR×EENS+∑n,j,t,wμn,j,t,wslmax0,Pn,j,t,wF2+Qn,j,t,wF2−Sn,jFmax+∑t,wμt,wsumax0,Ps,t,wU2+Qs,t,wU2−SsUmax+∑n,t,wμ¯n,t,wvmax0,Vn,t,w−Vnmax+μ¯n,t,wvmax0,Vnmin−Vn,t,w+∑t,wμ¯t,wpumax0,Ps,t,wU−Ps,t,1U−εF+μ¯t,wpumax0,−εF−Ps,t,wU+Ps,t,1U+∑n,l,t,wμn,l,t,wsdmax0,Pn,l,t,wDG2+Qn,l,t,wDG2−Sn,lDGmax+∑n,t,wμ¯n,t,wemax0,Enmin−Enini−∑τ=1tηnchPn,τ,wBch−1ηndisPn,τ,wBdis+∑n,wκn,wdr∑t∈TPn,t,wDR+∑n,t,wμ¯n,t,wemax0,Enini+∑τ=1tηnchPn,τ,wBch−1ηndisPn,τ,wBdis−Enmax.The solving process of the hybrid ALO + CSA algorithm for the proposed scheme is such that, initially, this solver determinesN random values (N is the population size) based on constraints (21)–(27) for the decision variables and the Lagrangian multipliers. Then, the values of the dependent variables are calculated from constraint (14) and the backward-forward power flow method. Next, the fitness function is determined for N values using the value of decision variables based on equation (28). This step is known as the initialization step. In the next steps, the decision variables are updated using the proposed algorithm based on the optimal value of the fitness function in the previous step, so that first the ALO process and then the CSA process are carried out. In these steps, the calculation of dependent variables and fitness function is based on the technique used in the initialization step. Finally, in this paper, it is assumed that the convergence conditions are obtained after repeating the updating steps of decision variables to a certain number called itermax. Finally, Figure 1 shows the flowchart of solving the proposed problem.Figure 1
Flowchart of the solution process.


## 3.1. Fuzzy Decision-Making Method

The problem presented in the previous section follows the Pareto optimization technique based on the sum of the weighted functions. In other words, for different values of the weight coefficients ،ωC, ωE, and ωR, different values are obtained for the Cost, EM, and EENS functions. Depicting the obtained points in 3D coordinates represents the Pareto front of the proposed scheme. Hence, there is a need to achieve an optimal point in this situation, which according to Ref. [47] is known as the best compromise solution. Therefore, in this paper, the fuzzy decision-making method is used to reach this goal. The implementation steps of this method are provided as follows [47]:Step 1. Determining the membership function: For each of the Cost, EM, and EENS functions, the value of the linear membership function (f^) is determined for different values of coefficients ωC, ωE, and ωR with the following equation:(20)f^m=1fm≤fmminfm−fmmaxfmmin−fmmaxfmmin≤fm≤fmmax0fmmax≤fmm=Cost,EM,EENS.In equation (20), the terms fmin and fmax denote the minimum and maximum values of an objective function, respectively. To calculate these parameters, the problem described by equations (1)–(19) is solved for the three cases ωC = 1, ωE = 1, and ωR = 1 [47]. More details will be discussed in Subsection A in Section 5.1Step 2. Determining the minimum membership function: In this step, the minimum value of the set f^Cost,f^EM,f^EENS is calculated for different values of ωC, ωE, and ωR [47]. It is worth mentioning that the value obtained in this section is considered as the term fi, where i represents the step of changes in the mentioned weight coefficients.Step 3. Determining the best compromise solutioI: The best compromise solution will be between the Cost, EM, and EENS functions corresponding to the maximum value of f [47].


## 3.2. Hybrid Meta-Heuristic Algorithm

The proposed problem by equations (1)–(19) is an MINLP problem. Therefore, in this paper, to compensate for the third research gap given in Section 1, the hybrid ALO [24] and CSA [48] (ALO + CSA) solver is used to achieve the optimal solution. Since decision variables in this algorithm are updated in two general processes, namely, ALO and CSA phases, it is expected that it can be a reliable optimal solution with a low standard deviation in the response. Further details on the capabilities of the algorithm are discussed in Section 4.2.A.To solve the proposed problem by evolutionary algorithms, even using the mentioned algorithm, the problem variables are divided into two general categories. The first category refers to the decision variables, which include variables of thePDG for the set of L - RES, QDG, PDR, xB, PBdis, PBch, and LNS, the values of which are determined using the hybrid ALO + CSA algorithm proportional to equations (21)–(27). The other category has dependent variables that include PDG for the set of RES, PU, QU, PF, QF, V, and φ. The value of PDG for the RES set is calculated based on the constraint (14), and the other dependent variables are determined by the AC-PF equations (3)–(7). In this paper, the backward-forward power flow method is used to solve the AC-PF problem [49].(21)Pn,l,t,wDG∈0,Sn,lDGmax∀n,l∈L−RES,t,w,(22)Qn,l,t,wDG∈0,Sn,lDGmax∀n,l,t,w,(23)Pn,t,wDR∈∀n,t,w,(24)xn,tB∈0,1∀n,t,(25)Pn,t,wBdis∈∀n,t,w,(26)Pn,t,wBch∈∀n,t,w,(27)Ln,t,wNS∈∀n,t,w.The penalty function method is employed in this study to estimate the limits of operation, flexibility, and reliability of the MG, (8)–(12); capacity limit of the DG, (13); DRP constraint, (16); and the limit of energy stored in the battery, (19). In this method, the penalty function for the constraints a ≤ b and a = b is μ.max (0, a − b) and κ.(a − b), which are added to the main objective function, i.e., to equation (1), for the mentioned constraints as given in equation (28) [25]. μ ≥ 0 and κ ∈ (−∞, +∞) represent Lagrange multipliers. Finally, the function modeled in equation (28) is known as the fitness function [50].(28)minF=ωC×Cost+ωE×EM+ωR×EENS+∑n,j,t,wμn,j,t,wslmax0,Pn,j,t,wF2+Qn,j,t,wF2−Sn,jFmax+∑t,wμt,wsumax0,Ps,t,wU2+Qs,t,wU2−SsUmax+∑n,t,wμ¯n,t,wvmax0,Vn,t,w−Vnmax+μ¯n,t,wvmax0,Vnmin−Vn,t,w+∑t,wμ¯t,wpumax0,Ps,t,wU−Ps,t,1U−εF+μ¯t,wpumax0,−εF−Ps,t,wU+Ps,t,1U+∑n,l,t,wμn,l,t,wsdmax0,Pn,l,t,wDG2+Qn,l,t,wDG2−Sn,lDGmax+∑n,t,wμ¯n,t,wemax0,Enmin−Enini−∑τ=1tηnchPn,τ,wBch−1ηndisPn,τ,wBdis+∑n,wκn,wdr∑t∈TPn,t,wDR+∑n,t,wμ¯n,t,wemax0,Enini+∑τ=1tηnchPn,τ,wBch−1ηndisPn,τ,wBdis−Enmax.The solving process of the hybrid ALO + CSA algorithm for the proposed scheme is such that, initially, this solver determinesN random values (N is the population size) based on constraints (21)–(27) for the decision variables and the Lagrangian multipliers. Then, the values of the dependent variables are calculated from constraint (14) and the backward-forward power flow method. Next, the fitness function is determined for N values using the value of decision variables based on equation (28). This step is known as the initialization step. In the next steps, the decision variables are updated using the proposed algorithm based on the optimal value of the fitness function in the previous step, so that first the ALO process and then the CSA process are carried out. In these steps, the calculation of dependent variables and fitness function is based on the technique used in the initialization step. Finally, in this paper, it is assumed that the convergence conditions are obtained after repeating the updating steps of decision variables to a certain number called itermax. Finally, Figure 1 shows the flowchart of solving the proposed problem.Figure 1
Flowchart of the solution process.


## 4. Numerical Results

### 4.1. Case Studies

The proposed scheme is implemented on a 69-bus radial MG as shown in Figure2 [51]. This network has base power and voltage of 1 MVA and 12.66 kV. The characteristics of distribution lines and substations in addition to peak load data are extracted from Ref. [51]. The allowable voltage range is [0.9, 1.1] p.u [52–56]. In this paper, it is assumed that the network has three types of consumers: residential, commercial, and industrial. Industrial consumers are at buses 49, 50, and 61, and commercial consumers are at buses 11, 12, 21, and 64. Other buses have residential consumers only. The daily load profile is equal to the product of the peak load and the daily load factor curve, which is plotted in Figure 3(a) for the mentioned consumers [57]. Also, the expected daily curve of energy prices will be as shown in Figure 3(b). In addition, it is assumed that each of the industrial consumers in buses 49 and 50 has a fuel cell (FC)-type DG with a capacity of 0.3 MVA. The industrial consumer in bus 61 has two DGs of FC type and microturbine (MT) with a capacity of 0.8 MVA and 0.7 MVA, respectively. Every commercial consumer needs a photovoltaic (PV) DG with a capacity of 0.2 MVA. Finally, two wind system (WS) DGs with a capacity of 0.3 MVA are located on buses 24 and 63. Note that the daily active power profile of RDG is equal to the product of its capacity and its daily power rate curve, which is shown for PV and WS in Figure 3(c) [58]. Coefficients of cost functions of fuel cost and emissions for different types of DGs are given in Table 2.Figure 2
The 69-bus Mg [51].Figure 3
Expected daily curve of (a) load factor [57], (b) energy price, and (c) RES power rate [58].
(a)(b)(c)Table 2
DGs characteristics.
DGTypeCoefficients of fuel cost functionPollution coefficients (kg/MWh)a ($)b ($/MWh)c ($/MWh2)NOxCO2SO21FC150180.010.00210.0003105.262MT200210.020.19950.0036723.933PV0000004WS000000This study assumes that the participation rate of commercial and residential consumers in the proposed DRP is 40%, but industrial loads have a participation rate of 30%. The network also has three batteries placed in buses 10, 24, and 63 with a capacity of 2.5 MWh and a charge/discharge efficiency of 95%. The battery charge/discharge rate is 0.5 MW and the minimum energy and initial energy are 0.25 MWh and 0.25 MWh, respectively. Furthermore, emission coefficients of NOX, SO2, and CO2 gases proportional to the energy received from the upstream network for the MG based on Ref. [42] are equal to 2.295 kg/MWh, 921.25 kg/MWh, and 3.583 kg/MWh, respectively. It is assumed that the FOR of the MG equipment is 1%. Finally, RWM generates 100 scenarios, in which the standard deviation for uncertainties of the load, energy cost, and active power of RDGs is set 10%. Moreover, to achieve high flexibility, εF is considered 0.05 p.u.


## 4.1. Case Studies

The proposed scheme is implemented on a 69-bus radial MG as shown in Figure2 [51]. This network has base power and voltage of 1 MVA and 12.66 kV. The characteristics of distribution lines and substations in addition to peak load data are extracted from Ref. [51]. The allowable voltage range is [0.9, 1.1] p.u [52–56]. In this paper, it is assumed that the network has three types of consumers: residential, commercial, and industrial. Industrial consumers are at buses 49, 50, and 61, and commercial consumers are at buses 11, 12, 21, and 64. Other buses have residential consumers only. The daily load profile is equal to the product of the peak load and the daily load factor curve, which is plotted in Figure 3(a) for the mentioned consumers [57]. Also, the expected daily curve of energy prices will be as shown in Figure 3(b). In addition, it is assumed that each of the industrial consumers in buses 49 and 50 has a fuel cell (FC)-type DG with a capacity of 0.3 MVA. The industrial consumer in bus 61 has two DGs of FC type and microturbine (MT) with a capacity of 0.8 MVA and 0.7 MVA, respectively. Every commercial consumer needs a photovoltaic (PV) DG with a capacity of 0.2 MVA. Finally, two wind system (WS) DGs with a capacity of 0.3 MVA are located on buses 24 and 63. Note that the daily active power profile of RDG is equal to the product of its capacity and its daily power rate curve, which is shown for PV and WS in Figure 3(c) [58]. Coefficients of cost functions of fuel cost and emissions for different types of DGs are given in Table 2.Figure 2
The 69-bus Mg [51].Figure 3
Expected daily curve of (a) load factor [57], (b) energy price, and (c) RES power rate [58].
(a)(b)(c)Table 2
DGs characteristics.
DGTypeCoefficients of fuel cost functionPollution coefficients (kg/MWh)a ($)b ($/MWh)c ($/MWh2)NOxCO2SO21FC150180.010.00210.0003105.262MT200210.020.19950.0036723.933PV0000004WS000000This study assumes that the participation rate of commercial and residential consumers in the proposed DRP is 40%, but industrial loads have a participation rate of 30%. The network also has three batteries placed in buses 10, 24, and 63 with a capacity of 2.5 MWh and a charge/discharge efficiency of 95%. The battery charge/discharge rate is 0.5 MW and the minimum energy and initial energy are 0.25 MWh and 0.25 MWh, respectively. Furthermore, emission coefficients of NOX, SO2, and CO2 gases proportional to the energy received from the upstream network for the MG based on Ref. [42] are equal to 2.295 kg/MWh, 921.25 kg/MWh, and 3.583 kg/MWh, respectively. It is assumed that the FOR of the MG equipment is 1%. Finally, RWM generates 100 scenarios, in which the standard deviation for uncertainties of the load, energy cost, and active power of RDGs is set 10%. Moreover, to achieve high flexibility, εF is considered 0.05 p.u.


## 5. Results

The proposed scheme along with the solution process is coded in MATLAB software, then the obtained numerical results are reported below.


### 5.1. Examining the Feasibility of the Proposed Solution Process

Table3 provides the results of the Pareto front for the proposed problem described by equations (1)–(19). In this table, four values of zero, 0.33, 0.5, and 1 are considered for each weighting coefficient ωC, ωE, and ωR. In the second to fourth rows of this table, the values of the objective functions Cost, EM, and EENS for the case studies are ωC = 1, ωE = 1, and ωR = 1, respectively. For ωC = 1, the minimum value of Cost ($ 873.26) is calculated because in this case, only the Cost function appears in the objective function (1). These conditions are available for EM and EENS for ωE = 1 and ωR = 1, respectively, where their minimum values are 47793 kg and 2.33 MWh, respectively. Also, the maximum values of the mentioned functions can be seen in these lines, so that the maximum values of Cost, EM, and EENS are $ 943.31, 55633 kg, and 8.67 MWh, respectively. Note that according to this table, it is observed that the direction of increase and decrease of one function is not in the same direction of other functions. For instance, the decrease of EENS is proportional to the increase in EM and Cost. The same is true for other functions. Eventually, to reach a compromise between these functions, this paper employs the fuzzy decision-making method. To this end, Table 4 reports the results of the best compromise solution obtained from various Pareto optimization techniques. In other words, to evaluate the Pareto technique, the results of other methods such as ε-constraint [47], the sum of normalized functions [41], and nondominated sorting genetic algorithm (NSGA) [59] are also given in this table. According to Table 4, the weighted functions sum method, described in Section 3.1, succeeds to obtain the minimum value for the three objective functions, namely, Cost, EM, and EENS. The ε-constraint method also provides the best compromise solution close to the results of the method presented in this paper. However, the best compromise point obtained from the other two methods has a significant difference compared to the method of weighted functions sum. The compromise point obtained in the proposed method is close to the minimum value of each objective function so that the distance between Cost at this point compared to its minimum value is approximately 4.87% ((879.84–873.26)/873.26). This value for EM and EENS power is roughly 0.75% and 16.7%, respectively.Table 3
The Prato front results for the proposed scheme.
ωCωEωRCost ($)EM (kg)EENS (MWh)100873.26556337.91010943.31477938.67001937.84534292.330.50.50890.23503628.060.500.5882.72540123.7800.50.5928.45489753.490.330.330.33895.68505533.61Table 4
The best compromise solution of the proposed scheme.
MethodCost ($)EM (kg)EENS (MWh)Sum of weighted function87984481532.72ε-constraint [47]88072482082.79Normalized objective function [41]88824490333.08Nondominated sorting genetic algorithm [59]88455487212.93Figure4 depicts the convergence results of the problem using the solvers of the DE [25], Teaching-Learning-based Optimization (TLBO) [60], CSA, ALO, and the suggested hybrid ALO + CSA. In these algorithms, the population size and maximum convergence iterations are 50 and 1000, respectively, and other adjustment parameters related to each solver are selected from Refs. [24, 25, 48, 60]. In addition, the problem is solved 20 times by each algorithm to calculate statistical indices such as standard deviation (SD) of the response. Finally, based on this figure, it can be seen that the proposed ALO + CSA algorithm obtains the minimum value for the objective function F, equation (1), compared to the other mentioned algorithms. Thus, according to Figure 4, it finds the best compromise point in the fewest number of convergence iterations (CI), which is 608. This convergence iteration corresponds to the minimum possible computational time (CT), which is 212.3 s. Other solvers require a CT higher than 225 s to achieve the best compromise solution in the proposed scheme, equations (1)–(19). In addition, the hybrid ALO + CSA algorithm has an SD of 0.97%, but this value for other algorithms is greater than 1.35%. This means that the proposed solver has a low dispersion in the final response relative to the other solvers, and has almost a unique response condition. Figure 5 shows the PDF of the objective functions Cost, EM, and EENS. In this figure, the PDF of each function is presented for RWM results, and also, the standard PDF determinates using the Gaussian method (GM). As can be seen from this figure, the minimum value of the Cost function reported in Table 4, which is $ 879.84, has a probability of about 8%. Also, the PDF of Cost is a nonnormal distribution, while the PDFs of EM and EENS have approximately normal distributions. Moreover, the minimum values of EM, i.e., 48153 kg, and EENS, i.e., 2.72 MWh, have a probability of 7% and 9%, respectively, according to Figure 5.Figure 4
Convergence curve of the proposed problem based on the different solvers.Figure 5
PDF of (a) Cost, (b) EM, and (c) EENS.
(a)(b)(c)


### 5.2. Evaluating the Performance of DGs and ALs in the MG

The expected daily curve of active and reactive power of DGs is presented in Figure6. Comparing Figures 3(c) and 6(a) and according to the data of Section 4.1, it is observed that PVs and WSs in all simulation hours inject active power equal to their maximum capacity into the MG proportional to climatic conditions. According to Table 2, their pollution coefficients and operating costs are zero; it is expected that they will inject high active power into the MG to minimize EM and Cost. Also, the minimization of the EENS is in proportion to the fact that local sources are responsible for supplying consumers [23]. Therefore, to minimize the EENS during N-1 contingency, RDGs need to inject high power into the MG. Besides, at hours 1:00 to 6:00, as the price of energy purchased from the upstream network, Figure 3(b), is lower than the fuel price of FCs and MTs, as shown in Table 2, NRDGs inject low active power into the MG during these hours to minimize environmental pollution and earthquake-induced shutdown. Nevertheless, at other times when the fuel cost of NRDGs is less than the energy price, they inject active power equal to their maximum capacity into the network. This performance of FCs and MTs from 7:00 to 24:00 is commensurate with the minimization of Cost, EM, and EENS functions. In addition, it can be seen from Figure 6(b) that DGs have allocated part of it to inject reactive power into the MG at times they have free capacity. FCs and MTs inject about 0.6 and 0.4 p.u. Reactive power into the MG during 1:00–6:00. At these hours, according to Figure 7, ALs perform charging operations and receive active power from the MG. Therefore, it is expected that there will be a high voltage drop in the MG during these hours. To compensate for this, NRDGs inject significant reactive power into the network. However, RDGs generally inject reactive power into the MG during most simulation hours. This performance is proportional to compensating for voltage drop caused by the consumption of ALs and passive loads in the MG.Figure 6
Expected daily curve of (a) DGs active power and (b) DGs reactive power.
(a)(b)Figure 7
Expected daily curve of (a) MG load and (b) battery active power.
(a)(b)Figure7 illustrates the expected daily performance curves of ALs such as DRPs and batteries. According to Figure 7(b), it can be observed that all residential, commercial, and industrial loads that have participated in the proposed DRP will increase their consumption between 1:00–7:00 and 23:00–00:00 because the energy price according to Figure 3(b) has a low value during these hours. Yet, in peak hours when the energy price according to Figure 3(b) has a high value of 30 $/MWh, the mentioned consumers operate in the low-consumption mode. This demand-side management approach is proportional to minimizing all three objective functions of operating costs, environmental pollution, and MG shutdown due to N – 1 contingency. In addition, battery performance is similar to that of consumers that participate in the proposed DRP so that they are charged during 1:00–6:00 and 23:00–00:00 to receive the energy needed for discharge during peak hours. Finally, it is noteworthy that the flexibility sources, including ALs and NRDGs, have a certain amount of active power at all simulation hours to improve network flexibility. This is because these sources must compensate for the prediction error in the power generation of RDGs to improve system flexibility. Since RDGs have a certain amount of active power at all simulation hours according to Figure 6(a), it is expected that the flexibility sources in the operation horizon will always have nonzero active power value.


### 5.3. Examining the Status of Technical and Economic Indices of the MG

In this section, two case studies have been investigated. These case studies refer to the power flow results of the MG (the network without DG and AL) in Case I and the proposed scheme in Case II. The results concerning technical, economic, and environmental indices for the mentioned case studies are summarized in Table5. Based on this table, by increasing network flexibility, which corresponds to a decrease in flexibility tolerance (εF), the values of Cost, EM, and EENS functions and operation indices such as energy losses (EL) and maximum voltage drop (MVD) increase, while the maximum overvoltage (MOV) is almost constant. In other words, improving the flexibility of the network is commensurate with the deterioration of the economic, environmental, reliability, and operation of the MG. Because the flexibility of the MG depends on the charging and discharging performance of ALs, it is expected that the amount of energy consumed by ALs and the consequent energy consumption of the MG increase in this situation. This is in line with the increased values of functions and operation indices. As another remark, note that the proposed scheme has been able to improve the economic conditions of the MG by about 22% ((1137.35–887.55)/113735) in the conditions of high flexibility (εF = 0.01), i.e., the Cost function is reduced by 22%. Also, the environmental status and reliability of the MG in Case II with εF = 0.01 improves compared to Case I by about 29.2% and 89.5%, respectively, according to Table 5. In terms of operation indices, EL and MVD, in this case, have decreased by about 26% and 39%, compared to Case (I). However, the MOV is roughly 0.012 p.u. in the above studies compared to the power flow studies, the amount of which is less than its permitted limit, i.e., 0.05 p.u. (1.05 − 1). Therefore, on average, about 33.5% of the operation status of the MG is improved by the proposed scheme compared to Case I.Table 5
Technical and economic results.
CaseIIIεF = 0.01εF = 0.03εF = 0.05Cost ($)1137.35887.55883.26879.84EM (kg)68528485124830848153EENS (MWh)27.82.932.812.72EL (MWh)5.213.863.753.67MVD (p.u)0.0920.0560.0540.053MOV (p.u)00.0120.0120.012


## 5.1. Examining the Feasibility of the Proposed Solution Process

Table3 provides the results of the Pareto front for the proposed problem described by equations (1)–(19). In this table, four values of zero, 0.33, 0.5, and 1 are considered for each weighting coefficient ωC, ωE, and ωR. In the second to fourth rows of this table, the values of the objective functions Cost, EM, and EENS for the case studies are ωC = 1, ωE = 1, and ωR = 1, respectively. For ωC = 1, the minimum value of Cost ($ 873.26) is calculated because in this case, only the Cost function appears in the objective function (1). These conditions are available for EM and EENS for ωE = 1 and ωR = 1, respectively, where their minimum values are 47793 kg and 2.33 MWh, respectively. Also, the maximum values of the mentioned functions can be seen in these lines, so that the maximum values of Cost, EM, and EENS are $ 943.31, 55633 kg, and 8.67 MWh, respectively. Note that according to this table, it is observed that the direction of increase and decrease of one function is not in the same direction of other functions. For instance, the decrease of EENS is proportional to the increase in EM and Cost. The same is true for other functions. Eventually, to reach a compromise between these functions, this paper employs the fuzzy decision-making method. To this end, Table 4 reports the results of the best compromise solution obtained from various Pareto optimization techniques. In other words, to evaluate the Pareto technique, the results of other methods such as ε-constraint [47], the sum of normalized functions [41], and nondominated sorting genetic algorithm (NSGA) [59] are also given in this table. According to Table 4, the weighted functions sum method, described in Section 3.1, succeeds to obtain the minimum value for the three objective functions, namely, Cost, EM, and EENS. The ε-constraint method also provides the best compromise solution close to the results of the method presented in this paper. However, the best compromise point obtained from the other two methods has a significant difference compared to the method of weighted functions sum. The compromise point obtained in the proposed method is close to the minimum value of each objective function so that the distance between Cost at this point compared to its minimum value is approximately 4.87% ((879.84–873.26)/873.26). This value for EM and EENS power is roughly 0.75% and 16.7%, respectively.Table 3
The Prato front results for the proposed scheme.
ωCωEωRCost ($)EM (kg)EENS (MWh)100873.26556337.91010943.31477938.67001937.84534292.330.50.50890.23503628.060.500.5882.72540123.7800.50.5928.45489753.490.330.330.33895.68505533.61Table 4
The best compromise solution of the proposed scheme.
MethodCost ($)EM (kg)EENS (MWh)Sum of weighted function87984481532.72ε-constraint [47]88072482082.79Normalized objective function [41]88824490333.08Nondominated sorting genetic algorithm [59]88455487212.93Figure4 depicts the convergence results of the problem using the solvers of the DE [25], Teaching-Learning-based Optimization (TLBO) [60], CSA, ALO, and the suggested hybrid ALO + CSA. In these algorithms, the population size and maximum convergence iterations are 50 and 1000, respectively, and other adjustment parameters related to each solver are selected from Refs. [24, 25, 48, 60]. In addition, the problem is solved 20 times by each algorithm to calculate statistical indices such as standard deviation (SD) of the response. Finally, based on this figure, it can be seen that the proposed ALO + CSA algorithm obtains the minimum value for the objective function F, equation (1), compared to the other mentioned algorithms. Thus, according to Figure 4, it finds the best compromise point in the fewest number of convergence iterations (CI), which is 608. This convergence iteration corresponds to the minimum possible computational time (CT), which is 212.3 s. Other solvers require a CT higher than 225 s to achieve the best compromise solution in the proposed scheme, equations (1)–(19). In addition, the hybrid ALO + CSA algorithm has an SD of 0.97%, but this value for other algorithms is greater than 1.35%. This means that the proposed solver has a low dispersion in the final response relative to the other solvers, and has almost a unique response condition. Figure 5 shows the PDF of the objective functions Cost, EM, and EENS. In this figure, the PDF of each function is presented for RWM results, and also, the standard PDF determinates using the Gaussian method (GM). As can be seen from this figure, the minimum value of the Cost function reported in Table 4, which is $ 879.84, has a probability of about 8%. Also, the PDF of Cost is a nonnormal distribution, while the PDFs of EM and EENS have approximately normal distributions. Moreover, the minimum values of EM, i.e., 48153 kg, and EENS, i.e., 2.72 MWh, have a probability of 7% and 9%, respectively, according to Figure 5.Figure 4
Convergence curve of the proposed problem based on the different solvers.Figure 5
PDF of (a) Cost, (b) EM, and (c) EENS.
(a)(b)(c)


## 5.2. Evaluating the Performance of DGs and ALs in the MG

The expected daily curve of active and reactive power of DGs is presented in Figure6. Comparing Figures 3(c) and 6(a) and according to the data of Section 4.1, it is observed that PVs and WSs in all simulation hours inject active power equal to their maximum capacity into the MG proportional to climatic conditions. According to Table 2, their pollution coefficients and operating costs are zero; it is expected that they will inject high active power into the MG to minimize EM and Cost. Also, the minimization of the EENS is in proportion to the fact that local sources are responsible for supplying consumers [23]. Therefore, to minimize the EENS during N-1 contingency, RDGs need to inject high power into the MG. Besides, at hours 1:00 to 6:00, as the price of energy purchased from the upstream network, Figure 3(b), is lower than the fuel price of FCs and MTs, as shown in Table 2, NRDGs inject low active power into the MG during these hours to minimize environmental pollution and earthquake-induced shutdown. Nevertheless, at other times when the fuel cost of NRDGs is less than the energy price, they inject active power equal to their maximum capacity into the network. This performance of FCs and MTs from 7:00 to 24:00 is commensurate with the minimization of Cost, EM, and EENS functions. In addition, it can be seen from Figure 6(b) that DGs have allocated part of it to inject reactive power into the MG at times they have free capacity. FCs and MTs inject about 0.6 and 0.4 p.u. Reactive power into the MG during 1:00–6:00. At these hours, according to Figure 7, ALs perform charging operations and receive active power from the MG. Therefore, it is expected that there will be a high voltage drop in the MG during these hours. To compensate for this, NRDGs inject significant reactive power into the network. However, RDGs generally inject reactive power into the MG during most simulation hours. This performance is proportional to compensating for voltage drop caused by the consumption of ALs and passive loads in the MG.Figure 6
Expected daily curve of (a) DGs active power and (b) DGs reactive power.
(a)(b)Figure 7
Expected daily curve of (a) MG load and (b) battery active power.
(a)(b)Figure7 illustrates the expected daily performance curves of ALs such as DRPs and batteries. According to Figure 7(b), it can be observed that all residential, commercial, and industrial loads that have participated in the proposed DRP will increase their consumption between 1:00–7:00 and 23:00–00:00 because the energy price according to Figure 3(b) has a low value during these hours. Yet, in peak hours when the energy price according to Figure 3(b) has a high value of 30 $/MWh, the mentioned consumers operate in the low-consumption mode. This demand-side management approach is proportional to minimizing all three objective functions of operating costs, environmental pollution, and MG shutdown due to N – 1 contingency. In addition, battery performance is similar to that of consumers that participate in the proposed DRP so that they are charged during 1:00–6:00 and 23:00–00:00 to receive the energy needed for discharge during peak hours. Finally, it is noteworthy that the flexibility sources, including ALs and NRDGs, have a certain amount of active power at all simulation hours to improve network flexibility. This is because these sources must compensate for the prediction error in the power generation of RDGs to improve system flexibility. Since RDGs have a certain amount of active power at all simulation hours according to Figure 6(a), it is expected that the flexibility sources in the operation horizon will always have nonzero active power value.


## 5.3. Examining the Status of Technical and Economic Indices of the MG

In this section, two case studies have been investigated. These case studies refer to the power flow results of the MG (the network without DG and AL) in Case I and the proposed scheme in Case II. The results concerning technical, economic, and environmental indices for the mentioned case studies are summarized in Table5. Based on this table, by increasing network flexibility, which corresponds to a decrease in flexibility tolerance (εF), the values of Cost, EM, and EENS functions and operation indices such as energy losses (EL) and maximum voltage drop (MVD) increase, while the maximum overvoltage (MOV) is almost constant. In other words, improving the flexibility of the network is commensurate with the deterioration of the economic, environmental, reliability, and operation of the MG. Because the flexibility of the MG depends on the charging and discharging performance of ALs, it is expected that the amount of energy consumed by ALs and the consequent energy consumption of the MG increase in this situation. This is in line with the increased values of functions and operation indices. As another remark, note that the proposed scheme has been able to improve the economic conditions of the MG by about 22% ((1137.35–887.55)/113735) in the conditions of high flexibility (εF = 0.01), i.e., the Cost function is reduced by 22%. Also, the environmental status and reliability of the MG in Case II with εF = 0.01 improves compared to Case I by about 29.2% and 89.5%, respectively, according to Table 5. In terms of operation indices, EL and MVD, in this case, have decreased by about 26% and 39%, compared to Case (I). However, the MOV is roughly 0.012 p.u. in the above studies compared to the power flow studies, the amount of which is less than its permitted limit, i.e., 0.05 p.u. (1.05 − 1). Therefore, on average, about 33.5% of the operation status of the MG is improved by the proposed scheme compared to Case I.Table 5
Technical and economic results.
CaseIIIεF = 0.01εF = 0.03εF = 0.05Cost ($)1137.35887.55883.26879.84EM (kg)68528485124830848153EENS (MWh)27.82.932.812.72EL (MWh)5.213.863.753.67MVD (p.u)0.0920.0560.0540.053MOV (p.u)00.0120.0120.012


## 6. Conclusion

The study proposed energy management of an MG penetrated by DGs and ALs with the help of a flexible-reliable operation strategy. The problem was established with three objective functions subject to constraints including AC-PF equations, operating constraints, flexibility and MG reliability, and operation models of DGs and ALs. The first function minimizes the expected operating cost imposed by MG and NRDGs. The second one attempts to find the minimum level of predicted pollutant emission and, finally, the third function minimizes the EENS associated with anN-1 event. The suggested design takes into account the operation and flexibility model of the MG in the problem constraints. The probabilistic programming assists to provide accurate modeling of uncertainties related to the load, electricity cost, output power of RDGs, and accessibility of MG devices. The weighted sum of functions-based Pareto optimization was then applied to formulate the design. The hybrid ALO-CSA algorithm finds the optimal solution of the problem in minimum convergence iterations and computational time when compared with NHEAs, so its convergence speed is higher than that of NHEAs. The standard deviation of the response is roughly 0.97%, so a unique response status is found. Therefore, the mentioned algorithm is able to extract a more accurate solution with a higher convergence speed than NHEAs because it has a more optimal point with a very low standard deviation and a lower calculation time than NHEAs. Moreover, the Pareto optimization method based on the weighted sum of functions has been able to obtain a compromise point for the proposed scheme where the operating cost, EENS, and environmental pollution are close to their minimum value. Thus, these indices are about 4.9%, 0.75% and 16.7% away from their minimum values, respectively. Furthermore, economic, environmental, reliability, and operational indicators are enhanced up to 22%, 29.2%, 89.5%, and 33.5% in an MG with high flexibility level (εF = 0.01) by applying energy management and demand-side management using an incentive-based DRP model in comparison to MG power flow studies. Therefore, according to the obtained results and comparing it with the research background, it can be seen that there are advantages for the proposed scheme, which are as follows: (1) The proposed scheme has been able to simultaneously improve the economic, technical, and environmental situations in microgrid compared with power flow studies. Also, in the technical discussion, it simultaneously improves several indicators, i.e., operation, reliability, and flexibility of the microgrid. (2) The proposed design based on hybrid algorithm has achieved an almost unique solution with very low dispersion compared with nonhybrid algorithms. It also has shorter computing time. (3) This scheme obtains the probabilistic model of microgrid operation, in which the changes of an index or variable are expressed in terms of probability.



---

*Source: 4188240-2022-09-26.xml*