Bi-Level Operation Scheduling of Distribution Systems with Multi-Microgrids Considering Uncertainties

A bi-level operation scheduling of distribution system operator (DSO) and multi-microgrids (MMGs) considering both the wholesale market and retail market is presented in this paper. To this end, the upper-level optimization problem minimizes the total costs from DSO’s point of view, while the profits of microgrids (MGs) are maximized in the lower-level optimization problem. Besides, a scenario-based stochastic programming framework using the heuristic moment matching (HMM) method is developed to tackle the uncertain nature of the problem. In this regard, the HMM technique is employed to model the scenario matrix with a reduced number of scenarios, which is effectively suitable to achieve the correlations among uncertainties. In order to solve the proposed non-linear bi-level model, Karush–Kuhn–Tucker (KKT) optimality conditions and linearization techniques are employed to transform the bi-level problem into a single-level mixed-integer linear programming (MILP) optimization problem. The effectiveness of the proposed model is demonstrated on a real-test MMG system.


Introduction
Nowadays, due to the different economical, technical, and environmental challenges in distribution systems, multi-microgrids (MMGs) including clusters of local loads and distributed energy resources (DERs) such as photovoltaic (PV) systems, micro-turbines (MTs), fuel cells (FCs), wind turbines (WTs), and energy storage systems (ESSs) are introduced as a promising solution [1]. In this framework, microgrid (MG) owners interact not only with the neighboring MGs, but also with the distribution system operator (DSO) to satisfy their own objective functions [2], which results in a bi-level operation scheduling problem for decision making of DSO and MMGs.
The bi-level operation scheduling of MMGs has been investigated to some extent. In [3], a bi-level operation model is presented for optimal scheduling of DSO and MGs. The upper-level Moreover, to tackle the uncertainty nature of input data such as solar irradiance, electrical load demand, and wholesale market prices, a scenario-based stochastic programming framework is explored. To this end, the HMM technique is implemented to build scenario matrix with reduced number of scenarios. As a certain level of correlation (or stochastic dependence) among input uncertain variables exists, the HMM method employs the matrix transformation to achieve the correlations.
Briefly, the main contributions of this paper are highlighted as follows: • Proposing a bi-level operation scheduling framework for decision making of DSO and MMGs in an uncertain environment, • Proposing a scenario matrix based on the HMM method to achieve the stochastic moments and correlations among the historical scenarios, • Transforming the non-linear bi-level optimization problem into the single-level MISOCP optimization problem through linearization techniques and KKT optimality conditions.
The remainder of the paper is organized as follows: In Section 2, the scenario matrix is modeled. The bi-level optimization model of the problem is formulated in Section 3. Case study and simulation results are discussed in Section 4. Finally, the conclusions are drawn in Section 5.

Uncertainties Modelling
Due to the uncertainty nature of input data in the proposed distribution management system (DMS), stochastic operation scheduling is developed. First, the uncertainty matrix based on the HMM method is modeled. Then, appropriate transformations are used to find the stochastic moments and correlations among the historical scenarios.

Modeling of Scenario Matrix
In order to cover the uncertainty of input data such as electrical load demand, solar irradiance, and wholesale market prices, a scenario-based stochastic programming approach is implemented. In this paper, the HMM method [23] is used to model scenario matrix, including reduced number of scenarios, that can approximate the stochastic nature of the historical scenarios. Although a larger number of scenarios leads to a more accurate approximation of the original model, this could ascend the computational burden rapidly. A sample of scenario matrix comprising of M scenarios of solar irradiance, electrical load demand, and market prices is presented in Figure 1. The HMM method employs the cubic transformation and the matrix transformation to achieve the stochastic moments and correlations, respectively. In this regard, as mentioned in [21], the first four stochastic moments i.e., expectation, standard deviation, skewness and kurtosis are selected to retain the stochastic features of scenarios. The HMM process can be discussed as follows: where, 1, 2, 3 i = refer to the three considered uncertain parameters and 1, 2, 3, 4 k = refer to the first four stochastic moments.
Step 2: Randomly generate Step 3: Employ the matrix transformation [23], as presented in (2), to satisfy the correlation where, L is a lower-triangle matrix of the correlation matrix that is determined by Cholesky decomposition [23].
Step 4: Utilize cubic transformation [23], as given in (3)  Step 5: calculate the correlation error ( c ε ) and the moment errors ( m ε ), as presented in (5) and (6), respectively. This algorithm converges to the solution when the amount of calculated errors are less than predefined thresholds given in (7).
Step 6: Calculate the scenario matrix M Ω containing scenarios of PV generation, electrical load demand, and market prices as follows [23]:

Problem Formulation
Firstly, the framework of the proposed bi-level optimization model considering correlations among uncertainties is summarized in Figure 2. Then, problem formulation of the bi-level optimization model and related solution methodology are described in the following subsections.

Upper-Level: Distribution System Operator (DSO) Decision Making
At this level the objective function is formulated to minimize the cost of DSO, which includes two cost terms as follows: where, the first term denotes the cost of selling/purchasing power to/from the day-ahead wholesale market at time t represented by (10).
The second term in (9) calculates the cost of exchanging power between DSO and MGs as follows:

