Two-Stage Optimal Microgrid Operation with a Risk-Based Hybrid Demand Response Program Considering Uncertainty

: Owing to the increasing utilization of renewable energy resources, distributed energy resources (DERs) become inevitably uncertain, and microgrid operators have di ﬃ culty in operating the power systems because of this uncertainty. In this study, we propose a two-stage optimization approach with a hybrid demand response program (DRP) considering a risk index for microgrids (MGs) under uncertainty. The risk-based hybrid DRP is presented to reduce both operational costs and uncertainty e ﬀ ect using demand response elasticity. The problem is formulated as a two-stage optimization that considers not only the expected operation costs but also risk expense of uncertainty. To address the optimization problem, an improved multi-layer artiﬁcial bee colony (IML-ABC) is incorporated into the MG operation. The e ﬀ ectiveness of the proposed approach is demonstrated through a numerical analysis based on a typical low-voltage grid-connected MG. As a result, the proposed approach can reduce the operation costs which are taken into account uncertainty in MG. Therefore, the two-stage optimal operation considering uncertainty has been su ﬃ ciently helpful for microgrid operators (MGOs) to make risk-based decisions.


Introduction
In a smart-grid environment, distributed energy resources (DERs), such as renewable energy, have gained more attention than traditional power generation units owing to the increasing trend to address environmental concerns [1,2]. Some of the most widely used renewable energy resources are wind turbine (WT) and photovoltaic (PV) systems because of the feasibility of existing related technologies. However, these forms of DERs are dependent on fluctuations and the unpredictive nature of wind and solar resources. Since the difficulty in managing distributed uncertainties by conventional power systems, microgrids (MGs) that can usefully control DERs have been introduced as a new concept [3]. MGs manage a cluster of loads and DERs, operating as a control to offer power locally. MGOs should be able to confirm the reliability of systems considering uncertain risk indices. In this regard, the accurate evaluation of MGs is a challenging task due to the uncertainties inherent in renewable energy. Moreover, methods are needed to mitigate the effects of uncertainty when scheduling optimal operations [4]. In this situation, using flexible energy, such as demand response (DR), can provide the required demand control for the system and can be used reliably in a relatively short time under the same conditions [5]. Thus, the DR program (DRP) for power systems is expected to advance steadily [6]. This is because power system infrastructure is focused on both stability and economics, and the DRP is a flexible and inexpensive resource for operating a system. MGOs can use the DRP to reduce the peak load, save on the power reserve, and ensure power reliability. The operators can also encourage customers to use less power during periods when demand reduction is required owing to uncertainties. Thus, customers engage in contracts with MGOs to reduce demand when requested, and MGOs offer incentive costs to customers that reflect the amount of demand reduced, enabling contracts to be maintained.
Various approaches have been used recently to address uncertain problems inherent in MG operation. To minimize operational costs under a deterministic and probabilistic environment, a stochastic approach for MG operation using energy storage system (ESS) was proposed in Reference [7]. In Reference [8], optimal planning for interconnected MGs under uncertainty in large-scale distribution systems was presented to improve reliability and economics. Further, to improve system operation and management efficiency, a stochastic resource planning strategy for MGs was introduced in Reference [9] to optimally manage resources required for both the generation and demand sides. In Reference [10], both the power generation of each unit and the exchange with upstream networks were assessed through the optimal operation of MGs. Artificial intelligence algorithms were used to forecast wind speeds and optimal set volumes for DERs, and ESS capacity was determined based on forecasted data to optimize the total operating costs by Motevasel et al. [11]. Moreover, a probabilistic methodology of uncertainties caused by wind and solar generation and load consumption for estimating spinning reserve requirement was presented in Reference [12]. The DRP was considered to control the frequency of a smart MG with renewable generation, and mixed-integer linear programming was used to solve the proposed model for wind-power generation under uncertainty [13]. A stochastic planning approach was suggested to model the probabilistic behavior of wind and DRs in an energy market [14]. The authors in Reference [15] studied the effects of demand-side management of appliances on the reliability, loss, and voltage profiles of power systems, and customer comfort was also considered. Using a risk-based stochastic optimization framework, the minimum required ESS to secure the desired voltage stability margin for distribution systems was computed by Jalali et al. [16]. The authors in Reference [17] presented a flexible risk control strategy with an ESS to help remedy the removal of line overload in post-contingency situations. To minimize operating costs of MGs using the genetic algorithm, the optimal energy and power capacity of energy storage systems were determined [18]. The operational costs of MGOs were reduced by selling remaining energy at a high load level by Samadi et al. [19]. This study addressed a daily power prediction module to provide MGs with solar power output data for DER scheduling. However, the reserve for compensating wind and PV power fluctuations was not considered within the daily DER schedule. The authors in Reference [20] solved the optimal power flow through particle swarm optimization on a system of MGs with WT. To date, several important studies on realizing optimal operation for MG systems accompanied with uncertainties have been conducted [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. However, the risk strategy of uncertainty was mostly not considered when computing the optimization problem. Although uncertainty is inevitable in MG operation, few studies have considered risk scenarios. The DRP can be used in various approaches to improve system reliability and reduce operational costs, but it has generally been used as a single solution. Therefore, it is necessary not only to operate MGs considering risk strategy but also to utilize the hybrid DRP according to the risk-averse tendency of MGOs.
In this study, we present a two-stage optimization model for optimal operation in grid-connected MG considering the DRP with a risk strategy. Load, PV, and WT are considered as uncertainty variables, and each unit is modeled through a certain probability distribution function (PDF). Monte Carlo simulation and k-means clustering are used implement the scenarios. The risk-based hybrid demand response program (RH-DRP) is proposed to reduce risk-based operational costs by determining optimal DR volume configuration in accordance with risk aversion parameter. In the two-stage optimization problem, the expected operational costs are calculated in first-stage and then, the risk-based operational costs are determined via the conditional value at risk (CVaR) index in second-stage. This is successfully addressed using the improved multi-layer artificial bee colony (IML-ABC), a modified conventional ABC algorithm that improves convergence speed. The main contributions of this paper can be summarized as follows: • The proposed two-stage optimization is developed to calculate the risk-based operation costs while considering uncertainty based on scenarios. It assists MGOs make decisions from more additional degree of freedom.

•
The effectiveness of applying the RH-DRP is evaluated in terms of the total risk-based operation costs reduction, and the superiority of the RH-DRP is demonstrated via comparison analysis with different DR strategies.

•
The IML-ABC shows better performance in searching for optimal solutions in comparison with the different optimal algorithms. Furthermore, it achieves simulation time reduction with the rapid convergence speed and can be applied for a variety of optimization problems.
The remainder of this paper is organized as follows. Section 2 introduces uncertainty modeling for WT, PV, and load. Section 3 shows the proposed RH-DRP strategy, and Section 4 provides formulations for the risk-based two-stage optimization problem. Section 5 presents the solution method used by the IML-ABC algorithm and summarizes the overall process of the optimal operation approach. Section 6 presents the numerical analysis results, and, finally, Section 7 concludes the paper.

Uncertainty Modeling
In this study, a probabilistic model is used to minimize operating costs in a microgrid considering uncertainty. Accurate prediction is not possible due to the stochastic behavior of wind and solar irradiance and it is always related to the uncertainty error of the plan for the next day. Thus, to consider more realistic compliance, we use probability density function (PDF) to model the behavior of wind, solar, and load to achieve optimal results considering uncertainties.

Wind Generation Modeling
The Weibull PDF has been regularly used to model wind speed at a forecasted time [21] and can be expressed as where d and C are the factors that characterize the Weibull PDF and determine the shape and scale, respectively, and v is the wind speed.
The output power of WT can be calculated using the WT power curve parameters: where P r , v ci , v r , and v co are the rated power, cut-in speed, rated speed, and cut-off speed of the WT, respectively.

Solar Generation Modeling
Solar generation output depends on the amount of sun irradiance. Forecasted solar power generation is commonly calculated using the beta PDF expressed as follows [22]: Energies 2020, 13, 6052 Here, si is the amount of solar irradiance (kWh/m 2 ), and a and b are the beta PDF parameters that can estimate the mean (µ s ) and standard deviation (σ s ) according to the solar irradiance data, respectively.
Depending on the characteristics of the PV panels, solar irradiance can be converted to solar power by where P pv (si) represents the amount of PV output power for irradiance, η pv is the efficiency of the PV panels, and A pv is the total surface area of the PV panels.

Load Modeling
Uncertainty in load consumption caused by the stochastic behavior of power consumers can be modeled as a typical probability function with a normal PDF [23].
where P load is load demand, and µ l and σ l are the mean value and standard deviation of the load, respectively.

Scenario Generation and Reduction
In stochastic cases, the high occurrence probability has the most effect on decision makers. However, high probability situations do not always arise, whereas low probability situations may occur at any time [24]. If the low probability scenario has a critical effect, it should be considered in the decision making. The MGO would then be confronted with a problem with uncertainty, such as renewable power and load consumption in MG operation. Our study conducted scenario modeling using Monte Carlo simulations to analyze various scenarios, including the risk situation for all low probabilities with a critical effect. The generated scenarios for the load demand (φ Ln ) and power generation of WT (φ WTn ) and PV (φ PVn ) are expressed as follows: φ WTn = (WT n , ψ WTn ), φ PVn = (PV n , ψ PVn ), where L n , WT n , and PV n are the n th scenario states of the load, wind power, and solar power, respectively; Ψ Ln , Ψ WTn , and Ψ PVn are the probability in the n th scenario of the load, wind power, and solar power, respectively; and S n is the possible states and probability of the n th scenario. The scenarios for load, wind, and PV are generated based on Equations (8)- (10). Subsequently, the integrated scenarios are modeled through Equation (11). To minimize the calculation time, in our work, k-means clustering method is utilized to divide the scenarios into groups so that the sum of the squared distances from the data objects in the group is minimal [25].
where S k-means denotes the groups divided using k-means clustering, and µ i is the mean of the points in S i .

Demand Response Strategy
Load management by adjusting consumer behavior has been implemented as DRP, which is aimed at peak shaving over high demand periods in the MGs operation [26]. The DRP includes modifying the electricity consumption patterns and incentives to promote change, and these incentives are primarily used when market prices are high or when system reliability is decreased. Hence, a new hybrid type of DRP, which considers the economics and risks of uncertainty, is proposed in this study. The RH-DRP can improve economics through load shifting in a day-ahead market and then, risk effects be reduced in response to power shortages because of uncertainties in MGs.

Demand Response Elasticity
Elasticity is an economic measure that assesses the percentage of change in demand due to price fluctuations [27]. In terms of electricity consumption, this percentage value changes as power demand varies with the changes in electricity market prices. This expected value is negative as higher electricity prices possibly results in load reductions. The elasticity of demand for electricity is calculated as follows: where E(t) is DR elasticity, and D(t) and MP(t) are the power demand capacity and electricity market price in time t, respectively. The two types of elasticity of demand are self-price elasticity and substitution elasticity. Self-price elasticity disregards period variation, but it considers the variation of consumption according to electricity price changes. Meanwhile, substitution elasticity is related to shifts in the power consumption of electricity within a day. Thus, it has a constraint that power capacity after demand response is equal to the initial power consumption.

Pattern of Each Power Consumer
The value of elasticity depends on the power consumption patterns of the consumer. Therefore, the electricity value can be categorized as follows: industrial, commercial, and residential. Industrial consumers generally have the largest power demand among the three categories. Industrial loads are closely related to running factories; thus, the possible capacities of industrial consumers to participate in the DR should be considered prior to contract signing with MGOs. Then, they would receive incentives that match their capacity to participate in the DRP, and this incentive should be worth more than the industrial consumers can gain from running a factory over a period. If the consumer cannot commit to the responsibility specified in the contract, they could be charged a penalty fee. The changed power demand capacity after the industrial DRP can be calculated using the following formula: where D new,i (t) and D old,i (t) are the power demand capacity after and before implementation of industrial DR, respectively; E i , MP 0 , π i (t), and pen(t) are the elasticity default value, reference electricity Energies 2020, 13, 6052 6 of 25 market price, incentive price per kilowatt-hour, and penalty fee for the industrial consumer at time t, respectively. Commercial and residential loads are smaller and more distributed than industrial loads. Their capacities are difficult to precontract prior to participation in a DRP; thus, the elasticity value is generally smaller than the elasticity value of industrial consumers. However, adjusting commercial and residential loads is easy as they do not suffer from economic damages (e.g., industrial consumers not being able to run their factories by participating in the DRP). Therefore, commercial and residential consumers are not obligated to pay a penalty fee and can appropriately participate in real-time transactions. The changed power demand capacity after commercial or residential DRP can be calculated using the formula as follows: where D new,cr (t) and D old,cr (t) are the power demand capacity after and before implementation of commercial and residential DR, respectively; π i (t) is incentive price per kilowatt-hour for commercial and residential consumers.

Risk-Based Hybrid Demand Response Program
In this study, we propose an RH-DRP, which not only reduces MG operational costs through economic DR in a day-ahead market but also decreases risk costs due to uncertainty by applying a risk-based DR in real-time operations. The two necessary conditions of the RH-DRP are the following: • The RH-DRP consists of economic and risk-based DR. Economic DR affects the economic improvement of MG operational scheduling and is applied via a contract with MGOs by considering the substitution elasticity in a day-ahead market. Meanwhile, risk-based DR reduces the risk costs of uncertainty and can be considered the self-price elasticity in real time without a precontract.

•
It determines the load capacity that each consumer can offer to participate in the DR. Thus, MGOs can reasonably request DR from consumers considering by the stability and economics of the power system, and the consumers should respond immediately.
As the economic DR considers substitution elasticity, it has constraint. The amount after load shifting is the same as the initial load. In a real-time market, risk-based DR participants increase system stability and reduce the risk costs of MG operations. If MGOs fail to control power shortages caused by uncertainty, a critical problem occurs, such as blackouts. Therefore, the reliable operation of MGs should be ensured when using flexible resources, such as DR. In particular, MGOs should provide additional compensation to consumers participating in risk-based DR. The incentive prices depend on the consumer type and DR participation. The formula for each incentive price is expressed as follows: where Inc e,c (t) and Inc r,c (t) are the incentive prices for the type of consumer c to participate in economic DR and risk-based DR at time t, respectively; π e,c (t) is the incentive price per kilowatt-hour for each consumer participating in economic DR at time t; π r,c (t), ac(t), and D r,c (t) are risk-based DR incentive prices per kilowatt-hour for each consumer, the additional compensation for participating in risk-based DR, and the capacity to participate in risk-based DR, which is required owing to the uncertainties for each consumer, respectively. Figure 1 illustrates the proposed two-stage optimal operation process. In the first-stage, operational costs are calculated by considering the expected scenario in a day-ahead market. This stage does not consider a probabilistic scheme and is associated with the first optimal scheduling using economic DR by the MGO. The second-stage computes power shortage volume caused by uncertainty, and then, risk-based DR is utilized to minimize risk index determined through CVaR. operational costs are calculated by considering the expected scenario in a day-ahead market. This stage does not consider a probabilistic scheme and is associated with the first optimal scheduling using economic DR by the MGO. The second-stage computes power shortage volume caused by uncertainty, and then, risk-based DR is utilized to minimize risk index determined through CVaR.

Objective Function
In the two-stage optimization model, the objective function for minimizing risk-based operational costs is as expressed follows: where ffirst(X) and fsecond(X) are the first-and second-stage objective functions, respectively, and X is the decision variable vector.

First-Stage Objective Function
The first-stage of the objective function minimizes the operation costs of the expected scenario using the following formula: where Cost(t) is the operational costs for MGs in the day-ahead market at time t.

Objective Function
In the two-stage optimization model, the objective function for minimizing risk-based operational costs is as expressed follows: where f first (X) and f second (X) are the first-and second-stage objective functions, respectively, and X is the decision variable vector.

First-Stage Objective Function
The first-stage of the objective function minimizes the operation costs of the expected scenario using the following formula: where Cost(t) is the operational costs for MGs in the day-ahead market at time t.
Operation costs are the sum of fuel costs for each generator, including start-up/shut-down costs, interactions between the utility and MG, and incentive prices of participating in economic DR. The detailed expression is as follows: Energies 2020, 13, 6052 8 of 25 Inc e,c (t), (20) where P Gi (t), P sj (t), P PVk (t), and P Wl (t) are the active power outputs of the ith diesel-generator (DG), jth storage device, kth PV panel, and lth wind generator at time t, respectively; B Gi (t), B sj (t), B PVk (t), and B Wl (t) are the bids of the generator, energy storage, solar energy, and wind energy at time t, respectively; S Gi and S sj are the start-up or shut-down costs for the ith generator and jth storage, respectively; P Grid (t) and B Grid (t) are the active power and bid price, which are bought and sold with the utility at time t, respectively; N G , N s , N pv , N w , and N c are the total number of generators, storage devices, PV, wind generator units, and economic DR consumers, respectively.
When calculating operational costs, X is considered for each unit of output power, amount of load reduction of DR, and on/off mode in a day-ahead market, which can be calculated as follows: where P g and U g are the active power and state vector of all units during time t, respectively, and n and T are the numbers of decision variables and periods, respectively.

Second-Stage Objective Function
In stochastic optimal operation for MGs, operational costs vary according to each scenario, and some scenarios are expected to be very low or even have a negative effect. To consider such risk-based scenarios in MG operation, a risk management criterion is included in the mathematical formulations, where risk management implies reducing the negative effect of uncertainty. Several indices for risk management are standard deviation, shortfall probability, value-at-risk (VaR), CVaR, and so on [28]. To evaluate the risk costs of power shortages caused by uncertainties in wind, PV, and load, we used the α confidence level CVaR (α-CVaR). The α-CVaR is determined using the expected costs of (1 − α) × 100% for all worst-case scenarios. The mathematical formula is expressed as follows: where CVaR α and C s are the risk costs calculated using the α-CVaR and the stochastic operational costs in scenario s, respectively; η is the smallest cost for a given confidence level α. VaR α is the α confidence level VaR [29], as shown in Equation (26). In risk-based MG operation, stochastic operational costs are considered power shortages due to uncertainties; they include the costs of purchasing real-time market and risk-based DR incentives, which are expressed by where P rGrid (t) and B rGrid (t) are the power and bid price purchased in the real-time market at time t, respectively. Figure 2 depicts the procedure to calculate α-CVaR according to power shortages. The power shortage for each scenario can be computed based on the difference from the expected scenario in the day-ahead model; then, the α% worst scenarios are determined to calculate the α-CVaR.
In risk-based MG operation, stochastic operational costs are considered power shortages due to uncertainties; they include the costs of purchasing real-time market and risk-based DR incentives, which are expressed by where PrGrid(t) and BrGrid(t) are the power and bid price purchased in the real-time market at time t, respectively. Figure 2 depicts the procedure to calculate α-CVaR according to power shortages. The power shortage for each scenario can be computed based on the difference from the expected scenario in the day-ahead model; then, the α% worst scenarios are determined to calculate the α-CVaR.
The second stage of the objective function is to minimize the α-CVaR in real time. The MGO can set parameter according to risk aversion tendency, and then the amount of power generation and risk-based DR capacity can be adjusted through this stage.
where β is the risk aversion parameter, which can generally be adjusted from 0 to 1.

Power Balance Constraints
The total power generation by each unit should satisfy and be equal to the power demand. ,exp ,exp ,exp The second stage of the objective function is to minimize the α-CVaR in real time. The MGO can set parameter according to risk aversion tendency, and then the amount of power generation and risk-based DR capacity can be adjusted through this stage.
where β is the risk aversion parameter, which can generally be adjusted from 0 to 1.

Power Balance Constraints
The total power generation by each unit should satisfy and be equal to the power demand.
Energies 2020, 13, 6052 10 of 25 Here, P DR , P Lm , and N m are the total power capacity of the participating RH-DRP, amount of power in the m th demand level, and total number of demands, respectively. Equations (31)-(33) represent the actual values of PV, WT, and load, which consist of the expected value and error.

Power Generation Constraints
The power output of each generator unit is limited by the lower and upper bounds.
Here, P G,min (t), P PV,min (t), P Wl,min (t), and P Grid,min (t) are the lower bounds of the active power generated by DG, PV, WT, and the transaction with utility, respectively; P G,max (t), P PV,max (t), P Wl,max (t), and P Grid,max (t) are the upper bounds of the active power of each unit at time t, respectively.

RH-DRP Constraints
The RH-DRP determines the maximum participating capacity under a contract with MGO, and the sum of each DR capacity is less than the maximum DR capacity.
where D max,c (t) is the maximum DR capacity of c at time t.

Energy Storage Constraints
Some limits on the charge and discharge capacities of ESS during each time interval exist, which can be calculated as follows: W ess,min ≤ W ess,t ≤ W ess,max , where W ess,t is the amount of energy in the battery at time t, η charge (η discharge ) is the efficiency of the ESS (dis)charging, P charge (P discharge ) is the permitted capacity of the ESS (dis)charging during a ∆t; W ess,min and W ess,max are the minimum and maximum limits on the amount of energy storage in ESS, respectively; P charge,max (P discharge,max ) is the maximum capacity of ESS (dis)charging during interval ∆t.

Improved Multi-Level Artificial Bee Colony Algorithm
There are many optimization methods when solving an optimization problem. Among them, the ABC algorithm is widely used in various fields recently [30]. Especially, the ABC algorithm has been widely used to solve the MG operation problem and it has been modified to improve performance in recent years [31][32][33]. Generally, the ABC algorithm is based on the food foraging behavior of bees. A colony of bees is divided into three groups: employed bees, onlooker bees, and scout bees. Each type performs its respective tasks to find the most abundant resources. In population-based optimization techniques, such as the ABC algorithm, finding the first global minimum value is critical. Depending on the closeness of the first global minimum value to the final optimization value, the convergence speed and optimization accuracy be improved. Accordingly, the population and scope of exploration are crucial to the optimization algorithm. Here, the IML-ABC algorithm is proposed by adjusting the number of bees and food sources based on the number of iterations for finding faster and more accurate optimal solution. The IML-ABC algorithm consists of five levels (i.e., adjusting, initialization, employed bees, onlooker bees, and scout bees): (i) Adjusting level: In the IML-ABC algorithm, the number of bees increases in the initial iteration and decreases according to the number of iterations to find the initial global minimum efficiently.
Here, FN n and FN are the number of food sources in the nth iteration and mean number of food sources, respectively. In addition, σ FN , λ, and iter are the rate of food source adjustment, value of dividing the phase of the iteration, and number of iterations, respectively. (ii) Initialization level: The IML-ABC algorithm makes a random initial population of food source positions. Each food source x i (i = 1, 2, . . . , SN) has a D-dimensional problem space and can be expressed as where x ij is the j th decision variable of the i th solution vector, and x j min and x j max represent the lower and upper limit values of the j components for the X i vector, respectively. (iii) Employed bees level: Each bee constantly explores to find a food source. When it finds the optimum solution, it selects a new and best position (v ij ) close to the reference position. In the IML-ABC algorithm, to determine a better global minimum solution, employed bees explore large ranges in initial iterations and search for the best solution in a narrow range as the number of iterations increases. This can be computed as follows: where v ij represents the new position of the food source i for the jth component, ϕ ij is a random number in the range [−1, 1], and R is the value of the exploration range reduction based on the number of iterations. (iv) Onlooker bees' level: An onlooker bee finds new positions that are closed to the old position and searches for food sources according to the probability value associated with the corresponding food source. Then, the greedy method is applied; the probability value of the selected food source can be expressed as follows: where P i and fit i are the probability value and fitness value of the ith food source evaluated by the ith employed bee, respectively; and F i is the value of the objective function for the X i solution. (v) Scout bee level: If solution of X i cannot be improved through the number of predetermined trials (limit), this solution is abandoned, and the corresponding employed bee is converted to a scout bee. This scout bee randomly attempts to find a new food source in the solution space to replace the abandoned solution using Equation (44). This procedure is repeated for several maximum cycles. Figure 3 summarizes the two-stage solution process for the optimal scheduling of MGs with uncertainties using RH-DRP. The procedure is performed sequentially as follows:

Improved Multi-Level Artificial Bee Colony Algorithm
Energies 2020, 13, x FOR PEER REVIEW 13 of 27 Figure 3. Overall process of the proposed optimal operation approach.

Input Data
The proposed two-stage risk-based optimal operation was examined on a low-voltage grid-connected MG (Figure 4). This MG system contained micro turbine (MT), fuel cell (FC), PV, WT, and ESS technologies. Table 1 depicts the unit data, including the bid, start-up/shut-down costs, and minimum/maximum power [21]. The total expected daily load was 1695 kWh, divided among industrial, commercial, and residential consumers, accounting for 50%, 33.33%, and 16.67% of the Calculate the MG input data based on fixed system information.
Step 1. Establish a stochastic model of uncertainties arising from renewable energies and load.
Step 2. Determine the RH-DRP participation power capacity for each consumer type. Step 3. Shift load considering the economic DR with demand elasticity.
Step 4. Calculate expected scenario operational costs through first-stage objective function.
Step 5. Generate the scenarios using Monte Carlo simulation and reduce the scenarios using the k-means clustering technique.
Step 7. Start the IML-ABC algorithm loop for the ith scenario.
Step 8. Find the best scheduling for each scenario with the goal to minimize operation costs.
Step 9. Array the scenarios in a probability distribution and calculate CVaR α .
Step 10. Output the risk-based optimal operation costs in the MG.

Input Data
The proposed two-stage risk-based optimal operation was examined on a low-voltage grid-connected MG (Figure 4). This MG system contained micro turbine (MT), fuel cell (FC), PV, WT, and ESS technologies. Table 1 depicts the unit data, including the bid, start-up/shut-down costs, and minimum/maximum power [21]. The total expected daily load was 1695 kWh, divided among industrial, commercial, and residential consumers, accounting for 50%, 33.33%, and 16.67% of the total initial demand, respectively [34]. Taking into account the increased load on uncertainty, we assumed that 20% of each group of consumers can participate in the DRP. Table 2 lists the elasticity value, incentive costs per power, and initial demand for each consumer [35].      Figure 5 shows the day-ahead market price for grid-connected MG operation. Here, the real-time market price is assumed to be 10% more expensive than the day-ahead market price, which is the APX hourly market price [36]. Table 3 presents the WT data and Weibull PDF factors. The solar panel was a 25 kW SOLAREX MSX (US) composed of 10 × 2.5 kW panels with 18.6% efficiency and 10 m 2 of total surface area [37]. In the ESS consisting of a nickel-metal hydride battery, charging and discharging efficiency was 95%, and the minimum and maximum amounts of storage were 5% and 100% of the battery capacity, respectively. Power generation units participate in the MG depending on their technical and economic features, and excess power is exchanged with the utility through the point of common coupling (PCC). Our operation model was performed with the proposed IML-ABC algorithm set to FN = 200 and iter max = 500. All cases were simulated in MATLAB (R2019a) on a laptop computer with a 2.9 GHz Intel Core i5-9400 CPU and 16 GB RAM.

Proposed Optimization Results
In the first-stage, the MGO uses economic DR among the RH-DRP to optimize MG operation in a day-ahead market with expected scenario. This stage was simulated by considering the hourly average values of load, wind speed, and solar irradiance. Table 4 presents the expected power production of PV and WT for each hour. Here, it is assumed that the MGO must buy all the power produced by the PV/WT at each time of the day. In the second stage, risk costs arising from uncertainty are minimized by the MGO. Figure 6 show the generation scenarios and reduction results using the PDF for each unit. The 4000 scenarios were generated by considering the uncertainty in the MG via Monte Carlo simulations. Subsequently, 20 scenarios were selected via k-means clustering considering the cluster number estimation and simulation times. From the 20 generated scenarios, risk-based operational costs were determined by considering α-CVaR using the RH-DRP.

Proposed Optimization Results
In the first-stage, the MGO uses economic DR among the RH-DRP to optimize MG operation in a day-ahead market with expected scenario. This stage was simulated by considering the hourly average values of load, wind speed, and solar irradiance. Table 4 presents the expected power production of PV and WT for each hour. Here, it is assumed that the MGO must buy all the power produced by the PV/WT at each time of the day. In the second stage, risk costs arising from uncertainty are minimized by the MGO. Figure 6 show the generation scenarios and reduction results using the PDF for each unit. The 4000 scenarios were generated by considering the uncertainty in the MG via Monte Carlo simulations. Subsequently, 20 scenarios were selected via k-means clustering considering the cluster number estimation and simulation times. From the 20 generated scenarios, risk-based operational costs were determined by considering α-CVaR using the RH-DRP.    Figure 7 illustrates the expected and risk-based scenario results of the daily load and power generation capacity of PV and WT. Here, the risk-based scenario is the average value comprising the 10% worst-case scenarios. The daily load increased by 162.69 kW, whereas the total power generated by PV and WT decreased by 9.62 kW and 5.25 kW, compared with the expected scenario, respectively.   Figure 8 shows the optimal results of expected operational costs and CVaR by changing the volume ratio of RH-DRP depending on the degree of risk aversion of MGO. We simulate when the risk aversion parameters are 0.1, 0.5, and 1 to consider various MGOs. MGO with β of 0.1 tends to ignore risk-based scenarios. MGO with β of 0.5 responds appropriately to risk-based scenarios, but do not attempt to avoid them entirely. MGO with β equal to 1 tries to completely avoid risk-based scenarios. When the volume of economic DR increases, the expected costs of MG operation decrease  Figure 8 shows the optimal results of expected operational costs and CVaR by changing the volume ratio of RH-DRP depending on the degree of risk aversion of MGO. We simulate when the risk aversion parameters are 0.1, 0.5, and 1 to consider various MGOs. MGO with β of 0.1 tends to ignore risk-based scenarios. MGO with β of 0.5 responds appropriately to risk-based scenarios, but do not attempt to avoid them entirely. MGO with β equal to 1 tries to completely avoid risk-based scenarios. When the volume of economic DR increases, the expected costs of MG operation decrease because the MGO can sell surplus power to the utility through load shifting at peak load times. However, risk-based DR capacity is not sufficient to address uncertainty problems; thus, the value of CVaR significantly increases as shown in Figure 8. Meanwhile, the high proportion of risk-based DR increases the expected costs and decreases the CVaR value. In addition, by increasing β, the MGO can decrease CVaR more effectively because they are more risk averse in MG operation.
Energies 2020, 13, x FOR PEER REVIEW 18 of 27 because the MGO can sell surplus power to the utility through load shifting at peak load times. However, risk-based DR capacity is not sufficient to address uncertainty problems; thus, the value of CVaR significantly increases as shown in Figure 8. Meanwhile, the high proportion of risk-based DR increases the expected costs and decreases the CVaR value. In addition, by increasing β, the MGO can decrease CVaR more effectively because they are more risk averse in MG operation.  Figure 9 shows risk-based operational costs depending on risk aversion parameter β (0.1, 0.5, or 1). As shown in Figure 9, when β is 0.1, 0.5, and 1, the optimal solution is that the ratio of economic and risk-based DR is consist of 17.5:2.5, 10:10, and 5:15, respectively. As observed in these  Figure 9 shows risk-based operational costs depending on risk aversion parameter β (0.1, 0.5, or 1). As shown in Figure 9, when β is 0.1, 0.5, and 1, the optimal solution is that the ratio of economic and risk-based DR is consist of 17.5:2.5, 10:10, and 5:15, respectively. As observed in these results, it should be noted that the MGO who can highly avoid risk should operate within the higher portion of the risk-based DR in RH-DRP. Table 5 summarizes the results of the risk-based optimal operation in MG. The lowest expected costs are determined when β is 0.1, but CVaR is very high, owing to insufficient management of risk problems. Meanwhile, the expected costs are slightly high when β is equal to 1, but the CVaR is significantly lower than in other cases. In other words, in a 10% worst scenario, in a low-voltage grid-connected MG operation, an MGO would incur operational costs of $622.8660 (β = 0.1), $338.3607 (β = 0.5), or $293.6971 (β = 1). As can be seen from these results, RH-DRP is properly utilized according to the degree of risk aversion of MGO, and it is particularly effective in reducing risk costs for MGO that tend to avoid risk situation caused by uncertainty.
Energies 2020, 13, x FOR PEER REVIEW  19 of 27 results, it should be noted that the MGO who can highly avoid risk should operate within the higher portion of the risk-based DR in RH-DRP. Table 5 summarizes the results of the risk-based optimal operation in MG. The lowest expected costs are determined when β is 0.1, but CVaR is very high, owing to insufficient management of risk problems. Meanwhile, the expected costs are slightly high when β is equal to 1, but the CVaR is significantly lower than in other cases. In other words, in a 10% worst scenario, in a low-voltage grid-connected MG operation, an MGO would incur operational costs of $622.8660 (β = 0.1), $338.3607 (β = 0.5), or $293.6971 (β = 1). As can be seen from these results, RH-DRP is properly utilized according to the degree of risk aversion of MGO, and it is particularly effective in reducing risk costs for MGO that tend to avoid risk situation caused by uncertainty.   Tables 6-8 present the optimal power scheduling in risk-based MG operation when β is equal to 0.1, 0.5, and 1, respectively. In early periods, the ESS charged and MT output are reduced due to cheap market prices. However, at peak load hours, the battery is fully discharged, and the MGO maximizes local generation to decrease operational costs to sell a considerable amount of energy to the main grid. Owing to low cost of power generation by the FC, the MGO then decides to use its maximum capacity for power supply. Based on the high operational cost of the MT, it is scheduled flexibly while considering the market price. Figure 10 shows RH-DRP participation in MG operation over time. The economic DR is primarily used from 9:00 to 17:00 to reduce the load in common. Through load reduction, the MGO can reduce operational costs by selling remaining power to the main grid when market prices are high. The risk-based DR responds to a power shortage caused by uncertainty and supports MG operation in terms of stability. In MG optimal operation, a larger β value confers a greater risk-based DR portion in RH-DRP. Therefore, it can be confirmed that RH-DRP is efficient for risk-based MG operations considering the risk aversion tendency.

Comparison of DR Analysis
In the comparison analysis, the effects of RH-DRP are verified by comparing with a different DR strategy for the typical low-voltage MGs operation. To evaluate these strategies, three cases are implemented as follows: All cases are simulated under the same constraints, but the DRP is different. In Case 1, risk-based costs are determined to purchase only real-time market prices without DRP. Case 2, which considers the prevalent economic DRP, is used to minimize operations costs but does not respond to the risk problem. Case 3 is applied with the proposed DRP strategy, which is the optimal risk-based operation in MG using RH-DRP. Figure 11 indicates the risk-based operational costs and CVaR at various risk parameter for each case. Case 1 shows high CVaR and risk-based operational costs as peak load reduction and risk control constraints are not considered. As Case 2 utilizes all DRP capacity as economic DR, it shows lower operational costs than Case 1; however, the greater the β value to operate the MG, the higher the risk-based operating costs owing to the lack of risk management capabilities. For Case 3, the CVaR is generally lower that of the other cases because the volume ratio of RH-DRP is appropriately adjusted according to β and risk management is effectively performed; furthermore, owing to proper risk management, the risk-based operational cost does not increase significantly, even when MG is operated with a large value of β. Therefore, when MGOs operate MGs considering uncertainty problems, Case 3 promotes optimal economic and stable operation. All cases are simulated under the same constraints, but the DRP is different. In Case 1, risk-based costs are determined to purchase only real-time market prices without DRP. Case 2, which considers the prevalent economic DRP, is used to minimize operations costs but does not respond to the risk problem. Case 3 is applied with the proposed DRP strategy, which is the optimal risk-based operation in MG using RH-DRP. Figure 11 indicates the risk-based operational costs and CVaR at various risk parameter for each case. Case 1 shows high CVaR and risk-based operational costs as peak load reduction and risk control constraints are not considered. As Case 2 utilizes all DRP capacity as economic DR, it shows lower operational costs than Case 1; however, the greater the β value to operate the MG, the higher the risk-based operating costs owing to the lack of risk management capabilities. For Case 3, the CVaR is generally lower that of the other cases because the volume ratio of RH-DRP is appropriately adjusted according to β and risk management is effectively performed; furthermore, owing to proper risk management, the risk-based operational cost does not increase significantly, even when MG is operated with a large value of β. Therefore, when MGOs operate MGs considering uncertainty problems, Case 3 promotes optimal economic and stable operation.

Performance Test for the IML-ABC
To demonstrate the superiority of IML-ABC, performance tests were evaluated without DRP to fairly compare the IML-ABC with various algorithms. Table 9 presents optimal results for the forecasted operation costs or each algorithm. The simulation process for different algorithms was repeated 20 times, with the results of the best, worst, and average solutions shown. The IML-ABC algorithm reached an optimal solution compared with the other algorithms. The IML-ABC performs better than the conventional ABC algorithm in terms of simulation time. Figure 12 illustrates the convergence performances of ABC and IML-ABC. The IML-ABC algorithm found the optimal solution during the first iteration, and, when compared with the ABC algorithm, it improved the initial convergence rate. The IML-ABC and ABC reached less than a 0.1% error value compared for the final solution at 168 and 423 iterations, respectively. Therefore, these results verify

Performance Test for the IML-ABC
To demonstrate the superiority of IML-ABC, performance tests were evaluated without DRP to fairly compare the IML-ABC with various algorithms. Table 9 presents optimal results for the forecasted operation costs or each algorithm. The simulation process for different algorithms was repeated 20 times, with the results of the best, worst, and average solutions shown. The IML-ABC algorithm reached an optimal solution compared with the other algorithms. The IML-ABC performs better than the conventional ABC algorithm in terms of simulation time. Figure 12 illustrates the convergence performances of ABC and IML-ABC. The IML-ABC algorithm found the optimal solution during the first iteration, and, when compared with the ABC algorithm, it improved the initial convergence rate. The IML-ABC and ABC reached less than a 0.1% error value compared for the final solution at 168 and 423 iterations, respectively. Therefore, these results verify that the simulation time necessary to find the optimal value is improved by reducing the number of unnecessary iterations. Energies 2020, 13, x FOR PEER REVIEW 24 of 27 that the simulation time necessary to find the optimal value is improved by reducing the number of unnecessary iterations.

Conclusions
In this study, we proposed two-stage optimal operation of MGs considering uncertainty problems with risk-based DR strategies. To account for uncertainty, WT, PV, and load were modeled as Weibull, Beta, and normal PDF, respectively. The proposed RH-DRP provides a hybrid solution using both economic and risk-based DR to reduce risk-based operational costs. The objective optimization problem consisted of two-stage optimization using the IML-ABC. In the first-stage, the expected operational costs were calculated by deterministic inputs; in the second stage, the risk-based operational costs considered CVaR. To verify the effectiveness of the proposed approach, simulation was conducted on a low-voltage grid-connected MG. The results demonstrate that reasonable operational costs can be determined under the risk aversion parameters compared to conventional optimization solutions. Moreover, our RH-DRP remarkably reduces risk-based operational costs via comparison with various DR strategies. The comparison with various other algorithms also confirmed the superiority of the proposed IML-ABC in identifying the optimal solution and reducing simulation times. Therefore, the proposed approach provides MGOs not only an additional degree of freedom in decision-making under uncertainty but also a solution to reduce risk-based operational costs. Our future work is under way to focus on not only reducing risk-based operating costs but also solving reliability about dynamic problems caused by uncertainty.

Conclusions
In this study, we proposed two-stage optimal operation of MGs considering uncertainty problems with risk-based DR strategies. To account for uncertainty, WT, PV, and load were modeled as Weibull, Beta, and normal PDF, respectively. The proposed RH-DRP provides a hybrid solution using both economic and risk-based DR to reduce risk-based operational costs. The objective optimization problem consisted of two-stage optimization using the IML-ABC. In the first-stage, the expected operational costs were calculated by deterministic inputs; in the second stage, the risk-based operational costs considered CVaR. To verify the effectiveness of the proposed approach, simulation was conducted on a low-voltage grid-connected MG. The results demonstrate that reasonable operational costs can be determined under the risk aversion parameters compared to conventional optimization solutions. Moreover, our RH-DRP remarkably reduces risk-based operational costs via comparison with various DR strategies. The comparison with various other algorithms also confirmed the superiority of the proposed IML-ABC in identifying the optimal solution and reducing simulation times. Therefore, the proposed approach provides MGOs not only an additional degree of freedom in decision-making under uncertainty but also a solution to reduce risk-based operational costs. Our future work is under way to focus on not only reducing risk-based operating costs but also solving reliability about dynamic problems caused by uncertainty.