Constraints
To achieve a stable operation of the DSO, the following constraints should be met.
• Bus voltage limits constraint: • Line current limits constraint: • Exchanged power limit with the wholesale electricity market: In order to deal with the limited capacity of sub-transmission transformer, the exchanged power between the wholesale market and DSO should be guaranteed as follows: • Exchanged power limit between DSO and MGs: Based on the aforementioned limitations in the contract for the exchanged power between DSO and MG operators, constraint (15) should be met.

Lower-Level: Multi-Microgrids (MMGs) Decision Making
The objective function of the lower-level optimization problem is formulated to maximize the profit of pth MG. For the sake of brevity in the formulation, the output power of PVs, MTs, and ESSs are aggregated as The considered objective function is as follows: OF =M ax Pr +Pr +Pr +Pr (16) Objective function (16) includes four profit terms represented in (17)- (20).

Constraints
To achieve the stable operation of the MMGs, the following constraints should be considered in the scheduling model.

•
Exchanged power limit between DSO and each MG: • Exchanged power limit between MGp and MGq: • Operation limit of MTs: At any time, the output power of MTs should be within an allowable range.
• Operation limits of ESSs: • Power balance of pth MG:

Solution Methodology
The amount of exchanged power between MMGs and DSO acts as the mutual variable between these two optimization levels. Thus, the proposed model in (1)-(30) is a bi-level optimization problem, in which suitable transformation methods should be employed to facilitate implementing the resulting model in commercial optimization packages. In this regard, as the lower-level problem is linear and convex, the KKT optimality conditions of the lower-level is utilized to transform the bi-level model into the single-level optimization problem. By implementing the complementary constraints, the optimization model becomes non-linear, which is easily linearized by big-M method as follows [16]: , , , , where, M is a large positive number. Finally, the resultant single-level MILP optimization problem can be solved by off-the-shelf software, such as GUROBI [24] and CPLEX [25].

Numerical Results and Discussion
In this section, the test system and working scenarios are introduced first. Then, simulation results of the proposed bi-level operation scheduling problem are discussed. Uncertainty modelling is implemented in MATLAB software and run on a computer with 2.6 GHz and 4 GB of RAM. Furthermore, the resultant single-level MILP problem is coded in the GAMS platform and solved by the CPLEX [25] solver.

Test System
The performance of the proposed bi-level model is investigated using an actual 84-bus 11.4 kV Taiwan Power Company (TPC) distribution system [1,26] with two substations and 11 feeders, as shown in Figure 3. As can be seen, this system includes three interconnected MGs. The DERs including PVs, MTs, and ESSs are utilized in the MMG system under study. Ten PV panels of 500 kW with similar type and size are considered, where the predicted values of PV output powers can be found in [27]. Seven identical 2 MW REDOX batteries [28], as ESSs, are installed at buses 3, 18, 40, 50, 57, 67, and 79. Moreover, eight similar Capstone MTs [29] with the maximum output power of 5 MW and electrical efficiency of 29% are considered. All the technical parameters and the necessary characteristics related to DERs are presented in Table 1 [1]. The hourly total load demand and load profile of each MG during the scheduling day are shown in Figure 4 [1]. Moreover, the wholesale market and retail market prices for the whole 24-h scheduling horizon can be extracted from [1,16].

Simulation Results
In order to investigate the effect of wholesale market prices on the simulation results, three different case studies are defined. The proposed prices in [1] are defined as case study 1, while 10% decrement and increment in these amounts are considered as case study 2 and case study 3, respectively.

Operation Scheduling
Optimal operation scheduling of PVs, MTs, and ESSs for scenario 1 in MG1, MG2, and MG3 are proposed in Figure 5. As can be seen, ESSs in MGs store power during off-peak periods and discharge the stored power at peak hours (i.e., 10:00-12:00 and 20:00-22:00).
Due to the lower generation cost of MTs than the wholesale market price during peak load periods, MTs work at their maximum allowable power generation level. Furthermore, based on the power generation capacities and load demand of each MG, MG1 has excess power generation, where MG2 and MG3 have shortage power. In this regard, Figure 6 shows the amounts of exchanged power among MGs through retail electricity market to balance the shortage and excess of power generations. It is observed that MG1 sells power to MG2 and MG3, while MG3 purchase power from MG2 to meet the supply-demand balance.

Considering Uncertainties
In order to investigate the effect of uncertainties on the performance of the proposed approach in depth, three scenarios for the number of installed PV systems are defined. The number of PVs for base case, i.e., 10 PVs, are considered as Scenario 2 (SC2), while 2 and 50 installed PVs are considered as Scenario 1 (SC1) and Scenario 3 (SC3), respectively.
First, the 3-year (2015-2018) historical hourly data related to PV generation, electrical load demand, and wholesale market prices is extracted from [30]. Then, the HMM method is utilized to generate sufficient number of scenarios over the range of 10 to 100. By comparing the historical data with generated scenarios in Figure 7, the effectiveness of the HMM method is demonstrated for SC1, SC2, and SC3. It can be seen that as the amount of scenarios increase from 10 to 100, the amount of errors related to the first four stochastic moments (i.e., expectation (M1), standard deviation (M2), skewness (M3) and kurtosis (M4)) between the historical data and generated scenarios are reduced significantly. Thus, these results confirm that the generated scenario matrix represents the stochastic features of historical scenarios accurately. Although the number of installed PVs is increased in SC3, the moment of errors are controlled within an acceptable range (less than 6%) shown in Figure 7c. Figure 8 shows the performance of the proposed approach regarding the number of scenarios. It can be seen that by increasing the number of scenarios, the total cost is decreased and the total execution time is increased. Moreover, when the number of scenarios reaches to a specific amount; thereafter, the amount of cost is decreased slightly and the execution time is increased significantly. Thus, 50 scenarios are selected as a trade-off between the execution time increment and cost reduction. Table 2 investigates the effect of correlations on the amount of total costs for the three defined case studies. Although the amounts related to the cost of DSO and profits of MGs are changed slightly, less than 1%, the correlations among uncertainties should be considered. It can be seen that by considering the correlations among variables, the amounts of MG profits are decreased, where the amounts of DSO costs are increased in the three case studies. This is because by considering correlations, in order to deal with the worst scenario, the operation costs are increased. Moreover, the profits of MGs are increased in SC3 comparing to SC2 due to the more output power of PVs in SC3. It can be concluded from Table 2 that the profits of MGs are decreased in case study 2 compared to case study 1. Because in case study 2 the wholesale market prices are decreased, thus more purchased power from wholesale market and less DER generations lead to less profits of MGs. On the other hand, the profits related to MGs are increased in case study 3 compared to case study 1 due to the wholesale market price increment. For the same reason, the cost of DSO is increased in case study 2 and decreased in case study 3 compared to case study 1.

Large-Scale Test Systems
In order to evaluate the efficiency and scalability of the proposed approach, larger-scale systems with 119 buses [31] and 1300 buses [32] are employed. The IEEE-119 bus system has 5 MGs with 22.809 MW and 17.041 MVAr loads, where the IEEE-1300 bus system includes 14 MGs with 132.07 MW and 115.26 MVAr loads. As can be seen in Table 3, the proposed approach finds the optimal solution point within an acceptable time, which indicates the scalability of the algorithm. Moreover, it can be concluded that the average operation cost is increased as the test system becomes larger. In order to investigate the performance of the proposed HMM method, the errors of the first four moments and correlation matrix have been calculated and compared with the relevant results extracted from the Latin hypercube sampling (LHS) method [33] and the importance sampling (IS) method [34] over the range of 10 to 100 generated scenarios. The results are presented in Figure 9. As can be seen, the proposed HMM method performs much better than the other considered methods in terms of expectation, standard deviation, skewness, kurtosis, and correlation matrix errors. Moreover, as the number of generated scenarios increases, the errors of the first moments decrease. In contrast, it can be seen that the number of generated scenarios does not strictly affect the error of the correlation matrix, because the random data in the transforming method adjust the correlations among the scenarios.

Conclusions
This paper proposed the optimal operation scheduling problem of DSO and MMGs in the presence of uncertainties. In the addressed bi-level model, the upper-level problem minimized the total costs from the DSO's perspective including the costs related to the exchanged power with the wholesale market and exchanged power between the DSO and MGs. For the lower-level problem, the profit of each MG including the profits related to the DERs, exchanged power between MGs and DSO, exchanged power among neighboring MGs, and the sold power to the loads was maximized. To tackle the uncertainty issue of input data such as electrical load demand, solar irradiance, and wholesale market prices, a scenario-based stochastic programming based on the HMM method was developed to build a scenario matrix with a reduced number of scenarios. The proposed approach effectively modeled the stochastic moments and correlation among the representative historical data. By applying KKT optimality conditions, the non-linear bi-level optimization problem was transformed into an MILP problem. To investigate the satisfactory performance of the proposed model, the 84-bus TPC distribution system was considered as a real case study. The simulation results revealed that the generated scenario matrix represents the stochastic features of historical scenarios accurately as the number of scenarios increases, the number of errors related to the first four stochastic moments between the historical data and generated scenarios is reduced significantly. It was shown that the correlations could affect the economy of the DSO and MGs by 1%. Besides, in case studies with the higher electricity market prices, MGs prefer to dispatch their owned DERs at the maximum allowable range, and then supply the probable shortage of power by purchasing the power from DSO and neighboring MGs. Moreover, the scalability of the proposed algorithm was evaluated by larger-scale systems with 119 buses and 1300 buses. Furthermore, the HMM method performed much better than LHS method and IS method over the range of 10 to 100 generated scenarios in terms of expectation, standard deviation, skewness, kurtosis, and correlation matrix errors.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
Index of hours.

MGp-DSO t Pr
Profit of exchanging power between DSO and MGp ($).

MGp-MGq t Pr
Profit of exchanging power between MGp and MGq ($).