# Two-Layer Coordinated Energy Management Method in the Smart Distribution Network including Multi-Microgrid Based on the Hybrid Flexible and Securable Operation Strategy

Authors: Mohammad Hosein Sabzalian; Sasan Pirouzi; Mauricio Aredes; Bruno Wanderley Franca; Ana Carolina Cunha
Journal: International Transactions on Electrical Energy Systems
 (2022)

---

## Abstract

With the advent of smart grid theory, distribution networks can include different microgrids (MGs). Therefore, to achieve the desired technical and economic objectives in these networks, there is a need for bilateral coordination between their operators. In the following, by defining an energy management problem for them, it is predicted that the mentioned goals can be achieved. Therefore, this paper presents the hybrid flexible-securable operation (HFSO) of a smart distribution network (SDN) with grid-connected multi-microgrids using a two-layer coordinated energy management strategy. In the first layer, the microgrid (MG) operator is coordinated with sources, storages, and demand response operators. This layer models the HFSO method in the grid-connected MGs, which is based on minimizing the difference between the sum of operating cost of nonrenewable distributed generations and cost of energy received from the SDN, and the sum of flexibility and security benefits. It is constrained to AC optimal power flow, flexibility and voltage security constraints, operation model of sources and storages, and demand response. The second layer concerns coordination between the MG operators and the SDN operator. Its formulation is the same as that of the first layer, except that the HFSO model is used in the SDN according to MGs power daily data obtained from the first layer problem. The strategy converts the mixed-integer nonlinear programming to linear one to obtain the optimal solution with low calculation time and error. Moreover, stochastic programming models the uncertainties of load, energy price, and renewable power. Eventually, numerical results confirm the capability of the scheme to improve technical and economic indices simultaneously. As a result, by expecting the optimal operation for sources, storage, and responsive loads, it succeeded to enhance energy loss, voltage profile, and voltage security of the mentioned networks by up to 30%, 22%, and 5%, respectively, compared to power flow studies. In addition, there was enhancement in economic and flexibility status of the SDN and MGs.


---

## Body

## 1. Introduction

### 1.1. Motivation

The distributed generations (DGs) which are growingly used recently to reduce environmental concerns are generally constituted by renewable energy sources (RESs) [1–4], nonrenewable energy sources (NRESs), and energy storage systems (ESSs) and loads participating in demand response programming (DRP) [5]. However, the mismanagement of these devices in the power system and especially in the distribution network will create new concerns such as increasing power loss and voltage deviation [6] as well as decreasing system flexibility and security [7]. Hence, different researchers have proposed roadmaps for the distribution grid in this condition, to achieve the optimal situation [8]. In fact, the smart grid concept has been used in [8] for this network. There are telecommunication platforms [9–11] in this network to obtain a smart energy management strategy. In other words, the power grid can include several microgrids (MGs) based on this concept, and then all the DGs, ESSs, and DRP-based customers in MG are controlled coordinately by using a local and a central manager or controller following information of wide-area monitoring systems (WAMSs) [12]. Therefore, the energy management approach requires a high potential for MGs and the smart distribution network (SDN) in various fields of flexibility, security, operation, etc. Overall, to satisfy these purposes, this paper proposes the two-layer energy management of the SDN covering the multi-microgrid (MMG) to achieve a flexible-securable operation in these grids with coordination between the SDN operator (SDNO) and MG operator (MGO).


### 1.2. Literature Review

There are different researches available in the field of energy management (EM) of distribution networks and MGs. A novel method is used in [13] to increase the distribution grid flexibility and reliability by determining the optimal scheduling of the DGs and ESSs. Furthermore, the optimal location of RESs, ESSs, and DRP is evaluated in [14] for the distribution system using a robust planning mechanism following the EM model. Another scheme is introduced in [15] to manage the charging and discharging modes of electric vehicles (EVs) based on the scheduling of RES generation power with coordination between residential and grid EM systems. In [16], the active power of virtual energy storage in the SDN is managed to obtain a daily flat load profile. The centralized stochastic EM model is formulated in [17] to develop a detailed and sequential procedure for aggregator operation to obtain low consumer cost and aggregator risk based on the data of day-ahead and real-time markets in active distribution networks. Besides, the demand-side management (DSM) as an approach to the EM is modeled in [18] to manage the energy of local sources and active loads in the SDN and also to improve the technical indices of the network such as bus voltages. Another solution to EM is the power management method by controlling the active and reactive power of local DGs and ESSs in the SDN simultaneously. This method is incorporated in the different researches to obtain some objectives such as the optimal operation in multiple cooperative MGs [19], flat voltage profile [20], securable distribution network [21], and minimum energy cost for EVs [22]. An SDN generally includes several MGs, in which case it is named a multi-microgrid (MMG) system. Hence, the EM system must be able to determine the optimal scheduling of all DGs, ESSs, and DRP in MGs and the SDN simultaneously. In this context, probabilistic and stochastic EM methods are introduced in [23, 24] for MMGs and interconnected MMGs with DGs, EVs, and ESSs, respectively. Furthermore, there are different EM strategies in MMGs such as EM based on a multi-agent system [25], EM according to a system architecture [26], and EM based on an aggregator to tackle the issue of the point of common coupling congestion [27]. Moreover, a multilayer energy management system is used in [28] for an SDN with MMGs, where the energy of local DGs, ESSs, and DRP-based customers in MGs is first managed, and thus the optimal energy scheduling is obtained for SDN devices based on optimal daily power of MGs.The coordinated operation of several MGs using the energy management system is discussed in [29]. The aim was to enhance reliability, reduce emissions, and optimize the operation of MGs. In the formulated problem, the estimated operating cost of MGs and nonrenewable energy systems together with expected energy not-supplied (EENS), emission, and voltage changes is minimized. Constraints of the problem include AC power flow equations associated with MGs, constraints related to reliability index, and modeling of resources and loads. A solution to the problem is provided by using a combined algorithm that adopts GWO and Teaching-Learning-Based Optimization (TLBO). A similar design has been employed in [30] and applied to an unbalanced model of multi-MGs. Reference [31] presents a robust decentralized energy management system to optimally schedule several MGs with hydrogen stations and electric vehicle parking lots. MGs enjoy a collaborative energy market in which hydrogen and electricity providers participate. Moreover, power to heat (P2H), power to hydrogen (P2H2), CHP, energy storage devices, and a demand response program are included to enhance the flexibility of the system. The authors of [32] optimally manage the energy of a combined hydrogen, heat, and power MG (CHHP-MG) that has hydrogen stations. The aim is to supply the demand for electric vehicles and micro-CHPS. There is a collaboration between the CHHP-MG and electricity and hydrogen markets so that the operation cost is decreased as much as possible. Power to heat (P2H) and power to hydrogen (P2H2) are also included to meet the needs of heat and hydrogen. In an attempt to enhance the flexibility of the integrated system and reach a low-carbon MG, a storage system with various energy types is adopted together with heat and power demand response (HPDR) programs. Reference [33] presents an MG scheduling model to optimally coordinate hydrogen storage with demand response, energy storage, and market mechanisms. To minimize the operating costs of the MG, the design uses stochastic scheduling constrained by risk.


### 1.3. Research Gaps

There are three main research gaps in the literature based on Table1 regarding the energy management of MMG, which are summarized as follows:(i)
In most research works, SDN and microgrid models are considered separately. Additionally, in most studies, direct coordination of sources, storages, and responsive loads with the SDN operator has been taken into account. This strategy refers to single-layer energy management [13–15]. According to this strategy, if the SDN operation model is examined simultaneously in the presence of MGs, the volume of information received from sources, storages, and responsive loads in the SDN operator will be very huge, which complicates the operator's decision. Moreover, in such a strategy, the goals of the SDN operator in the objective function have been considered in most research because it is only in charge of network management. However, to solve this issue, multilayered energy management can be employed so that the sources, storages, and responsive loads are coordinated with the MG operator, and the MG operators are connected only to the SDN operator. Therefore, it is expected that the volume of information required for the SDN operator will decrease, compared to that of single-layer energy management, and its processing speed will increase. In this case, the objectives of different networks can also be satisfied.(ii)
In most studies, simultaneous modeling of one or two indices has been considered; e.g., in [13–20, 22–28], only the operation index has been modeled. Yet the network has different technical and economic indices, and the improvement of one index does not guarantee the enhancement of the status of another index. For example, to reduce the operating cost of the network, it is necessary to inject high energy sources into the network, which leads to overvoltage. Hence, to tackle this challenge, it is essential to simultaneously model different indices in the network operation problem.(iii)
Generally, a nonlinear and non-convex formulation is used for the SDN energy management model in most studies. However, this formulation is solved by numerical methods or evolutionary algorithms with high calculation time. Furthermore, the best solutions are locally optimal due to non-convex forms [34–37].Table 1
Taxonomy of recent research works.
Ref.Energy management asImproved index ofProblem modelIntegratedMultilayerFLSEOPNonlinearLinear[13]✓✗✗✗✓✓✗[14]✓✗✗✗✓✓✗[15]✓✗✗✗✓✓✗[16]✓✗✗✗✓✓✗[17]✓✗✗✗✓✓✗[18]✓✗✗✗✓✓✗[19]✓✗✗✗✓✓✗[20]✓✗✗✗✓✓✗[21]✓✗✗✓✓✓✗[22]✓✗✗✗✓✗✓[23]✓✗✗✗✓✓✗[24]✓✗✗✗✓✓✗[25]✓✗✗✗✓✓✗[26]✓✗✗✗✓✓✗[27]✓✗✗✗✓✓✗[28]✗✓✗✗✓✓✗[29]✓✗✗✗✓✓✗[30]✓✗✗✗✓✓✗[31]✓✗✗✗✓✓✗[32]✓✗✗✗✓✓✗[33]✓✗✗✗✓✗✓PM✗✓✓✓✓✗✓FL: flexibility; SE: security energy; OP: operation; PM: proposed model.


### 1.4. Contributions

To overcome the first and second issues, this paper introduces the hybrid flexible-securable operation (HFSO) based on a two-layer coordinated energy management strategy (CEMS) for an SDN with grid-connected MMGs as illustrated in Figure1. In the proposed strategy, the operator of local DGs, ESSs, and DRP in the grid-connected MG is coordinated with the MGO in the first layer of CEMS. The MGO and SDN devices are also operated coordinately with SDNO in the second layer of the CEMS. Hence, the first layer of the CEMS will consist of the HFSO mechanism for MGs so that the difference between the sum of NRES operation and MG demanded energy costs and the sum of system flexibility and security benefits is minimized. This mechanism is subject to the AC power flow equations and system operation limits related to MG; network flexibility and voltage security constraints; and formulation of local DGs, ESSs, and DRP. Furthermore, the second layer of CEMS refers to the HFSO formulation in the SDN considering MG power daily data obtained from the first layer model. It should be noted that the HFSO model in the second layer possesses the same formulation as the first layer of CEMS. The problem model in two layers is mixed-integer nonlinear programming (MINLP). Therefore, to cope with the third issue, this paper applies the linear format to achieve an optimal solution satisfying low calculation time and error in comparison with the originally proposed problem. The new form follows the Taylor series approach to linearize AC power flow equations and uses a polygon to linearize the circular inequalities. In the next step, scenario-based stochastic programming (SBSP) is incorporated to model the uncertainties of load, energy price, and maximum active power of RESs via coupling the Roulette Wheel Mechanism (RWM) and Kantorovich method. The latter is realized through a few scenario samples with high probability. Note that, to obtain this condition, it is considered that telecommunication platforms [38–40] are installed in distribution networks and MGs.Figure 1
The proposed two-layer CEMS-based HFSO.Overall, the main contributions of this paper can be summarized as follows:(i)
Using a linear format for the two-layer coordinated energy management strategy to model the hybrid flexible-securable operation in the SDN with grid-connected MMGs, to reduce the amount of information in SDNO, and to improve the simplicity of system management.(ii)
Removing the coordination between the operator of local DGs, ESSs, and DRP in the grid-connected MGs with an SDNO by establishing coordination between MGOs and the SDNO.(iii)
Modeling flexibility and security benefits in the objective function and formulating the energy related to these indices in the proposed problem to investigate system flexibility, security, and operation simultaneously.


### 1.5. Paper Organization

The rest of the paper is organized as follows: Section2 presents the nonlinear formulation of the HSFO in the SDN including MMG based on two-layer CEMS. Section 3 describes the stochastic linear model of the proposed problem. Sections 4 and 5 address the numerical simulations and the main conclusions of the paper, respectively.


## 1.1. Motivation

The distributed generations (DGs) which are growingly used recently to reduce environmental concerns are generally constituted by renewable energy sources (RESs) [1–4], nonrenewable energy sources (NRESs), and energy storage systems (ESSs) and loads participating in demand response programming (DRP) [5]. However, the mismanagement of these devices in the power system and especially in the distribution network will create new concerns such as increasing power loss and voltage deviation [6] as well as decreasing system flexibility and security [7]. Hence, different researchers have proposed roadmaps for the distribution grid in this condition, to achieve the optimal situation [8]. In fact, the smart grid concept has been used in [8] for this network. There are telecommunication platforms [9–11] in this network to obtain a smart energy management strategy. In other words, the power grid can include several microgrids (MGs) based on this concept, and then all the DGs, ESSs, and DRP-based customers in MG are controlled coordinately by using a local and a central manager or controller following information of wide-area monitoring systems (WAMSs) [12]. Therefore, the energy management approach requires a high potential for MGs and the smart distribution network (SDN) in various fields of flexibility, security, operation, etc. Overall, to satisfy these purposes, this paper proposes the two-layer energy management of the SDN covering the multi-microgrid (MMG) to achieve a flexible-securable operation in these grids with coordination between the SDN operator (SDNO) and MG operator (MGO).


## 1.2. Literature Review

There are different researches available in the field of energy management (EM) of distribution networks and MGs. A novel method is used in [13] to increase the distribution grid flexibility and reliability by determining the optimal scheduling of the DGs and ESSs. Furthermore, the optimal location of RESs, ESSs, and DRP is evaluated in [14] for the distribution system using a robust planning mechanism following the EM model. Another scheme is introduced in [15] to manage the charging and discharging modes of electric vehicles (EVs) based on the scheduling of RES generation power with coordination between residential and grid EM systems. In [16], the active power of virtual energy storage in the SDN is managed to obtain a daily flat load profile. The centralized stochastic EM model is formulated in [17] to develop a detailed and sequential procedure for aggregator operation to obtain low consumer cost and aggregator risk based on the data of day-ahead and real-time markets in active distribution networks. Besides, the demand-side management (DSM) as an approach to the EM is modeled in [18] to manage the energy of local sources and active loads in the SDN and also to improve the technical indices of the network such as bus voltages. Another solution to EM is the power management method by controlling the active and reactive power of local DGs and ESSs in the SDN simultaneously. This method is incorporated in the different researches to obtain some objectives such as the optimal operation in multiple cooperative MGs [19], flat voltage profile [20], securable distribution network [21], and minimum energy cost for EVs [22]. An SDN generally includes several MGs, in which case it is named a multi-microgrid (MMG) system. Hence, the EM system must be able to determine the optimal scheduling of all DGs, ESSs, and DRP in MGs and the SDN simultaneously. In this context, probabilistic and stochastic EM methods are introduced in [23, 24] for MMGs and interconnected MMGs with DGs, EVs, and ESSs, respectively. Furthermore, there are different EM strategies in MMGs such as EM based on a multi-agent system [25], EM according to a system architecture [26], and EM based on an aggregator to tackle the issue of the point of common coupling congestion [27]. Moreover, a multilayer energy management system is used in [28] for an SDN with MMGs, where the energy of local DGs, ESSs, and DRP-based customers in MGs is first managed, and thus the optimal energy scheduling is obtained for SDN devices based on optimal daily power of MGs.The coordinated operation of several MGs using the energy management system is discussed in [29]. The aim was to enhance reliability, reduce emissions, and optimize the operation of MGs. In the formulated problem, the estimated operating cost of MGs and nonrenewable energy systems together with expected energy not-supplied (EENS), emission, and voltage changes is minimized. Constraints of the problem include AC power flow equations associated with MGs, constraints related to reliability index, and modeling of resources and loads. A solution to the problem is provided by using a combined algorithm that adopts GWO and Teaching-Learning-Based Optimization (TLBO). A similar design has been employed in [30] and applied to an unbalanced model of multi-MGs. Reference [31] presents a robust decentralized energy management system to optimally schedule several MGs with hydrogen stations and electric vehicle parking lots. MGs enjoy a collaborative energy market in which hydrogen and electricity providers participate. Moreover, power to heat (P2H), power to hydrogen (P2H2), CHP, energy storage devices, and a demand response program are included to enhance the flexibility of the system. The authors of [32] optimally manage the energy of a combined hydrogen, heat, and power MG (CHHP-MG) that has hydrogen stations. The aim is to supply the demand for electric vehicles and micro-CHPS. There is a collaboration between the CHHP-MG and electricity and hydrogen markets so that the operation cost is decreased as much as possible. Power to heat (P2H) and power to hydrogen (P2H2) are also included to meet the needs of heat and hydrogen. In an attempt to enhance the flexibility of the integrated system and reach a low-carbon MG, a storage system with various energy types is adopted together with heat and power demand response (HPDR) programs. Reference [33] presents an MG scheduling model to optimally coordinate hydrogen storage with demand response, energy storage, and market mechanisms. To minimize the operating costs of the MG, the design uses stochastic scheduling constrained by risk.


## 1.3. Research Gaps

There are three main research gaps in the literature based on Table1 regarding the energy management of MMG, which are summarized as follows:(i)
In most research works, SDN and microgrid models are considered separately. Additionally, in most studies, direct coordination of sources, storages, and responsive loads with the SDN operator has been taken into account. This strategy refers to single-layer energy management [13–15]. According to this strategy, if the SDN operation model is examined simultaneously in the presence of MGs, the volume of information received from sources, storages, and responsive loads in the SDN operator will be very huge, which complicates the operator's decision. Moreover, in such a strategy, the goals of the SDN operator in the objective function have been considered in most research because it is only in charge of network management. However, to solve this issue, multilayered energy management can be employed so that the sources, storages, and responsive loads are coordinated with the MG operator, and the MG operators are connected only to the SDN operator. Therefore, it is expected that the volume of information required for the SDN operator will decrease, compared to that of single-layer energy management, and its processing speed will increase. In this case, the objectives of different networks can also be satisfied.(ii)
In most studies, simultaneous modeling of one or two indices has been considered; e.g., in [13–20, 22–28], only the operation index has been modeled. Yet the network has different technical and economic indices, and the improvement of one index does not guarantee the enhancement of the status of another index. For example, to reduce the operating cost of the network, it is necessary to inject high energy sources into the network, which leads to overvoltage. Hence, to tackle this challenge, it is essential to simultaneously model different indices in the network operation problem.(iii)
Generally, a nonlinear and non-convex formulation is used for the SDN energy management model in most studies. However, this formulation is solved by numerical methods or evolutionary algorithms with high calculation time. Furthermore, the best solutions are locally optimal due to non-convex forms [34–37].Table 1
Taxonomy of recent research works.
Ref.Energy management asImproved index ofProblem modelIntegratedMultilayerFLSEOPNonlinearLinear[13]✓✗✗✗✓✓✗[14]✓✗✗✗✓✓✗[15]✓✗✗✗✓✓✗[16]✓✗✗✗✓✓✗[17]✓✗✗✗✓✓✗[18]✓✗✗✗✓✓✗[19]✓✗✗✗✓✓✗[20]✓✗✗✗✓✓✗[21]✓✗✗✓✓✓✗[22]✓✗✗✗✓✗✓[23]✓✗✗✗✓✓✗[24]✓✗✗✗✓✓✗[25]✓✗✗✗✓✓✗[26]✓✗✗✗✓✓✗[27]✓✗✗✗✓✓✗[28]✗✓✗✗✓✓✗[29]✓✗✗✗✓✓✗[30]✓✗✗✗✓✓✗[31]✓✗✗✗✓✓✗[32]✓✗✗✗✓✓✗[33]✓✗✗✗✓✗✓PM✗✓✓✓✓✗✓FL: flexibility; SE: security energy; OP: operation; PM: proposed model.


## 1.4. Contributions

To overcome the first and second issues, this paper introduces the hybrid flexible-securable operation (HFSO) based on a two-layer coordinated energy management strategy (CEMS) for an SDN with grid-connected MMGs as illustrated in Figure1. In the proposed strategy, the operator of local DGs, ESSs, and DRP in the grid-connected MG is coordinated with the MGO in the first layer of CEMS. The MGO and SDN devices are also operated coordinately with SDNO in the second layer of the CEMS. Hence, the first layer of the CEMS will consist of the HFSO mechanism for MGs so that the difference between the sum of NRES operation and MG demanded energy costs and the sum of system flexibility and security benefits is minimized. This mechanism is subject to the AC power flow equations and system operation limits related to MG; network flexibility and voltage security constraints; and formulation of local DGs, ESSs, and DRP. Furthermore, the second layer of CEMS refers to the HFSO formulation in the SDN considering MG power daily data obtained from the first layer model. It should be noted that the HFSO model in the second layer possesses the same formulation as the first layer of CEMS. The problem model in two layers is mixed-integer nonlinear programming (MINLP). Therefore, to cope with the third issue, this paper applies the linear format to achieve an optimal solution satisfying low calculation time and error in comparison with the originally proposed problem. The new form follows the Taylor series approach to linearize AC power flow equations and uses a polygon to linearize the circular inequalities. In the next step, scenario-based stochastic programming (SBSP) is incorporated to model the uncertainties of load, energy price, and maximum active power of RESs via coupling the Roulette Wheel Mechanism (RWM) and Kantorovich method. The latter is realized through a few scenario samples with high probability. Note that, to obtain this condition, it is considered that telecommunication platforms [38–40] are installed in distribution networks and MGs.Figure 1
The proposed two-layer CEMS-based HFSO.Overall, the main contributions of this paper can be summarized as follows:(i)
Using a linear format for the two-layer coordinated energy management strategy to model the hybrid flexible-securable operation in the SDN with grid-connected MMGs, to reduce the amount of information in SDNO, and to improve the simplicity of system management.(ii)
Removing the coordination between the operator of local DGs, ESSs, and DRP in the grid-connected MGs with an SDNO by establishing coordination between MGOs and the SDNO.(iii)
Modeling flexibility and security benefits in the objective function and formulating the energy related to these indices in the proposed problem to investigate system flexibility, security, and operation simultaneously.


## 1.5. Paper Organization

The rest of the paper is organized as follows: Section2 presents the nonlinear formulation of the HSFO in the SDN including MMG based on two-layer CEMS. Section 3 describes the stochastic linear model of the proposed problem. Sections 4 and 5 address the numerical simulations and the main conclusions of the paper, respectively.


## 2. HFSO Approach in the SDN Using Two-Layer CEMS

The HFSO model in the SDN is presented in this section according to the two-layer CEMS. The first layer problem refers to the optimal energy management of each MG based on minimizing the difference between the MG operation cost and the sum of flexibility and security benefits subject to MG power flow equations and limitations, renewable and flexible sources constraints, and security and flexibility formulations. Moreover, the proposed energy management strategy in the second layer is used for the SDN with the same purposes and constraints as in the MGs.In this article, it is assumed that the necessary telecommunication infrastructure [41–43] exists to establish a smart network. Therefore, in this situation, the problem can be implemented in distribution network operators and MGs. It is also assumed that the operation is carried out according to the normal conditions (without events) of the networks.


### 2.1. MMG Local Operation (First Layer)

In this section, the CEMS models each MG of the MMG scheme to obtain a flexible and securable local operation approach. Hence, this method is formulated as follows.


#### 2.1.1. Objective Function

The proposed problem in this section performs the minimization of the difference between the expected MG operation cost and the sum of expected flexibility and security benefits from the MG operator (MGO) viewpoint. This statement is expressed in (1) in a mathematical format, where the first part of this equation refers to the MG operation cost including the cost of energy received from the upstream network and NRES fuel cost [28]. The second part of (1) denotes the MG flexibility benefit that is equal to the product of the flexibility incentive price (FIP) and flexibility energy (FE) from all FSs, i.e., DRP, ESS, and NRES. Finally, the security benefit of MG is expressed in the third term of (1), depending on the security inactive price (SIP) and the security energy (SE) from all sources and active loads.In this paper, since an economic model for MGs and SDNs is considered, flexibility, voltage security, and network operation indices are formulated in the economic model as were cost and profit. Therefore, the objective function of the proposed problem considers the mentioned indices in an integrated manner. Note, however, that this type of modeling is based on multi-objective optimization based on a normalized objective function, in which all parts of the objective function will have the same dimension [20, 22, 34] similar to the objective function (1) in which all parts are in dollars.(1)min∑w∈Sπw∑t∈Tρt,wPms,t,wMS+∑n∈MNβnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Security benefit.


#### 2.1.2. MG Constraints

The constraints of this section include AC power flow equations, expressed in (2)–(7), and system operation limits, described in (8)–(11) [22]. It is noteworthy that (2)–(7) represent active and reactive power balance in MG buses [44], active and reactive power from MG lines, and values of voltage angle and magnitude in the slack bus, respectively. In this model, variables PMS and QMS denote, respectively, MG station active and reactive powers, existing only in the slack bus. The value of these variables is zero in other buses. In addition, the system operation limits consist of voltage limits in MG buses (8), MG line capacity limit (9), and MG station capacity and power factor limits (10) and (11). It should be noted that limitation (11) presents an equivalent equation for the power factor limit of the MG station so that the reactive power of this station can change in the defined range in (11) by considering PFmin as the minimum allowed power factor to satisfy this limitation.(2)Pn,t,wMS+Pn,t,wR+Pn,t,wNR+Pn,t,wdis−Pn,t,wch−∑j∈MNAn,jPn,j,t,wL=Ln,t,wP−Ln,t,wDR,∀n∈MN,t∈ST,w∈S,(3)Qn,t,wMS+Qn,t,wNR+Qn,t,wS−∑j∈MNAn,jQn,j,t,wL=Ln,t,wQ,∀n∈MN,t∈ST,w∈S,(4)Pn,j,t,wL=Gn,jVn,t,w−Vn,t,wVj,t,wGn,jcosδn,t,w−δj,t,w+Bn,jsinδn,t,w−δj,t,w,∀n,j∈MN,t∈ST,w∈S,(5)Qn,j,t,wL=−Bn,jVn,t,w+Vn,t,wVj,t,wBn,jcosδn,t,w−δj,t,w−Gn,jsinδn,t,w−δj,t,w,∀n,j∈MN,t∈ST,w∈S,(6)δn,t,w=δref,∀n=ms,t∈ST,w∈S,(7)Vn,t,w=Vref,∀n=ms,t∈ST,w∈S,(8)Vmin≤Vn,t,w≤Vmax,∀n∈MN,t∈ST,w∈S,(9)Pn,j,t,wL2+Qn,j,t,wL2≤Sn,jL,max2,∀n,j∈MN,t∈ST,w∈S,(10)Pn,t,wMS2+Qn,t,wMS2≤SnG,max2,∀n=ms,t∈ST,w∈S,(11)−tancos−1PFm×Pn,t,wMS≤Qn,t,wMS≤tancos−1PFm×Pn,t,wMS,∀n=ms,t∈ST,w∈S.


#### 2.1.3. DRP Constraints

The incentive model of DRP is considered in this paper, based on the fact that all loads can shift part of their consumption in peak load times to off-peak load hours, according to the electrical energy price. This statement is modeled as (12) and (13), where (12) points to the fact that the consumption energy of loads with/without DRP should be equal [45]. The limitation of the proposed DRP is also expressed in (13), where ξ defines the rate of participation of load in the DRP, which varies between 0 and 1 [45]. Finally, the DRP benefit is used in the MG operation cost model in the first part of (1), as it considers the total demanded power from the SDN, where this power depends on the operation of sources and active loads.Note that in the proposed model for DRP, consumers participating in DRP reduce their energy consumption during the hours with high energy prices and receive this amount of reduced energy during the hours with inexpensive energy from the grid. Accordingly, constraint (13) expresses the controllable range of DRP power. Constraint (12) also guarantees that all reduced energy during the hours with expensive energy prices must be provided during the hours with inexpensive energy prices. Following this DRP method, a residential consumer can be an example, such that high energy consumption loads such as laundry, ironing, exercising with exercise equipment, and other items during off-peak hours can be shifted to intervals with cheap energy prices.(12)∑t∈STLn,t,wP−Ln,t,wDR=∑t∈STLn,t,wP,∀n∈MN,w∈S,(13)−ξ×Ln,t,wP≤Ln,t,wDR≤ξ×Ln,t,wP,∀n=s,t∈ST,w∈S.


#### 2.1.4. DG Constraints

Generally, DG includes two types of energy sources that are renewable and nonrenewable, i.e., RESs and NRESs. The main purpose of RESs is the generation of active power according to the standard IEEE 1547 [46–49], and the output power of RESs depends on natural phenomena such as wind speed and solar radiation. Therefore, this statement is modeled as (14) considering PR,max varying in each simulation time and scenario sample. However, the NRES model includes its capacity limit as (15) which is constant for each time and scenario [28].(14)0≤Pn,t,wR≤Pn,t,wR,max,∀n∈MN,t∈ST,w∈S,(15)Pn,t,wNR2+Qn,t,wNR2≤SnNR,max2,∀n∈MN,t∈ST,w∈S.


#### 2.1.5. ESS Constraints

The operation model of ESS follows (16)–(21) [50–52], with constraints (16) and (17) referring, respectively, to discharge and charge rates of ESS considering binary variable x as ESS operation status; i.e., it acts in discharge mode if x = 1; otherwise, it operates in charge mode. Moreover, (18) and (19) calculate the stored energy in ESS at different simulation times, where it should vary in the defined range based on (20). Finally, the ESS charger capacity limit is modeled as (21) indicating a circular inequality with a radius of SS,max. In addition, the ESS operation cost/benefit as DRP benefit is considered in MG operation cost formulation.(16)0≤Pn,t,wdis≤DRn,t,wxn,t,∀n∈MN,t∈ST,w∈S,(17)0≤Pn,t,wch≤CRn,t,w1−xn,t,∀n∈MN,t∈ST,w∈S,(18)En,t,w=En,t−1,w+ηchPn,t,wch−1ηdisPn,t,wdis,∀n∈MN,t∈ST,w∈S,(19)En,t,w=Enini,∀n∈MN,t=0,w∈S,(20)Enmin≤En,t,w≤Enmax,∀n∈MN,t∈ST,w∈S,(21)Pn,t,wdis−Pn,t,wch2+Qn,t,wS2≤SnS,max2,∀n∈MN,t∈ST,w∈S.


#### 2.1.6. Security Constraints

Voltage security is one of the important indices in the distribution network or low voltage system [15]. This index tends to reduce in heavily loading conditions, while it increases or improves with the rise in the injection power by all sources and active loads in MG. Generally, the voltage stability margin (VSM) method based on the voltage stability index (SI-index) is used to investigate the network voltage security [21]. This method can also determine, firstly, the weakest bus of the network with the minimum SI or the worst SI (WSI); then, the value of WSI is calculated by (22), and constraint (23) should be satisfied. In addition, voltage security will be improved if the injection power rises by local sources, i.e., RES, NRES, ESS, and DRP. Therefore, security energy is equal to the sum of all injection power in the operation horizon as shown in the following equation:(22)WSIwb,t,w=.Vwb−1,t,w4−4Vwb−1,t,w2Rwb−1,wbPwb−1,wb,t,wL+Xwb−1,wbQwb−1,wb,t,wL−4Xwb−1,wbPwb−1,wb,t,wL−Rwb−1,wbQwb−1,wb,t,wL2,∀t∈ST,w∈S,(23)WSIwb,t,w≥WSImin,∀t∈ST,w∈S,(24)SwE=∑n∈MN∑t∈STLn,t,wDR+Pn,t,wR+Pn,t,wNR+Pn,t,wdis−Pn,t,wch,∀w∈S.


#### 2.1.7. Flexibility Constraints

It should be noted that flexibility in the power system is defined as “the modification of generation injection and/or consumption patterns in reaction to an external price or activation signal to provide a service within the electrical system” according to [53]. Hence, each flexible source, i.e., NRES, ESS, and DRP, includes both upward and downward flexibility services. There is upward flexibility power (F+) in scenario w for each FS if the difference between the power of FS in scenarios w and 1 (refers to the scenario that has forecasted values for uncertain parameters) is positive; otherwise, it is considered as downward flexibility power (F−) with constraints (25)–(27) [32]. Moreover, the flexibility energy is equal to the sum of upward and downward flexibility power of all FSs, as seen in (28).It is noteworthy that the flexibility quantity is equal to the ability of flexibility sources to compensate for power fluctuations of RES in different scenarios compared to its predicted value. Therefore, two upward and downward functions are obtained for flexibility sources. Upward (downward) mode occurs when the source must compensate for the RES power deficit (excess) relative to its predicted value. Since, in both modes of operation, the difference in active power for the flexibility source is measured, the quantity obtained is in watts; that is, it is a power quantity. Hence, the sum of them will represent the flexibility energy. Of course, the availability of RES sources is considered in this paper, where this is expected to happen in the future.(25)Fn,t,wDR+−Fn,t,wDR−=Ln,t,wDR−Ln,t,1DR,∀n∈MN,t∈ST,w∈S,Fn,t,wDR+,Fn,t,wDR−≥0,(26)Fn,t,wNR+−Fn,t,wNR−=Pn,t,wNR−Pn,t,1NR,∀n∈MN,t∈ST,w∈S,Fn,t,wNR+,Fn,t,wNR−≥0,(27)Fn,t,wS+−Fn,t,wS−=Pn,t,wdis−Pn,t,wch−Pn,t,1dis−Pn,t,1ch,∀n∈MN,t∈ST,w∈S,Fn,t,wS+,Fn,t,wS−≥0,(28)FwE=∑n∈MN∑t∈STFn,t,wDR++Fn,t,wDR−+Fn,t,wNR++Fn,t,wNR−+Fn,t,wS++Fn,t,wS−,∀w∈S.


### 2.2. SDN Operation (Second Layer)

The second layer of the proposed CEMS refers to the optimal scheduling management of SDN devices by the SDN operator (SDNO) according to the demand of MGs and the SDN conditions. It should be noted that the problem model in this section minimizes the difference between the expected SDN operation cost and the sum of expected flexibility and security benefits subject to the network, constraints of RESs and FSs, and security and flexibility equations. Therefore, the model of the second layer, expressed in (29)–(33), is the same as the first one. However, the difference remains in the nodal active and reactive power balance and security energy equations which include the MGs demand as (30)–(32), and the terms of ds, DN, PDS, and QDS are substituted by ms, MN, PMS, and QMS, respectively.(29)min∑w∈Sπw∑t∈Tρt,wPds,t,wDS+∑n∈SNβnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Security benefit,subject to(30)Pn,t,wDS+Pn,t,wR+Pn,t,wNR+Pn,t,wdis−Pn,t,wch−∑j∈DNAn,jPn,j,t,wL=Ln,t,wP−Ln,t,wDR+∑ms∈MLCn,msPms,t,wMS.∀n∈MN,t∈ST,w∈S,(31)Qn,t,wDS+Qn,t,wNR+Qn,t,wS−∑j∈DNAn,jQn,j,t,wL=Ln,t,wQ+∑ms∈MLCn,msQms,t,wMS,∀n∈MN,t∈ST,w∈S,(32)SwE=∑n∈DN∑t∈STLn,t,wDR+Pn,t,wR+Pn,t,wNR+Pn,t,wdis−Pn,t,wch−∑ms∈MLCn,msPms,t,wMS,∀w∈S.(33)Constraints4−23and25−28substitutingdsbyms,DNbyMN,PDSbyPMS,andQDSbyQMS.


## 2.1. MMG Local Operation (First Layer)

In this section, the CEMS models each MG of the MMG scheme to obtain a flexible and securable local operation approach. Hence, this method is formulated as follows.


### 2.1.1. Objective Function

The proposed problem in this section performs the minimization of the difference between the expected MG operation cost and the sum of expected flexibility and security benefits from the MG operator (MGO) viewpoint. This statement is expressed in (1) in a mathematical format, where the first part of this equation refers to the MG operation cost including the cost of energy received from the upstream network and NRES fuel cost [28]. The second part of (1) denotes the MG flexibility benefit that is equal to the product of the flexibility incentive price (FIP) and flexibility energy (FE) from all FSs, i.e., DRP, ESS, and NRES. Finally, the security benefit of MG is expressed in the third term of (1), depending on the security inactive price (SIP) and the security energy (SE) from all sources and active loads.In this paper, since an economic model for MGs and SDNs is considered, flexibility, voltage security, and network operation indices are formulated in the economic model as were cost and profit. Therefore, the objective function of the proposed problem considers the mentioned indices in an integrated manner. Note, however, that this type of modeling is based on multi-objective optimization based on a normalized objective function, in which all parts of the objective function will have the same dimension [20, 22, 34] similar to the objective function (1) in which all parts are in dollars.(1)min∑w∈Sπw∑t∈Tρt,wPms,t,wMS+∑n∈MNβnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Security benefit.


### 2.1.2. MG Constraints

The constraints of this section include AC power flow equations, expressed in (2)–(7), and system operation limits, described in (8)–(11) [22]. It is noteworthy that (2)–(7) represent active and reactive power balance in MG buses [44], active and reactive power from MG lines, and values of voltage angle and magnitude in the slack bus, respectively. In this model, variables PMS and QMS denote, respectively, MG station active and reactive powers, existing only in the slack bus. The value of these variables is zero in other buses. In addition, the system operation limits consist of voltage limits in MG buses (8), MG line capacity limit (9), and MG station capacity and power factor limits (10) and (11). It should be noted that limitation (11) presents an equivalent equation for the power factor limit of the MG station so that the reactive power of this station can change in the defined range in (11) by considering PFmin as the minimum allowed power factor to satisfy this limitation.(2)Pn,t,wMS+Pn,t,wR+Pn,t,wNR+Pn,t,wdis−Pn,t,wch−∑j∈MNAn,jPn,j,t,wL=Ln,t,wP−Ln,t,wDR,∀n∈MN,t∈ST,w∈S,(3)Qn,t,wMS+Qn,t,wNR+Qn,t,wS−∑j∈MNAn,jQn,j,t,wL=Ln,t,wQ,∀n∈MN,t∈ST,w∈S,(4)Pn,j,t,wL=Gn,jVn,t,w−Vn,t,wVj,t,wGn,jcosδn,t,w−δj,t,w+Bn,jsinδn,t,w−δj,t,w,∀n,j∈MN,t∈ST,w∈S,(5)Qn,j,t,wL=−Bn,jVn,t,w+Vn,t,wVj,t,wBn,jcosδn,t,w−δj,t,w−Gn,jsinδn,t,w−δj,t,w,∀n,j∈MN,t∈ST,w∈S,(6)δn,t,w=δref,∀n=ms,t∈ST,w∈S,(7)Vn,t,w=Vref,∀n=ms,t∈ST,w∈S,(8)Vmin≤Vn,t,w≤Vmax,∀n∈MN,t∈ST,w∈S,(9)Pn,j,t,wL2+Qn,j,t,wL2≤Sn,jL,max2,∀n,j∈MN,t∈ST,w∈S,(10)Pn,t,wMS2+Qn,t,wMS2≤SnG,max2,∀n=ms,t∈ST,w∈S,(11)−tancos−1PFm×Pn,t,wMS≤Qn,t,wMS≤tancos−1PFm×Pn,t,wMS,∀n=ms,t∈ST,w∈S.


### 2.1.3. DRP Constraints

The incentive model of DRP is considered in this paper, based on the fact that all loads can shift part of their consumption in peak load times to off-peak load hours, according to the electrical energy price. This statement is modeled as (12) and (13), where (12) points to the fact that the consumption energy of loads with/without DRP should be equal [45]. The limitation of the proposed DRP is also expressed in (13), where ξ defines the rate of participation of load in the DRP, which varies between 0 and 1 [45]. Finally, the DRP benefit is used in the MG operation cost model in the first part of (1), as it considers the total demanded power from the SDN, where this power depends on the operation of sources and active loads.Note that in the proposed model for DRP, consumers participating in DRP reduce their energy consumption during the hours with high energy prices and receive this amount of reduced energy during the hours with inexpensive energy from the grid. Accordingly, constraint (13) expresses the controllable range of DRP power. Constraint (12) also guarantees that all reduced energy during the hours with expensive energy prices must be provided during the hours with inexpensive energy prices. Following this DRP method, a residential consumer can be an example, such that high energy consumption loads such as laundry, ironing, exercising with exercise equipment, and other items during off-peak hours can be shifted to intervals with cheap energy prices.(12)∑t∈STLn,t,wP−Ln,t,wDR=∑t∈STLn,t,wP,∀n∈MN,w∈S,(13)−ξ×Ln,t,wP≤Ln,t,wDR≤ξ×Ln,t,wP,∀n=s,t∈ST,w∈S.


### 2.1.4. DG Constraints

Generally, DG includes two types of energy sources that are renewable and nonrenewable, i.e., RESs and NRESs. The main purpose of RESs is the generation of active power according to the standard IEEE 1547 [46–49], and the output power of RESs depends on natural phenomena such as wind speed and solar radiation. Therefore, this statement is modeled as (14) considering PR,max varying in each simulation time and scenario sample. However, the NRES model includes its capacity limit as (15) which is constant for each time and scenario [28].(14)0≤Pn,t,wR≤Pn,t,wR,max,∀n∈MN,t∈ST,w∈S,(15)Pn,t,wNR2+Qn,t,wNR2≤SnNR,max2,∀n∈MN,t∈ST,w∈S.


### 2.1.5. ESS Constraints

The operation model of ESS follows (16)–(21) [50–52], with constraints (16) and (17) referring, respectively, to discharge and charge rates of ESS considering binary variable x as ESS operation status; i.e., it acts in discharge mode if x = 1; otherwise, it operates in charge mode. Moreover, (18) and (19) calculate the stored energy in ESS at different simulation times, where it should vary in the defined range based on (20). Finally, the ESS charger capacity limit is modeled as (21) indicating a circular inequality with a radius of SS,max. In addition, the ESS operation cost/benefit as DRP benefit is considered in MG operation cost formulation.(16)0≤Pn,t,wdis≤DRn,t,wxn,t,∀n∈MN,t∈ST,w∈S,(17)0≤Pn,t,wch≤CRn,t,w1−xn,t,∀n∈MN,t∈ST,w∈S,(18)En,t,w=En,t−1,w+ηchPn,t,wch−1ηdisPn,t,wdis,∀n∈MN,t∈ST,w∈S,(19)En,t,w=Enini,∀n∈MN,t=0,w∈S,(20)Enmin≤En,t,w≤Enmax,∀n∈MN,t∈ST,w∈S,(21)Pn,t,wdis−Pn,t,wch2+Qn,t,wS2≤SnS,max2,∀n∈MN,t∈ST,w∈S.


### 2.1.6. Security Constraints

Voltage security is one of the important indices in the distribution network or low voltage system [15]. This index tends to reduce in heavily loading conditions, while it increases or improves with the rise in the injection power by all sources and active loads in MG. Generally, the voltage stability margin (VSM) method based on the voltage stability index (SI-index) is used to investigate the network voltage security [21]. This method can also determine, firstly, the weakest bus of the network with the minimum SI or the worst SI (WSI); then, the value of WSI is calculated by (22), and constraint (23) should be satisfied. In addition, voltage security will be improved if the injection power rises by local sources, i.e., RES, NRES, ESS, and DRP. Therefore, security energy is equal to the sum of all injection power in the operation horizon as shown in the following equation:(22)WSIwb,t,w=.Vwb−1,t,w4−4Vwb−1,t,w2Rwb−1,wbPwb−1,wb,t,wL+Xwb−1,wbQwb−1,wb,t,wL−4Xwb−1,wbPwb−1,wb,t,wL−Rwb−1,wbQwb−1,wb,t,wL2,∀t∈ST,w∈S,(23)WSIwb,t,w≥WSImin,∀t∈ST,w∈S,(24)SwE=∑n∈MN∑t∈STLn,t,wDR+Pn,t,wR+Pn,t,wNR+Pn,t,wdis−Pn,t,wch,∀w∈S.


### 2.1.7. Flexibility Constraints

It should be noted that flexibility in the power system is defined as “the modification of generation injection and/or consumption patterns in reaction to an external price or activation signal to provide a service within the electrical system” according to [53]. Hence, each flexible source, i.e., NRES, ESS, and DRP, includes both upward and downward flexibility services. There is upward flexibility power (F+) in scenario w for each FS if the difference between the power of FS in scenarios w and 1 (refers to the scenario that has forecasted values for uncertain parameters) is positive; otherwise, it is considered as downward flexibility power (F−) with constraints (25)–(27) [32]. Moreover, the flexibility energy is equal to the sum of upward and downward flexibility power of all FSs, as seen in (28).It is noteworthy that the flexibility quantity is equal to the ability of flexibility sources to compensate for power fluctuations of RES in different scenarios compared to its predicted value. Therefore, two upward and downward functions are obtained for flexibility sources. Upward (downward) mode occurs when the source must compensate for the RES power deficit (excess) relative to its predicted value. Since, in both modes of operation, the difference in active power for the flexibility source is measured, the quantity obtained is in watts; that is, it is a power quantity. Hence, the sum of them will represent the flexibility energy. Of course, the availability of RES sources is considered in this paper, where this is expected to happen in the future.(25)Fn,t,wDR+−Fn,t,wDR−=Ln,t,wDR−Ln,t,1DR,∀n∈MN,t∈ST,w∈S,Fn,t,wDR+,Fn,t,wDR−≥0,(26)Fn,t,wNR+−Fn,t,wNR−=Pn,t,wNR−Pn,t,1NR,∀n∈MN,t∈ST,w∈S,Fn,t,wNR+,Fn,t,wNR−≥0,(27)Fn,t,wS+−Fn,t,wS−=Pn,t,wdis−Pn,t,wch−Pn,t,1dis−Pn,t,1ch,∀n∈MN,t∈ST,w∈S,Fn,t,wS+,Fn,t,wS−≥0,(28)FwE=∑n∈MN∑t∈STFn,t,wDR++Fn,t,wDR−+Fn,t,wNR++Fn,t,wNR−+Fn,t,wS++Fn,t,wS−,∀w∈S.


## 2.1.1. Objective Function

The proposed problem in this section performs the minimization of the difference between the expected MG operation cost and the sum of expected flexibility and security benefits from the MG operator (MGO) viewpoint. This statement is expressed in (1) in a mathematical format, where the first part of this equation refers to the MG operation cost including the cost of energy received from the upstream network and NRES fuel cost [28]. The second part of (1) denotes the MG flexibility benefit that is equal to the product of the flexibility incentive price (FIP) and flexibility energy (FE) from all FSs, i.e., DRP, ESS, and NRES. Finally, the security benefit of MG is expressed in the third term of (1), depending on the security inactive price (SIP) and the security energy (SE) from all sources and active loads.In this paper, since an economic model for MGs and SDNs is considered, flexibility, voltage security, and network operation indices are formulated in the economic model as were cost and profit. Therefore, the objective function of the proposed problem considers the mentioned indices in an integrated manner. Note, however, that this type of modeling is based on multi-objective optimization based on a normalized objective function, in which all parts of the objective function will have the same dimension [20, 22, 34] similar to the objective function (1) in which all parts are in dollars.(1)min∑w∈Sπw∑t∈Tρt,wPms,t,wMS+∑n∈MNβnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Security benefit.


## 2.1.2. MG Constraints

The constraints of this section include AC power flow equations, expressed in (2)–(7), and system operation limits, described in (8)–(11) [22]. It is noteworthy that (2)–(7) represent active and reactive power balance in MG buses [44], active and reactive power from MG lines, and values of voltage angle and magnitude in the slack bus, respectively. In this model, variables PMS and QMS denote, respectively, MG station active and reactive powers, existing only in the slack bus. The value of these variables is zero in other buses. In addition, the system operation limits consist of voltage limits in MG buses (8), MG line capacity limit (9), and MG station capacity and power factor limits (10) and (11). It should be noted that limitation (11) presents an equivalent equation for the power factor limit of the MG station so that the reactive power of this station can change in the defined range in (11) by considering PFmin as the minimum allowed power factor to satisfy this limitation.(2)Pn,t,wMS+Pn,t,wR+Pn,t,wNR+Pn,t,wdis−Pn,t,wch−∑j∈MNAn,jPn,j,t,wL=Ln,t,wP−Ln,t,wDR,∀n∈MN,t∈ST,w∈S,(3)Qn,t,wMS+Qn,t,wNR+Qn,t,wS−∑j∈MNAn,jQn,j,t,wL=Ln,t,wQ,∀n∈MN,t∈ST,w∈S,(4)Pn,j,t,wL=Gn,jVn,t,w−Vn,t,wVj,t,wGn,jcosδn,t,w−δj,t,w+Bn,jsinδn,t,w−δj,t,w,∀n,j∈MN,t∈ST,w∈S,(5)Qn,j,t,wL=−Bn,jVn,t,w+Vn,t,wVj,t,wBn,jcosδn,t,w−δj,t,w−Gn,jsinδn,t,w−δj,t,w,∀n,j∈MN,t∈ST,w∈S,(6)δn,t,w=δref,∀n=ms,t∈ST,w∈S,(7)Vn,t,w=Vref,∀n=ms,t∈ST,w∈S,(8)Vmin≤Vn,t,w≤Vmax,∀n∈MN,t∈ST,w∈S,(9)Pn,j,t,wL2+Qn,j,t,wL2≤Sn,jL,max2,∀n,j∈MN,t∈ST,w∈S,(10)Pn,t,wMS2+Qn,t,wMS2≤SnG,max2,∀n=ms,t∈ST,w∈S,(11)−tancos−1PFm×Pn,t,wMS≤Qn,t,wMS≤tancos−1PFm×Pn,t,wMS,∀n=ms,t∈ST,w∈S.


## 2.1.3. DRP Constraints

The incentive model of DRP is considered in this paper, based on the fact that all loads can shift part of their consumption in peak load times to off-peak load hours, according to the electrical energy price. This statement is modeled as (12) and (13), where (12) points to the fact that the consumption energy of loads with/without DRP should be equal [45]. The limitation of the proposed DRP is also expressed in (13), where ξ defines the rate of participation of load in the DRP, which varies between 0 and 1 [45]. Finally, the DRP benefit is used in the MG operation cost model in the first part of (1), as it considers the total demanded power from the SDN, where this power depends on the operation of sources and active loads.Note that in the proposed model for DRP, consumers participating in DRP reduce their energy consumption during the hours with high energy prices and receive this amount of reduced energy during the hours with inexpensive energy from the grid. Accordingly, constraint (13) expresses the controllable range of DRP power. Constraint (12) also guarantees that all reduced energy during the hours with expensive energy prices must be provided during the hours with inexpensive energy prices. Following this DRP method, a residential consumer can be an example, such that high energy consumption loads such as laundry, ironing, exercising with exercise equipment, and other items during off-peak hours can be shifted to intervals with cheap energy prices.(12)∑t∈STLn,t,wP−Ln,t,wDR=∑t∈STLn,t,wP,∀n∈MN,w∈S,(13)−ξ×Ln,t,wP≤Ln,t,wDR≤ξ×Ln,t,wP,∀n=s,t∈ST,w∈S.


## 2.1.4. DG Constraints

Generally, DG includes two types of energy sources that are renewable and nonrenewable, i.e., RESs and NRESs. The main purpose of RESs is the generation of active power according to the standard IEEE 1547 [46–49], and the output power of RESs depends on natural phenomena such as wind speed and solar radiation. Therefore, this statement is modeled as (14) considering PR,max varying in each simulation time and scenario sample. However, the NRES model includes its capacity limit as (15) which is constant for each time and scenario [28].(14)0≤Pn,t,wR≤Pn,t,wR,max,∀n∈MN,t∈ST,w∈S,(15)Pn,t,wNR2+Qn,t,wNR2≤SnNR,max2,∀n∈MN,t∈ST,w∈S.


## 2.1.5. ESS Constraints

The operation model of ESS follows (16)–(21) [50–52], with constraints (16) and (17) referring, respectively, to discharge and charge rates of ESS considering binary variable x as ESS operation status; i.e., it acts in discharge mode if x = 1; otherwise, it operates in charge mode. Moreover, (18) and (19) calculate the stored energy in ESS at different simulation times, where it should vary in the defined range based on (20). Finally, the ESS charger capacity limit is modeled as (21) indicating a circular inequality with a radius of SS,max. In addition, the ESS operation cost/benefit as DRP benefit is considered in MG operation cost formulation.(16)0≤Pn,t,wdis≤DRn,t,wxn,t,∀n∈MN,t∈ST,w∈S,(17)0≤Pn,t,wch≤CRn,t,w1−xn,t,∀n∈MN,t∈ST,w∈S,(18)En,t,w=En,t−1,w+ηchPn,t,wch−1ηdisPn,t,wdis,∀n∈MN,t∈ST,w∈S,(19)En,t,w=Enini,∀n∈MN,t=0,w∈S,(20)Enmin≤En,t,w≤Enmax,∀n∈MN,t∈ST,w∈S,(21)Pn,t,wdis−Pn,t,wch2+Qn,t,wS2≤SnS,max2,∀n∈MN,t∈ST,w∈S.


## 2.1.6. Security Constraints

Voltage security is one of the important indices in the distribution network or low voltage system [15]. This index tends to reduce in heavily loading conditions, while it increases or improves with the rise in the injection power by all sources and active loads in MG. Generally, the voltage stability margin (VSM) method based on the voltage stability index (SI-index) is used to investigate the network voltage security [21]. This method can also determine, firstly, the weakest bus of the network with the minimum SI or the worst SI (WSI); then, the value of WSI is calculated by (22), and constraint (23) should be satisfied. In addition, voltage security will be improved if the injection power rises by local sources, i.e., RES, NRES, ESS, and DRP. Therefore, security energy is equal to the sum of all injection power in the operation horizon as shown in the following equation:(22)WSIwb,t,w=.Vwb−1,t,w4−4Vwb−1,t,w2Rwb−1,wbPwb−1,wb,t,wL+Xwb−1,wbQwb−1,wb,t,wL−4Xwb−1,wbPwb−1,wb,t,wL−Rwb−1,wbQwb−1,wb,t,wL2,∀t∈ST,w∈S,(23)WSIwb,t,w≥WSImin,∀t∈ST,w∈S,(24)SwE=∑n∈MN∑t∈STLn,t,wDR+Pn,t,wR+Pn,t,wNR+Pn,t,wdis−Pn,t,wch,∀w∈S.


## 2.1.7. Flexibility Constraints

It should be noted that flexibility in the power system is defined as “the modification of generation injection and/or consumption patterns in reaction to an external price or activation signal to provide a service within the electrical system” according to [53]. Hence, each flexible source, i.e., NRES, ESS, and DRP, includes both upward and downward flexibility services. There is upward flexibility power (F+) in scenario w for each FS if the difference between the power of FS in scenarios w and 1 (refers to the scenario that has forecasted values for uncertain parameters) is positive; otherwise, it is considered as downward flexibility power (F−) with constraints (25)–(27) [32]. Moreover, the flexibility energy is equal to the sum of upward and downward flexibility power of all FSs, as seen in (28).It is noteworthy that the flexibility quantity is equal to the ability of flexibility sources to compensate for power fluctuations of RES in different scenarios compared to its predicted value. Therefore, two upward and downward functions are obtained for flexibility sources. Upward (downward) mode occurs when the source must compensate for the RES power deficit (excess) relative to its predicted value. Since, in both modes of operation, the difference in active power for the flexibility source is measured, the quantity obtained is in watts; that is, it is a power quantity. Hence, the sum of them will represent the flexibility energy. Of course, the availability of RES sources is considered in this paper, where this is expected to happen in the future.(25)Fn,t,wDR+−Fn,t,wDR−=Ln,t,wDR−Ln,t,1DR,∀n∈MN,t∈ST,w∈S,Fn,t,wDR+,Fn,t,wDR−≥0,(26)Fn,t,wNR+−Fn,t,wNR−=Pn,t,wNR−Pn,t,1NR,∀n∈MN,t∈ST,w∈S,Fn,t,wNR+,Fn,t,wNR−≥0,(27)Fn,t,wS+−Fn,t,wS−=Pn,t,wdis−Pn,t,wch−Pn,t,1dis−Pn,t,1ch,∀n∈MN,t∈ST,w∈S,Fn,t,wS+,Fn,t,wS−≥0,(28)FwE=∑n∈MN∑t∈STFn,t,wDR++Fn,t,wDR−+Fn,t,wNR++Fn,t,wNR−+Fn,t,wS++Fn,t,wS−,∀w∈S.


## 2.2. SDN Operation (Second Layer)

The second layer of the proposed CEMS refers to the optimal scheduling management of SDN devices by the SDN operator (SDNO) according to the demand of MGs and the SDN conditions. It should be noted that the problem model in this section minimizes the difference between the expected SDN operation cost and the sum of expected flexibility and security benefits subject to the network, constraints of RESs and FSs, and security and flexibility equations. Therefore, the model of the second layer, expressed in (29)–(33), is the same as the first one. However, the difference remains in the nodal active and reactive power balance and security energy equations which include the MGs demand as (30)–(32), and the terms of ds, DN, PDS, and QDS are substituted by ms, MN, PMS, and QMS, respectively.(29)min∑w∈Sπw∑t∈Tρt,wPds,t,wDS+∑n∈SNβnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Security benefit,subject to(30)Pn,t,wDS+Pn,t,wR+Pn,t,wNR+Pn,t,wdis−Pn,t,wch−∑j∈DNAn,jPn,j,t,wL=Ln,t,wP−Ln,t,wDR+∑ms∈MLCn,msPms,t,wMS.∀n∈MN,t∈ST,w∈S,(31)Qn,t,wDS+Qn,t,wNR+Qn,t,wS−∑j∈DNAn,jQn,j,t,wL=Ln,t,wQ+∑ms∈MLCn,msQms,t,wMS,∀n∈MN,t∈ST,w∈S,(32)SwE=∑n∈DN∑t∈STLn,t,wDR+Pn,t,wR+Pn,t,wNR+Pn,t,wdis−Pn,t,wch−∑ms∈MLCn,msPms,t,wMS,∀w∈S.(33)Constraints4−23and25−28substitutingdsbyms,DNbyMN,PDSbyPMS,andQDSbyQMS.


## 3. Linear Stochastic Operation of SDN

### 3.1. The Linear Model of HFSO in the SDN

The proposed problems (1)–(28) and (29)–(33) have a non-convex NLP model due to nonlinear constraints (4), (5), (9), (10), (15), (21), and (22) and non-convex equations (4) and (5) [54–56]. Hence, this model can find a locally optimal solution in the best condition [57] and follows numerical and evolutionary methods such as Newton–Raphson technique and genetic algorithm [57]. These methods will lead to long calculation time for the HFSO in the SDN. Hence, the proposed NLP model is converted to a linear approximation model in the next step, to obtain a global optimal solution with low calculation error and time compared to the original model of the CEMS-based HFSO. Finally, the details of the proposed techniques for the linearization of different constraints are expressed as follows:(1)
The Linear Model of the Power Flow Equations. For the linearization of constraints (4) and (5) based on medium and low voltage network structure, it is assumed that the voltage angle difference across an MG or distribution line (δn—δj) is less than 0.105 radian and voltage (V) can be formulated as 1 + ΔV, where ΔV (ΔV << 1 p.u) refers to the voltage deviation. Note that the voltage magnitude should vary between 0.9 and 1.05 p.u.; hence, it is close to 1 p.u., and the proposed formulation has low calculation error. Therefore, the terms cos(δn—δj) and sin(δn—δj) are almost equal to 1 and (δn - δj), and V4, V2, and VnVj are reformulated, respectively, as 1 + 4ΔV, 1 + 2ΔV, and 1 + ΔVn + ΔVj whereas ΔV4, ΔV2, ΔVnΔVj, and ΔV×(δn - δj) are considered equal to zero due to their very small value [22].(2)
The Linear Format of Circular Limits. It is noted that constraints (9), (10), (15), and (21) are circular inequalities and they can be linearized according to a polygon approximation method based on Figure 2. Details are presented in [34]. Accordingly, each side of the polygon is a straight line, and its equation can be obtained from the tangent to the circle at a specific point as depicted in Figure 2 [34].(3)
Linearization of the WSI Equation. Table 2 presents the values of terms 2 and 3 for the WSI formulation (2), based on data of the distribution line between buses 17 and 18 for a 33-bus distribution network [58]. Accordingly, the value of term 3 is much less than that of term 2 for different values of PL, QL, and V. Moreover, the difference between the values of term 2 for voltage values of 0.95 and 1 p.u. or 0.9. and 0.95 is almost equal to 5%. Therefore, for the linearization of constraint (22), term 3 is removed from this equation. V in term 2 is assumed to be 0.95 p.u., and V4 is rewritten as 1 + 4ΔV, to obtain a linear model for WSI considering low calculation error.Figure 2
The proposed linearization method for circular inequality [34].Table 2
The values of the second and third WSI formulation in different values ofPL, QL, and V.
(PL, QL) in p.u.(0.1, 0)(0.07, 0.07)Term2323V = 0.9 p.u.0.00155 × 10−70.00182 × 10−8V = 0.95 p.u.0.001585 × 10−70.00192 × 10−8V = 1 p.u.0.001665 × 10−70.00202 × 10−8(PL, QL) in p.u.(0, 0.1)(0.07, −0.07)Term2323V = 0.9 p.u.0.00128 × 10−72 × 10−41.2 × 10−6V = 0.95 p.u.0.001268 × 10−72.1 × 10−41.2 × 10−6V = 1 p.u.0.001338 × 10−72.2 × 10−41.2 × 10−6Therefore, the linear model of the proposed strategy in the SDN can be written as follows.


#### 3.1.1. First Layer Problem Model

As seen in problems (34)–(43), the objective function of (34) is the same as (1). However, constraints (35)–(37) present, respectively, the linear equations for (4), (5), and (8) based on the first proposed linearization method. In (37), the parameter ΔVmax is equal to (Vmax–Vmin)/nl, where nl is the total number of linearization segments of the voltage magnitude term. Constraints (38)–(41) also express the linear format of circular inequalities (9), (10), (15), and (21) according to the proposed method in Figure 2. Finally, the equivalent linear equation with WSI formula (22) is formulated as (42), based on the third proposed linearization technique, and constraint (43) selects the linear formulation in the problem (1)–(28).(34)min∑w∈Sπw∑t∈Tρt,wPms,t,wMS+∑n∈MNβnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Security benefit,subject to(35)Pn,j,t,wL=Gn,jΔVn,t,w−ΔVj,t,w−Bn,jδn,t,w−δj,t,w,∀n,j∈MN,t∈ST,w∈S,(36)Qn,j,t,wL=−Bn,jΔVn,t,w−ΔVj,t,w−Gn,jδn,t,w−δj,t,w,∀n,j∈MN,t∈ST,w∈S,(37)Vmin−1≤ΔVn,t,w≤Vmax−1,∀n∈MN,t∈ST,w∈S,(38)Pn,j,t,wLcosk×Δα+Qn,j,t,wLsink×Δα≤Sn,jL,max,∀n,j∈MN,t∈ST,w∈S,k∈K=1,2,…,nk,Δα=2πnk,(39)Pn,t,wMScosk×Δα+Qn,t,wMSsink×Δα≤SnG,max,∀n=ms,t∈ST,w∈S,k∈K,(40)Pn,t,wNRcosk×Δα+Qn,t,wNRsink×Δα≤SnNR,max,∀n∈MN,t∈ST,w∈S,k∈K,(41)Pn,t,wdis−Pn,t,wchcosk×Δα+Qn,t,wSsink×Δα≤SnS,max,∀n∈MN,t∈ST,w∈S,k∈K,(42)WSIwb,t,w=1+4ΔVwb,t,w−40.952Rwb−1,wbPwb−1,wb,t,wL+Xwb−1,wbQwb−1,wb,t,wL,∀t∈ST,w∈S.(43)Constraints2,3,6,7,11−14,16−20,23−28.As the proposed NLP and LP models can generally obtain the local optimal solution and the model status is not the same for both models, it is not reasonable to compare them to determine the calculation error. Therefore, the power flow analysis based on the NLP and LP formulae for the 33-bus distribution network [58] is presented in this section to obtain the calculation error of different variables. The decision variables in this problem according to the two proposed models include constant values [22], and it is also considered that the proposed linearization methods of 1 and 3 can result in considerable calculation error if the proposed linear format is not suitable. However, the computation error due to the proposed linearization model 2 can reach a low value if more sides are used for the polygon. Finally, based on the expressed results in Table 3, the computation error for active and reactive power is, respectively, about 0.5% and 0.6% and is close, respectively, to 0.3% and 0.4% for voltage magnitude as well as voltage angle and WSI in the LP format considering the NLP model. Moreover, the LP’s calculation time is much less than NLP’s computational time. Therefore, the proposed NLP model can be substituted by the proposed LP method due to low calculation time and error in the LP format.Table 3
Comparison of AC power flow analysis in the NLP and LP model.
ModelNLPLPDeviation (%)Station active power (p.u.)3.9173.8970.51Station reactive power (p.u.)2.4352.4210.57Mean of voltage magnitude (p.u.)0.9490.9520.31Mean of voltage angle (rad)4.918 × 10−54.936 × 10−50.38WSI (without unit)0.6960.6930.42Calculation time (s)2.9880.731—


#### 3.1.2. Second Layer Problem Model

The linear programming of this section is presented in (44)–(46), including the objective function of (29) with the same equation as (44). Moreover, the constraints of the linear HFSO in the SDN contain linear equations (29)–(33) and (34)–(43) that are used in (45) and (46), respectively.(44)min∑w∈Sπw∑t∈Tρt,wPms,t,wDS+∑n∈MNβnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Security benefit,subject to(45)Constraints30−32.(46)Constraints6,7,11−14,16−20,23,25−28,35−42with substitutingdsbyms,DNbyMN,PDSbyPMS,andQDSbyQMS.


### 3.2. Stochastic CEMS-Based HFSO

In the proposed two-layer CEMS for HFSO of SDN, the parameters of active and reactive load,LP and LQ; energy price, ρ; and maximum active power of RES, PR,max, are considered uncertainties. Therefore, scenario-based stochastic programming (SBSP) is used in this paper to model these parameters based on Roulette Wheel Mechanism (RWM) and to generate scenario samples. The Kantorovich method is applied to reduce the generated scenario samples by RWM with a higher probability in comparison with removed scenarios [59]. The RWM generates scenarios according to normal distribution for load and price forecasting error and based on beta/Weibull distribution, which is suitable for RES type of solar/wind. Finally, more details of the proposed method are presented in [59].


### 3.3. Solution Method

In the proposed two-layer CEMS of the SDN, in this paper, the second layer problem depends on the values ofPMS and QMS. Furthermore, the optimal values of PMS and QMS obtained by the first layer problem depend on SDN operation limits. Hence, the master/slave decomposition (MSD) method is used in this section to obtain the optimal solution for the proposed strategy at a low calculation time [60]. The MMG operation problem is solved in the master problem (MP), and SDN scheduling is investigated in the slave problem (SP). Finally, the algorithm process is summarized as follows:(i)
Step 1 (initial MP): In this section, problem (34)–(43) is solved to determine the optimal values of PMS and QMS, PMS⌢v and QMS^v. v is the iteration and is equal to 1 in this step, and the symbol “ ^ ” represents the optimal value of a variable.(ii)
Step 2 (SP): This problem is formulated as follows:(47)min∑w∈Sπw∑t∈Tρt,wPms,t,wDS+∑n∈,βnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Security benefit+ω⋅x−x^,subject to(48)Constraints45−46.Model (47) and (48) is the same as the problem (44)–(46), except that ω⋅x−x^ is added to the objective function (44). In (47), x represents the variables PMS and QMS, and x^ is the optimal value of these variables, which can be calculated from MP. In addition, ω is a constant coefficient, i.e., 103. The SP may not reach the optimal solution without the expression ω⋅x−x^ forx^ values due to considering SDN operation limits. To compensate for this, the phraseω⋅x−x^ is added to the objective function of the second layer problem [60]. Therefore, in this case, the determined value of x in this step is equal to the optimal values of PMS and QMS in iteration v. Note also that the expression ω⋅x−x^ is nonlinear. To linearize it, the term ω⋅x++x− replacesω.x−x^ in (47), where x+ and x− are auxiliary variables with constraints x++x−=x−x^, 0≤x+≤M.z, and 0≤x−≤M⋅1−z being determined by problem (47) and (48). Moreover, M is a large fixed number, i.e., 106, and z is an auxiliary binary variable. Finally, the optimal values of the variables calculated in this step are denoted by the symbol “∼.”(iii)
Step 3 (convergence analysis): The problems proposed in Section 3.1 meet the convergence conditions if the following equation is established:(49)maxx˜v−x^v,∀n,t,w≤ε,∀x=PMSQMS,whereε is the permissible computational error or tolerance of the MSD method. It must be said that if the relationship is not established, Step 4 must be performed, and v tends to v + 1.(iv)
Step 4 (MP): This problem is formulated as follows:(50)min∑w∈Sπw∑t∈Tρt,wPms,t,wMS+∑n∈MNβnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Securitybenefit.Subject to(51)Constraints35−43,(52)x≤x˜i,∀i=1,2,…,v−1.Equations (50) and (51) are initial MP, and (52) is a constraint based on previous SP results, which can be called “MSD cut.” Finally, the flowchart of the MSD for the proposed problem is presented in Figure 3.Figure 3
MSD flowchart to solve the proposed problem.


## 3.1. The Linear Model of HFSO in the SDN

The proposed problems (1)–(28) and (29)–(33) have a non-convex NLP model due to nonlinear constraints (4), (5), (9), (10), (15), (21), and (22) and non-convex equations (4) and (5) [54–56]. Hence, this model can find a locally optimal solution in the best condition [57] and follows numerical and evolutionary methods such as Newton–Raphson technique and genetic algorithm [57]. These methods will lead to long calculation time for the HFSO in the SDN. Hence, the proposed NLP model is converted to a linear approximation model in the next step, to obtain a global optimal solution with low calculation error and time compared to the original model of the CEMS-based HFSO. Finally, the details of the proposed techniques for the linearization of different constraints are expressed as follows:(1)
The Linear Model of the Power Flow Equations. For the linearization of constraints (4) and (5) based on medium and low voltage network structure, it is assumed that the voltage angle difference across an MG or distribution line (δn—δj) is less than 0.105 radian and voltage (V) can be formulated as 1 + ΔV, where ΔV (ΔV << 1 p.u) refers to the voltage deviation. Note that the voltage magnitude should vary between 0.9 and 1.05 p.u.; hence, it is close to 1 p.u., and the proposed formulation has low calculation error. Therefore, the terms cos(δn—δj) and sin(δn—δj) are almost equal to 1 and (δn - δj), and V4, V2, and VnVj are reformulated, respectively, as 1 + 4ΔV, 1 + 2ΔV, and 1 + ΔVn + ΔVj whereas ΔV4, ΔV2, ΔVnΔVj, and ΔV×(δn - δj) are considered equal to zero due to their very small value [22].(2)
The Linear Format of Circular Limits. It is noted that constraints (9), (10), (15), and (21) are circular inequalities and they can be linearized according to a polygon approximation method based on Figure 2. Details are presented in [34]. Accordingly, each side of the polygon is a straight line, and its equation can be obtained from the tangent to the circle at a specific point as depicted in Figure 2 [34].(3)
Linearization of the WSI Equation. Table 2 presents the values of terms 2 and 3 for the WSI formulation (2), based on data of the distribution line between buses 17 and 18 for a 33-bus distribution network [58]. Accordingly, the value of term 3 is much less than that of term 2 for different values of PL, QL, and V. Moreover, the difference between the values of term 2 for voltage values of 0.95 and 1 p.u. or 0.9. and 0.95 is almost equal to 5%. Therefore, for the linearization of constraint (22), term 3 is removed from this equation. V in term 2 is assumed to be 0.95 p.u., and V4 is rewritten as 1 + 4ΔV, to obtain a linear model for WSI considering low calculation error.Figure 2
The proposed linearization method for circular inequality [34].Table 2
The values of the second and third WSI formulation in different values ofPL, QL, and V.
(PL, QL) in p.u.(0.1, 0)(0.07, 0.07)Term2323V = 0.9 p.u.0.00155 × 10−70.00182 × 10−8V = 0.95 p.u.0.001585 × 10−70.00192 × 10−8V = 1 p.u.0.001665 × 10−70.00202 × 10−8(PL, QL) in p.u.(0, 0.1)(0.07, −0.07)Term2323V = 0.9 p.u.0.00128 × 10−72 × 10−41.2 × 10−6V = 0.95 p.u.0.001268 × 10−72.1 × 10−41.2 × 10−6V = 1 p.u.0.001338 × 10−72.2 × 10−41.2 × 10−6Therefore, the linear model of the proposed strategy in the SDN can be written as follows.


### 3.1.1. First Layer Problem Model

As seen in problems (34)–(43), the objective function of (34) is the same as (1). However, constraints (35)–(37) present, respectively, the linear equations for (4), (5), and (8) based on the first proposed linearization method. In (37), the parameter ΔVmax is equal to (Vmax–Vmin)/nl, where nl is the total number of linearization segments of the voltage magnitude term. Constraints (38)–(41) also express the linear format of circular inequalities (9), (10), (15), and (21) according to the proposed method in Figure 2. Finally, the equivalent linear equation with WSI formula (22) is formulated as (42), based on the third proposed linearization technique, and constraint (43) selects the linear formulation in the problem (1)–(28).(34)min∑w∈Sπw∑t∈Tρt,wPms,t,wMS+∑n∈MNβnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Security benefit,subject to(35)Pn,j,t,wL=Gn,jΔVn,t,w−ΔVj,t,w−Bn,jδn,t,w−δj,t,w,∀n,j∈MN,t∈ST,w∈S,(36)Qn,j,t,wL=−Bn,jΔVn,t,w−ΔVj,t,w−Gn,jδn,t,w−δj,t,w,∀n,j∈MN,t∈ST,w∈S,(37)Vmin−1≤ΔVn,t,w≤Vmax−1,∀n∈MN,t∈ST,w∈S,(38)Pn,j,t,wLcosk×Δα+Qn,j,t,wLsink×Δα≤Sn,jL,max,∀n,j∈MN,t∈ST,w∈S,k∈K=1,2,…,nk,Δα=2πnk,(39)Pn,t,wMScosk×Δα+Qn,t,wMSsink×Δα≤SnG,max,∀n=ms,t∈ST,w∈S,k∈K,(40)Pn,t,wNRcosk×Δα+Qn,t,wNRsink×Δα≤SnNR,max,∀n∈MN,t∈ST,w∈S,k∈K,(41)Pn,t,wdis−Pn,t,wchcosk×Δα+Qn,t,wSsink×Δα≤SnS,max,∀n∈MN,t∈ST,w∈S,k∈K,(42)WSIwb,t,w=1+4ΔVwb,t,w−40.952Rwb−1,wbPwb−1,wb,t,wL+Xwb−1,wbQwb−1,wb,t,wL,∀t∈ST,w∈S.(43)Constraints2,3,6,7,11−14,16−20,23−28.As the proposed NLP and LP models can generally obtain the local optimal solution and the model status is not the same for both models, it is not reasonable to compare them to determine the calculation error. Therefore, the power flow analysis based on the NLP and LP formulae for the 33-bus distribution network [58] is presented in this section to obtain the calculation error of different variables. The decision variables in this problem according to the two proposed models include constant values [22], and it is also considered that the proposed linearization methods of 1 and 3 can result in considerable calculation error if the proposed linear format is not suitable. However, the computation error due to the proposed linearization model 2 can reach a low value if more sides are used for the polygon. Finally, based on the expressed results in Table 3, the computation error for active and reactive power is, respectively, about 0.5% and 0.6% and is close, respectively, to 0.3% and 0.4% for voltage magnitude as well as voltage angle and WSI in the LP format considering the NLP model. Moreover, the LP’s calculation time is much less than NLP’s computational time. Therefore, the proposed NLP model can be substituted by the proposed LP method due to low calculation time and error in the LP format.Table 3
Comparison of AC power flow analysis in the NLP and LP model.
ModelNLPLPDeviation (%)Station active power (p.u.)3.9173.8970.51Station reactive power (p.u.)2.4352.4210.57Mean of voltage magnitude (p.u.)0.9490.9520.31Mean of voltage angle (rad)4.918 × 10−54.936 × 10−50.38WSI (without unit)0.6960.6930.42Calculation time (s)2.9880.731—


### 3.1.2. Second Layer Problem Model

The linear programming of this section is presented in (44)–(46), including the objective function of (29) with the same equation as (44). Moreover, the constraints of the linear HFSO in the SDN contain linear equations (29)–(33) and (34)–(43) that are used in (45) and (46), respectively.(44)min∑w∈Sπw∑t∈Tρt,wPms,t,wDS+∑n∈MNβnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Security benefit,subject to(45)Constraints30−32.(46)Constraints6,7,11−14,16−20,23,25−28,35−42with substitutingdsbyms,DNbyMN,PDSbyPMS,andQDSbyQMS.


## 3.1.1. First Layer Problem Model

As seen in problems (34)–(43), the objective function of (34) is the same as (1). However, constraints (35)–(37) present, respectively, the linear equations for (4), (5), and (8) based on the first proposed linearization method. In (37), the parameter ΔVmax is equal to (Vmax–Vmin)/nl, where nl is the total number of linearization segments of the voltage magnitude term. Constraints (38)–(41) also express the linear format of circular inequalities (9), (10), (15), and (21) according to the proposed method in Figure 2. Finally, the equivalent linear equation with WSI formula (22) is formulated as (42), based on the third proposed linearization technique, and constraint (43) selects the linear formulation in the problem (1)–(28).(34)min∑w∈Sπw∑t∈Tρt,wPms,t,wMS+∑n∈MNβnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Security benefit,subject to(35)Pn,j,t,wL=Gn,jΔVn,t,w−ΔVj,t,w−Bn,jδn,t,w−δj,t,w,∀n,j∈MN,t∈ST,w∈S,(36)Qn,j,t,wL=−Bn,jΔVn,t,w−ΔVj,t,w−Gn,jδn,t,w−δj,t,w,∀n,j∈MN,t∈ST,w∈S,(37)Vmin−1≤ΔVn,t,w≤Vmax−1,∀n∈MN,t∈ST,w∈S,(38)Pn,j,t,wLcosk×Δα+Qn,j,t,wLsink×Δα≤Sn,jL,max,∀n,j∈MN,t∈ST,w∈S,k∈K=1,2,…,nk,Δα=2πnk,(39)Pn,t,wMScosk×Δα+Qn,t,wMSsink×Δα≤SnG,max,∀n=ms,t∈ST,w∈S,k∈K,(40)Pn,t,wNRcosk×Δα+Qn,t,wNRsink×Δα≤SnNR,max,∀n∈MN,t∈ST,w∈S,k∈K,(41)Pn,t,wdis−Pn,t,wchcosk×Δα+Qn,t,wSsink×Δα≤SnS,max,∀n∈MN,t∈ST,w∈S,k∈K,(42)WSIwb,t,w=1+4ΔVwb,t,w−40.952Rwb−1,wbPwb−1,wb,t,wL+Xwb−1,wbQwb−1,wb,t,wL,∀t∈ST,w∈S.(43)Constraints2,3,6,7,11−14,16−20,23−28.As the proposed NLP and LP models can generally obtain the local optimal solution and the model status is not the same for both models, it is not reasonable to compare them to determine the calculation error. Therefore, the power flow analysis based on the NLP and LP formulae for the 33-bus distribution network [58] is presented in this section to obtain the calculation error of different variables. The decision variables in this problem according to the two proposed models include constant values [22], and it is also considered that the proposed linearization methods of 1 and 3 can result in considerable calculation error if the proposed linear format is not suitable. However, the computation error due to the proposed linearization model 2 can reach a low value if more sides are used for the polygon. Finally, based on the expressed results in Table 3, the computation error for active and reactive power is, respectively, about 0.5% and 0.6% and is close, respectively, to 0.3% and 0.4% for voltage magnitude as well as voltage angle and WSI in the LP format considering the NLP model. Moreover, the LP’s calculation time is much less than NLP’s computational time. Therefore, the proposed NLP model can be substituted by the proposed LP method due to low calculation time and error in the LP format.Table 3
Comparison of AC power flow analysis in the NLP and LP model.
ModelNLPLPDeviation (%)Station active power (p.u.)3.9173.8970.51Station reactive power (p.u.)2.4352.4210.57Mean of voltage magnitude (p.u.)0.9490.9520.31Mean of voltage angle (rad)4.918 × 10−54.936 × 10−50.38WSI (without unit)0.6960.6930.42Calculation time (s)2.9880.731—


## 3.1.2. Second Layer Problem Model

The linear programming of this section is presented in (44)–(46), including the objective function of (29) with the same equation as (44). Moreover, the constraints of the linear HFSO in the SDN contain linear equations (29)–(33) and (34)–(43) that are used in (45) and (46), respectively.(44)min∑w∈Sπw∑t∈Tρt,wPms,t,wDS+∑n∈MNβnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Security benefit,subject to(45)Constraints30−32.(46)Constraints6,7,11−14,16−20,23,25−28,35−42with substitutingdsbyms,DNbyMN,PDSbyPMS,andQDSbyQMS.


## 3.2. Stochastic CEMS-Based HFSO

In the proposed two-layer CEMS for HFSO of SDN, the parameters of active and reactive load,LP and LQ; energy price, ρ; and maximum active power of RES, PR,max, are considered uncertainties. Therefore, scenario-based stochastic programming (SBSP) is used in this paper to model these parameters based on Roulette Wheel Mechanism (RWM) and to generate scenario samples. The Kantorovich method is applied to reduce the generated scenario samples by RWM with a higher probability in comparison with removed scenarios [59]. The RWM generates scenarios according to normal distribution for load and price forecasting error and based on beta/Weibull distribution, which is suitable for RES type of solar/wind. Finally, more details of the proposed method are presented in [59].


## 3.3. Solution Method

In the proposed two-layer CEMS of the SDN, in this paper, the second layer problem depends on the values ofPMS and QMS. Furthermore, the optimal values of PMS and QMS obtained by the first layer problem depend on SDN operation limits. Hence, the master/slave decomposition (MSD) method is used in this section to obtain the optimal solution for the proposed strategy at a low calculation time [60]. The MMG operation problem is solved in the master problem (MP), and SDN scheduling is investigated in the slave problem (SP). Finally, the algorithm process is summarized as follows:(i)
Step 1 (initial MP): In this section, problem (34)–(43) is solved to determine the optimal values of PMS and QMS, PMS⌢v and QMS^v. v is the iteration and is equal to 1 in this step, and the symbol “ ^ ” represents the optimal value of a variable.(ii)
Step 2 (SP): This problem is formulated as follows:(47)min∑w∈Sπw∑t∈Tρt,wPms,t,wDS+∑n∈,βnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Security benefit+ω⋅x−x^,subject to(48)Constraints45−46.Model (47) and (48) is the same as the problem (44)–(46), except that ω⋅x−x^ is added to the objective function (44). In (47), x represents the variables PMS and QMS, and x^ is the optimal value of these variables, which can be calculated from MP. In addition, ω is a constant coefficient, i.e., 103. The SP may not reach the optimal solution without the expression ω⋅x−x^ forx^ values due to considering SDN operation limits. To compensate for this, the phraseω⋅x−x^ is added to the objective function of the second layer problem [60]. Therefore, in this case, the determined value of x in this step is equal to the optimal values of PMS and QMS in iteration v. Note also that the expression ω⋅x−x^ is nonlinear. To linearize it, the term ω⋅x++x− replacesω.x−x^ in (47), where x+ and x− are auxiliary variables with constraints x++x−=x−x^, 0≤x+≤M.z, and 0≤x−≤M⋅1−z being determined by problem (47) and (48). Moreover, M is a large fixed number, i.e., 106, and z is an auxiliary binary variable. Finally, the optimal values of the variables calculated in this step are denoted by the symbol “∼.”(iii)
Step 3 (convergence analysis): The problems proposed in Section 3.1 meet the convergence conditions if the following equation is established:(49)maxx˜v−x^v,∀n,t,w≤ε,∀x=PMSQMS,whereε is the permissible computational error or tolerance of the MSD method. It must be said that if the relationship is not established, Step 4 must be performed, and v tends to v + 1.(iv)
Step 4 (MP): This problem is formulated as follows:(50)min∑w∈Sπw∑t∈Tρt,wPms,t,wMS+∑n∈MNβnPn,t,wNR︷Operation cost−FIP×FwE︷Flexibility benefit−SIP×SwE︷Securitybenefit.Subject to(51)Constraints35−43,(52)x≤x˜i,∀i=1,2,…,v−1.Equations (50) and (51) are initial MP, and (52) is a constraint based on previous SP results, which can be called “MSD cut.” Finally, the flowchart of the MSD for the proposed problem is presented in Figure 3.Figure 3
MSD flowchart to solve the proposed problem.


## 4. Numerical Results and Discussion

### 4.1. Case Study

The proposed two-layer CEMS-based HFSO is applied on the 33-bus distribution network shown in Figure4 [58]. This network includes radial MG1 (15-bus), MG2 (13-bus), and MG3 (14-bus) located, respectively, on buses 22, 14, and 33. One should note that the line, station, and peak load data of these networks are presented in [28], and the load value for simulation hours is equal to the multiplication of the peak load value and daily load factor curve plotted in Figure 5 [61]. In addition, there are several NRES diesel generators, and the RES includes photovoltaic and wind systems, battery energy storage in the proposed SDN and MGs, their location and capacity, and other data described in [28]. The daily forecasted power percentage of RESs is based on Figure 5 [61]. In this paper, it is assumed that the power factor of MGs load is 0.9, and all loads can participate in the DRP by considering ξ = 0.3. Finally, the daily energy price is presented in [20], and the range of voltage, flexibility and security incentive prices, and minimum allowed value of WSI are assumed, respectively, to be [0.9 p.u, 1.1 p.u] [62], 10 and 10 $/MWh, and 0.8.Figure 4
The scheme of the 33-bus distribution network including MMG [58].Figure 5
Forecasted power percentage of RES and network load [61].In the SDN, the slack bus (bus 1) is assumed to have a voltage magnitude of 1 p.u., and the voltage angle is zero. Each MG also consists of sources, ESSs, and DRP, which are expected to be managed so that each MG in the bus connected to the SDN can make the voltage magnitude of the connecting bus equal to 1 p.u. Thus, buses 14, 22, and 33, on which MG2, MG1, and MG3 are placed, respectively, are considered as the slack buses for the mentioned MGs, with a voltage magnitude of 1 p.u. However, the voltage angle of these buses is equal to that obtained from the calculations in the proposed scheme.


### 4.2. Results

The proposed linear form of the two-layer CEMS-based HFSO is simulated in GAMS23.5.2, and thus it is solved by the CPLEX method [63] by considering 180 linearization segments for circular constraints. The RWM generates 1000 scenario samples, and then the Kantorovich method reduces this number to 20 scenario samples with high probability.


#### 4.2.1. Comparison of the convergence results of the proposed scheme for linear and nonlinear models

The MSD convergence for the proposed scheme in the MINLP and MILP models is investigated in Table4. In the MINLP model, BARON, BONMIN, DISOPT, KNITRO, and OQNLP solvers are used [63]. According to Table 4, among these solvers, DISOPT and OQNLP failed to find the optimal solution to the proposed scheme. Other algorithms also have different optimum points, meaning that the values of their objective functions are not the same. Therefore, MINLP solvers cannot achieve the unique optimal solution, which reduces the reliability. However, among the mentioned algorithms, it can be seen that BONMIN has a more desirable situation than other solvers because it has the shortest computational time (3027.11 s) and fewest convergence iterations for the MSD (12 iterations), and calculates the minimum objective function ($2621.73). Then, BONMIN, CBC, and CPLEX solvers are used to solve the proposed MILP model [42]. According to Table 4, it can be seen that these algorithms obtain a unique optimal solution; i.e., they reach a constant value for the objective function in the three algorithms, in the least number of MSD convergence iterations (8 iterations). Thus, these results have higher reliability. They also have a lower computational time than solvers in the MINLP model, where CPLEX has the shortest possible time (102.68 s). Therefore, the MILP model with CPLEX solver is suitable for solving the proposed problem in MSD format (Section 3.3), which has low computational time and low computational error based on Section 3.1.Table 4
Convergence results of the proposed problem based on MSD obtained by different solvers.
ModelSolverObjective function ($) in problem ofConvergence iteration forCalculation timesModel statusFirst layerSecond layerFirst layer problemSecond layer problemMSDMINLPBARON2702.651589.349284143273.43Locally optimalBONMIN2621.731537.617568123027.11Locally optimalDISOPT——————InfeasibleKNITRO2768.281611.5910396173347.38Locally optimalOQNLP——————InfeasibleMILPBONMIN2182.8861275.2274698129.57OptimalCBC2182.8861275.2264608113.83OptimalCPLEX2182.8861275.2255518102.68Optimal


#### 4.2.2. Optimal Scheduling of MG Devices

The results of the HFSO in each MG based on the first layer model of CEMS are presented in this section to determine the optimal scheduling of MG sources and active loads. All these numerical results are plotted in Figure6. Accordingly, the RESs, i.e., photovoltaic and wind systems, inject their maximum active power into MGs due to their zero operation cost so that the power level of total RESs in the MG1 and MG2 is generally less than network load in simulation for all hours. However, it is close to the total active load at the periods of 6 : 00–13 : 00 and 4 : 00–13 : 00 for MG1 and MG2, respectively. RES active power is also close to/greater than the network load at the period of 4 : 00–6:00/6 : 00–16 : 00 for MG3. Moreover, the NRESs, i.e., diesel generators existing in MG1 and MG2, are generally switched off from 1 : 00 to 6 : 00, because the fuel price of NRESs is greater than the energy price in this period. But, in the other hours, it injects high active power into MGs due to the small fuel price considering energy price and large capacity in comparison with other flexible sources. Finally, the active loads, i.e., ESS and DRP, are charged/discharged in the low/high energy price hours or period or 1 : 00–7:00/18 : 00–22 : 00, to obtain minimum energy cost for MGs. However, they are discharged at some hours of the period of medium energy price or 8 : 00–17 : 00 and 23 : 00–00 : 00, to achieve the maximum security energy or benefit. These loads are charged at other times of this period to provide the required energy for their discharge operation. As the last point, the MG devices operate according to Figure 6 to achieve the minimum energy cost received from the upstream network or SDN and to obtain maximum security and flexibility for MG by satisfying the MG operation limits such as voltage, capacity of MG lines, and power factor limitations.Figure 6
Daily expected total active power curve of different devices in (a) MG1, (b) MG2, and (c) MG3.
(a)(b)(c)


#### 4.2.3. Optimal Operation of the SDN

The daily expected active and reactive power of MGs absorbed/injected from/into the SDN according to the first layer CEMS-based HFSO are shown in Figure7. Based on Figure 7(a), MG1 receives the active power from the SDN at the period 1 : 00–7:00 due to low energy price in this interval, but it injects 3 MW active power into the SDN in other hours, arising from high injection power of MG local sources, i.e., NRES, RES, ESS, and DRP, according to Figure 6(a). However, the power generation capacity of these sources in MG2 and MG3 is lower than that in MG1. Hence, these MGs cannot inject the active power into the SDN at more simulation times as seen in Figure 7(a). Moreover, the ESS and NRES can provide the reactive load of MG1. Thus, MG1 does not receive/deliver reactive power from/to the SDN to obtain unit power factor and suitable regulation of bus voltages based on Figure 7(b). However, this condition cannot happen in MG2 and MG3, due to the low reactive power generation capacity of ESS and NRES in these MGs. Therefore, they absorb the reactive power from the SDN regarding Figure 7(b).Figure 7
Daily expected power curve of different MGs in the SDN: (a) active power; (b) reactive power.
(a)(b)In addition, the optimal operation of the SDN devices based on HFSO matching the second layer CEMS is illustrated in Figure8, where the same result as Figure 6 can be observed to obtain the minimum cost of energy received from the upstream grid or the sub-transmission network and high flexibility and security to the SDN subject to SDN operation limits, (8)–(11). Figure 9 shows the daily scheduling active and reactive power of the SDN station regarding the proposed strategy.Figure 8
Daily expected total active power curve of different devices in the SDN.Figure 9
Daily expected active/reactive power curve of the SDN station.As can be seen in Figure9, the SDN absorbs the active power from the sub-transmission network from 1 : 00 to 7 : 00 due to low energy prices at these hours. However, it has been able to inject active power of 1 to 2.5 MW into the upstream network at other hours as MGs and SDN devices deliver high active power to the SDN in the period of 8 : 00–00 : 00 according to Figures 7(a) and 8. Moreover, since the reactive power generation capacity of ESS and NRES in the SDN is low, this network cannot obtain a unity power factor and generally receives reactive power from the upstream grid based on Figure 9.


#### 4.2.4. Capabilities of HFSO-based two-layer CEMS

This section considers the following cases to investigate the capability of the proposed scheme:(i)
Case I: power flow analysis.(ii)
Case II: proposed scheme considering SIP and FIP = 0.(iii)
Case III: proposed scheme considering SIP = 0 and FIP = 10 $/MWh.(iv)
Case IV: proposed scheme considering SIP = 10 $/MWh and FIP = 0.(v)
Case V: proposed scheme considering SIP and FIP = 10 $/MWh.Tables5 and 6 express, respectively, the technical and economic capabilities of the proposed strategy. Based on Table 5, there is high energy loss in power flow studies, in which the presence of sources, ESSs, and DRP is not taken into account in MGs and SDN. According to Figures 5 and 7, the total energy demand of all loads in MG1–3 and SDN is 126.4364 MWh, 121.8664 MWh, 83.7832 MWh, and 56.6 MWh, respectively. Therefore, energy loss in the networks in Case I is roughly 5.35%, 4.83%, 2.15%, and 4.85% of their consumption energy. Nevertheless, by proper management of sources, ESSs, and DRP in other study cases (II–V), energy loss reduces compared to Case I. In the proposed scheme (Case V), energy loss in MG1–3 and SDN has been reduced by about 32.9% ((6.76–4.533)/6.76), 31%, 36%, and 32% compared to Case I. Such conditions are also true for the maximum voltage drop of the network. The maximum voltage drop in Case I for the mentioned MGs is higher than the permissible range of 0.1 p.u. (1–0.9). However, its value in the proposed scheme (Case V) has been enhanced by about 22.5%, 63.2%, 80.9%, and 49.5% compared to Case I for MG1–3 and SDN. Concerning voltage security, the mentioned networks in Case I have the sum value of WSI close to their minimum value of 19.2 (24 × 0.8). Yet the proposed scheme succeeds to improve this index by about 14.66%, 5.18%, 16.11%, and 19.37% for MG1–3 and SDN. It is worth noting that Case V in terms of different indices enjoys a lower percentage of improvement compared to Cases II-IV. For instance, in terms of energy loss, Case II provides more suitable results than Case V because, in this case, minimization of the expected cost of energy purchase from the upstream network and the expected operation cost of sources is only considered. To reduce energy cost, energy consumption and energy loss need to be small in the network as per (1). Hence, in Case II, less energy loss resulted compared to Case V in which the economic status of the operation, flexibility, and network security indices is taken into account simultaneously. Concerning voltage security and maximum voltage drop, Case IV provides more desirable results than Case V because objective function (1) in this case study only considers the improvement of operation and network security statuses. Consequently, network operators encourage sources, ESSs, and DRP to inject more energy into MGs and SDN. As a result, by supplying more security energy (SE) than other cases, this case succeeds to obtain a more suitable WSI and voltage drop than Case V and even other case studies. However, it should be noted that these circumstances are reached for an increase in the maximum overvoltage of the networks, but their values are smaller than the permissible limit of 0.1 p.u. (1.1−1). Concerning flexibility, Case III is superior to Case V because it takes into account the improvement of operation and flexibility indices in the objective function (1). By supplying more flexibility energy (FE), Case III provides more flexibility than other case studies. Furthermore, Cases II-V obtain suitable power factors for distribution substations of different networks so that the distribution substation in MG1 has always a power factor of unity. This parameter in Case I is smaller than the allowable value of 0.9 p.u. In the following, as seen in Table 5, MG1, MG2, and SDN have achieved suitable WSI by obtaining higher security energy (SE). However, MG3 needs a low SE to reach the suitable WSI. This fact is due to the high voltage drop in MG1, MG2, and SDN compared to MG3. Hence, MG3 elements inject low SE into MG3 compared to other networks. Finally, it should be noted that the active power generation capacity of flexible sources, i.e., NRES, ESS, and DRP, in MG1 is higher than that of other grids based on Section 4.1. Therefore, the FE in this network is greater than that of MG1, MG2, and SDN.Table 5
The values of different technical indices in MGs and SDN.
MG1IndexCase ICase IICase IIICase IVCase VOperationDaily energy loss (MWh)6.764.434.644.614.533Maximum voltage drop/overvoltage (p.u.)0.102/00.071/0.0450.080/0.0430.070/0.050.079/0.046Minimum power factor of station0.861111SecurityPoor voltage bus1413131313Sum of WSI19.3322.02121.98522.21222.164SE (MWh)—155.26154.11168.47164.013FlexibilityFE (MWh)—54.28165.84553.49663.976MG2IndexCase ICase IICase IIICase IVCase VOperationDaily energy loss (MWh)5.893.9744.1124.1064.067Maximum voltage drop/overvoltage (p.u.)0.143/00.051/0.0360.055/0.0350.050/0.0400.054/0.038Minimum power factor of station0.880.900.900.900.90SecurityPoor voltage bus1212121212Sum of WSI19.2720.09120.03520.31420.268SE (MWh)—117.562116.792126.113123.849FlexibilityFE (MWh)—25.19233.45824.27730.953MG3IndexCase ICase IICase IIICase IVCase VOperationDaily energy loss (MWh)1.8061.0581.2371.2181.156Maximum voltage drop/overvoltage (p.u.)0.136/00.022/0.0150.028/0.0130.021/0.0170.026/0.016Minimum power factor of station0.870.900.900.900.90SecurityPoor voltage bus1414141414Sum of WSI19.4822.17522.05722.84622.619SE (MWh)—67.31566.10478.97475.94FlexibilityFE (MWh)—16.3321.7215.5119.60SDNIndexCase ICase IICase IIICase IVCase VOperationDaily energy loss (MWh)2.751.731.9741.9521.87Maximum voltage drop/overvoltage (p.u.)0.097/00.046/0.0390.051/0.0380.045/0.0420.049/0.041Minimum power factor of station0.830.930.920.920.92SecurityPoor voltage bus1831313131Sum of WSI19.6722.9822.8323.7123.48SE (MWh)—88.37186.92299.21396.975FlexibilityFE (MWh)—14.54518.38713.79416.311Table 6
The values of different economic indices in MGs and SDN.
CaseIIIIIIIVVMG1Energy cost ($)2836.11−1145.37−1094.22−1109.23−1090.732NRES operation cost ($)02638.122678.432704.522749.906Security benefit ($)—001684.71640.13Flexibility benefit ($)—0658.450639.76Profit ($)−2836.11−1492.75−925.7689.41620.716MG2Energy cost ($)1022.43−213.18−200.33−204.27−198.324NRES operation cost ($)0982.341001.121015.461032.385Security benefit ($)—001261.131238.49Flexibility benefit ($)—0334.580309.53Profit ($)−1022.43−769.16−466.21449.94713.96MG3Energy cost ($)473.6399.45105.08103.96106.79NRES operation cost ($)00000Security benefit ($)—00789.74759.40Flexibility benefit ($)—0217.20195.60Profit ($)−473.63−99.45112.12685.78848.21SDNEnergy cost ($)1347.41−1085.29−1030.35−1045.74−1022.362NRES operation cost ($)0812.31833.64856.27880Security benefit ($)—00992.13969.75Flexibility benefit ($)—0183.870163.11Profit ($)−1347.41273.29380.581181.61275.22Table6 presents the economic indices of the proposed strategy that makes significant profit to all grids. The proposed strategy’s profit is equal to the difference between the sum of flexibility and security benefits and the sum of NRES operation and network energy costs. In Table 6, the cost of energy for this network is negative in Cases II-V. This means that the network is selling energy to upstream network in most of operating hours, as shown in Figure 7. Therefore, it earns income from the upstream network by selling energy. In Table 6, the energy cost is generally negative for this network in Cases II-V, meaning the revenue in this condition. Moreover, due to flexibility, security, and energy management, the sum of the revenue is greater than the NRES operation cost for all networks in the proposed scheme (Case V). It is worth noting that, by considering the improvement of operation, flexibility, and voltage security indices of the network, the proposed scheme (Case V) makes a higher profit for the mentioned networks compared to Cases II-IV. This is in accordance with the advantage of the suggested scheme compared to the literature listed in Table 1 or Section 1.


## 4.1. Case Study

The proposed two-layer CEMS-based HFSO is applied on the 33-bus distribution network shown in Figure4 [58]. This network includes radial MG1 (15-bus), MG2 (13-bus), and MG3 (14-bus) located, respectively, on buses 22, 14, and 33. One should note that the line, station, and peak load data of these networks are presented in [28], and the load value for simulation hours is equal to the multiplication of the peak load value and daily load factor curve plotted in Figure 5 [61]. In addition, there are several NRES diesel generators, and the RES includes photovoltaic and wind systems, battery energy storage in the proposed SDN and MGs, their location and capacity, and other data described in [28]. The daily forecasted power percentage of RESs is based on Figure 5 [61]. In this paper, it is assumed that the power factor of MGs load is 0.9, and all loads can participate in the DRP by considering ξ = 0.3. Finally, the daily energy price is presented in [20], and the range of voltage, flexibility and security incentive prices, and minimum allowed value of WSI are assumed, respectively, to be [0.9 p.u, 1.1 p.u] [62], 10 and 10 $/MWh, and 0.8.Figure 4
The scheme of the 33-bus distribution network including MMG [58].Figure 5
Forecasted power percentage of RES and network load [61].In the SDN, the slack bus (bus 1) is assumed to have a voltage magnitude of 1 p.u., and the voltage angle is zero. Each MG also consists of sources, ESSs, and DRP, which are expected to be managed so that each MG in the bus connected to the SDN can make the voltage magnitude of the connecting bus equal to 1 p.u. Thus, buses 14, 22, and 33, on which MG2, MG1, and MG3 are placed, respectively, are considered as the slack buses for the mentioned MGs, with a voltage magnitude of 1 p.u. However, the voltage angle of these buses is equal to that obtained from the calculations in the proposed scheme.


## 4.2. Results

The proposed linear form of the two-layer CEMS-based HFSO is simulated in GAMS23.5.2, and thus it is solved by the CPLEX method [63] by considering 180 linearization segments for circular constraints. The RWM generates 1000 scenario samples, and then the Kantorovich method reduces this number to 20 scenario samples with high probability.


### 4.2.1. Comparison of the convergence results of the proposed scheme for linear and nonlinear models

The MSD convergence for the proposed scheme in the MINLP and MILP models is investigated in Table4. In the MINLP model, BARON, BONMIN, DISOPT, KNITRO, and OQNLP solvers are used [63]. According to Table 4, among these solvers, DISOPT and OQNLP failed to find the optimal solution to the proposed scheme. Other algorithms also have different optimum points, meaning that the values of their objective functions are not the same. Therefore, MINLP solvers cannot achieve the unique optimal solution, which reduces the reliability. However, among the mentioned algorithms, it can be seen that BONMIN has a more desirable situation than other solvers because it has the shortest computational time (3027.11 s) and fewest convergence iterations for the MSD (12 iterations), and calculates the minimum objective function ($2621.73). Then, BONMIN, CBC, and CPLEX solvers are used to solve the proposed MILP model [42]. According to Table 4, it can be seen that these algorithms obtain a unique optimal solution; i.e., they reach a constant value for the objective function in the three algorithms, in the least number of MSD convergence iterations (8 iterations). Thus, these results have higher reliability. They also have a lower computational time than solvers in the MINLP model, where CPLEX has the shortest possible time (102.68 s). Therefore, the MILP model with CPLEX solver is suitable for solving the proposed problem in MSD format (Section 3.3), which has low computational time and low computational error based on Section 3.1.Table 4
Convergence results of the proposed problem based on MSD obtained by different solvers.
ModelSolverObjective function ($) in problem ofConvergence iteration forCalculation timesModel statusFirst layerSecond layerFirst layer problemSecond layer problemMSDMINLPBARON2702.651589.349284143273.43Locally optimalBONMIN2621.731537.617568123027.11Locally optimalDISOPT——————InfeasibleKNITRO2768.281611.5910396173347.38Locally optimalOQNLP——————InfeasibleMILPBONMIN2182.8861275.2274698129.57OptimalCBC2182.8861275.2264608113.83OptimalCPLEX2182.8861275.2255518102.68Optimal


### 4.2.2. Optimal Scheduling of MG Devices

The results of the HFSO in each MG based on the first layer model of CEMS are presented in this section to determine the optimal scheduling of MG sources and active loads. All these numerical results are plotted in Figure6. Accordingly, the RESs, i.e., photovoltaic and wind systems, inject their maximum active power into MGs due to their zero operation cost so that the power level of total RESs in the MG1 and MG2 is generally less than network load in simulation for all hours. However, it is close to the total active load at the periods of 6 : 00–13 : 00 and 4 : 00–13 : 00 for MG1 and MG2, respectively. RES active power is also close to/greater than the network load at the period of 4 : 00–6:00/6 : 00–16 : 00 for MG3. Moreover, the NRESs, i.e., diesel generators existing in MG1 and MG2, are generally switched off from 1 : 00 to 6 : 00, because the fuel price of NRESs is greater than the energy price in this period. But, in the other hours, it injects high active power into MGs due to the small fuel price considering energy price and large capacity in comparison with other flexible sources. Finally, the active loads, i.e., ESS and DRP, are charged/discharged in the low/high energy price hours or period or 1 : 00–7:00/18 : 00–22 : 00, to obtain minimum energy cost for MGs. However, they are discharged at some hours of the period of medium energy price or 8 : 00–17 : 00 and 23 : 00–00 : 00, to achieve the maximum security energy or benefit. These loads are charged at other times of this period to provide the required energy for their discharge operation. As the last point, the MG devices operate according to Figure 6 to achieve the minimum energy cost received from the upstream network or SDN and to obtain maximum security and flexibility for MG by satisfying the MG operation limits such as voltage, capacity of MG lines, and power factor limitations.Figure 6
Daily expected total active power curve of different devices in (a) MG1, (b) MG2, and (c) MG3.
(a)(b)(c)


### 4.2.3. Optimal Operation of the SDN

The daily expected active and reactive power of MGs absorbed/injected from/into the SDN according to the first layer CEMS-based HFSO are shown in Figure7. Based on Figure 7(a), MG1 receives the active power from the SDN at the period 1 : 00–7:00 due to low energy price in this interval, but it injects 3 MW active power into the SDN in other hours, arising from high injection power of MG local sources, i.e., NRES, RES, ESS, and DRP, according to Figure 6(a). However, the power generation capacity of these sources in MG2 and MG3 is lower than that in MG1. Hence, these MGs cannot inject the active power into the SDN at more simulation times as seen in Figure 7(a). Moreover, the ESS and NRES can provide the reactive load of MG1. Thus, MG1 does not receive/deliver reactive power from/to the SDN to obtain unit power factor and suitable regulation of bus voltages based on Figure 7(b). However, this condition cannot happen in MG2 and MG3, due to the low reactive power generation capacity of ESS and NRES in these MGs. Therefore, they absorb the reactive power from the SDN regarding Figure 7(b).Figure 7
Daily expected power curve of different MGs in the SDN: (a) active power; (b) reactive power.
(a)(b)In addition, the optimal operation of the SDN devices based on HFSO matching the second layer CEMS is illustrated in Figure8, where the same result as Figure 6 can be observed to obtain the minimum cost of energy received from the upstream grid or the sub-transmission network and high flexibility and security to the SDN subject to SDN operation limits, (8)–(11). Figure 9 shows the daily scheduling active and reactive power of the SDN station regarding the proposed strategy.Figure 8
Daily expected total active power curve of different devices in the SDN.Figure 9
Daily expected active/reactive power curve of the SDN station.As can be seen in Figure9, the SDN absorbs the active power from the sub-transmission network from 1 : 00 to 7 : 00 due to low energy prices at these hours. However, it has been able to inject active power of 1 to 2.5 MW into the upstream network at other hours as MGs and SDN devices deliver high active power to the SDN in the period of 8 : 00–00 : 00 according to Figures 7(a) and 8. Moreover, since the reactive power generation capacity of ESS and NRES in the SDN is low, this network cannot obtain a unity power factor and generally receives reactive power from the upstream grid based on Figure 9.


### 4.2.4. Capabilities of HFSO-based two-layer CEMS

This section considers the following cases to investigate the capability of the proposed scheme:(i)
Case I: power flow analysis.(ii)
Case II: proposed scheme considering SIP and FIP = 0.(iii)
Case III: proposed scheme considering SIP = 0 and FIP = 10 $/MWh.(iv)
Case IV: proposed scheme considering SIP = 10 $/MWh and FIP = 0.(v)
Case V: proposed scheme considering SIP and FIP = 10 $/MWh.Tables5 and 6 express, respectively, the technical and economic capabilities of the proposed strategy. Based on Table 5, there is high energy loss in power flow studies, in which the presence of sources, ESSs, and DRP is not taken into account in MGs and SDN. According to Figures 5 and 7, the total energy demand of all loads in MG1–3 and SDN is 126.4364 MWh, 121.8664 MWh, 83.7832 MWh, and 56.6 MWh, respectively. Therefore, energy loss in the networks in Case I is roughly 5.35%, 4.83%, 2.15%, and 4.85% of their consumption energy. Nevertheless, by proper management of sources, ESSs, and DRP in other study cases (II–V), energy loss reduces compared to Case I. In the proposed scheme (Case V), energy loss in MG1–3 and SDN has been reduced by about 32.9% ((6.76–4.533)/6.76), 31%, 36%, and 32% compared to Case I. Such conditions are also true for the maximum voltage drop of the network. The maximum voltage drop in Case I for the mentioned MGs is higher than the permissible range of 0.1 p.u. (1–0.9). However, its value in the proposed scheme (Case V) has been enhanced by about 22.5%, 63.2%, 80.9%, and 49.5% compared to Case I for MG1–3 and SDN. Concerning voltage security, the mentioned networks in Case I have the sum value of WSI close to their minimum value of 19.2 (24 × 0.8). Yet the proposed scheme succeeds to improve this index by about 14.66%, 5.18%, 16.11%, and 19.37% for MG1–3 and SDN. It is worth noting that Case V in terms of different indices enjoys a lower percentage of improvement compared to Cases II-IV. For instance, in terms of energy loss, Case II provides more suitable results than Case V because, in this case, minimization of the expected cost of energy purchase from the upstream network and the expected operation cost of sources is only considered. To reduce energy cost, energy consumption and energy loss need to be small in the network as per (1). Hence, in Case II, less energy loss resulted compared to Case V in which the economic status of the operation, flexibility, and network security indices is taken into account simultaneously. Concerning voltage security and maximum voltage drop, Case IV provides more desirable results than Case V because objective function (1) in this case study only considers the improvement of operation and network security statuses. Consequently, network operators encourage sources, ESSs, and DRP to inject more energy into MGs and SDN. As a result, by supplying more security energy (SE) than other cases, this case succeeds to obtain a more suitable WSI and voltage drop than Case V and even other case studies. However, it should be noted that these circumstances are reached for an increase in the maximum overvoltage of the networks, but their values are smaller than the permissible limit of 0.1 p.u. (1.1−1). Concerning flexibility, Case III is superior to Case V because it takes into account the improvement of operation and flexibility indices in the objective function (1). By supplying more flexibility energy (FE), Case III provides more flexibility than other case studies. Furthermore, Cases II-V obtain suitable power factors for distribution substations of different networks so that the distribution substation in MG1 has always a power factor of unity. This parameter in Case I is smaller than the allowable value of 0.9 p.u. In the following, as seen in Table 5, MG1, MG2, and SDN have achieved suitable WSI by obtaining higher security energy (SE). However, MG3 needs a low SE to reach the suitable WSI. This fact is due to the high voltage drop in MG1, MG2, and SDN compared to MG3. Hence, MG3 elements inject low SE into MG3 compared to other networks. Finally, it should be noted that the active power generation capacity of flexible sources, i.e., NRES, ESS, and DRP, in MG1 is higher than that of other grids based on Section 4.1. Therefore, the FE in this network is greater than that of MG1, MG2, and SDN.Table 5
The values of different technical indices in MGs and SDN.
MG1IndexCase ICase IICase IIICase IVCase VOperationDaily energy loss (MWh)6.764.434.644.614.533Maximum voltage drop/overvoltage (p.u.)0.102/00.071/0.0450.080/0.0430.070/0.050.079/0.046Minimum power factor of station0.861111SecurityPoor voltage bus1413131313Sum of WSI19.3322.02121.98522.21222.164SE (MWh)—155.26154.11168.47164.013FlexibilityFE (MWh)—54.28165.84553.49663.976MG2IndexCase ICase IICase IIICase IVCase VOperationDaily energy loss (MWh)5.893.9744.1124.1064.067Maximum voltage drop/overvoltage (p.u.)0.143/00.051/0.0360.055/0.0350.050/0.0400.054/0.038Minimum power factor of station0.880.900.900.900.90SecurityPoor voltage bus1212121212Sum of WSI19.2720.09120.03520.31420.268SE (MWh)—117.562116.792126.113123.849FlexibilityFE (MWh)—25.19233.45824.27730.953MG3IndexCase ICase IICase IIICase IVCase VOperationDaily energy loss (MWh)1.8061.0581.2371.2181.156Maximum voltage drop/overvoltage (p.u.)0.136/00.022/0.0150.028/0.0130.021/0.0170.026/0.016Minimum power factor of station0.870.900.900.900.90SecurityPoor voltage bus1414141414Sum of WSI19.4822.17522.05722.84622.619SE (MWh)—67.31566.10478.97475.94FlexibilityFE (MWh)—16.3321.7215.5119.60SDNIndexCase ICase IICase IIICase IVCase VOperationDaily energy loss (MWh)2.751.731.9741.9521.87Maximum voltage drop/overvoltage (p.u.)0.097/00.046/0.0390.051/0.0380.045/0.0420.049/0.041Minimum power factor of station0.830.930.920.920.92SecurityPoor voltage bus1831313131Sum of WSI19.6722.9822.8323.7123.48SE (MWh)—88.37186.92299.21396.975FlexibilityFE (MWh)—14.54518.38713.79416.311Table 6
The values of different economic indices in MGs and SDN.
CaseIIIIIIIVVMG1Energy cost ($)2836.11−1145.37−1094.22−1109.23−1090.732NRES operation cost ($)02638.122678.432704.522749.906Security benefit ($)—001684.71640.13Flexibility benefit ($)—0658.450639.76Profit ($)−2836.11−1492.75−925.7689.41620.716MG2Energy cost ($)1022.43−213.18−200.33−204.27−198.324NRES operation cost ($)0982.341001.121015.461032.385Security benefit ($)—001261.131238.49Flexibility benefit ($)—0334.580309.53Profit ($)−1022.43−769.16−466.21449.94713.96MG3Energy cost ($)473.6399.45105.08103.96106.79NRES operation cost ($)00000Security benefit ($)—00789.74759.40Flexibility benefit ($)—0217.20195.60Profit ($)−473.63−99.45112.12685.78848.21SDNEnergy cost ($)1347.41−1085.29−1030.35−1045.74−1022.362NRES operation cost ($)0812.31833.64856.27880Security benefit ($)—00992.13969.75Flexibility benefit ($)—0183.870163.11Profit ($)−1347.41273.29380.581181.61275.22Table6 presents the economic indices of the proposed strategy that makes significant profit to all grids. The proposed strategy’s profit is equal to the difference between the sum of flexibility and security benefits and the sum of NRES operation and network energy costs. In Table 6, the cost of energy for this network is negative in Cases II-V. This means that the network is selling energy to upstream network in most of operating hours, as shown in Figure 7. Therefore, it earns income from the upstream network by selling energy. In Table 6, the energy cost is generally negative for this network in Cases II-V, meaning the revenue in this condition. Moreover, due to flexibility, security, and energy management, the sum of the revenue is greater than the NRES operation cost for all networks in the proposed scheme (Case V). It is worth noting that, by considering the improvement of operation, flexibility, and voltage security indices of the network, the proposed scheme (Case V) makes a higher profit for the mentioned networks compared to Cases II-IV. This is in accordance with the advantage of the suggested scheme compared to the literature listed in Table 1 or Section 1.


## 4.2.1. Comparison of the convergence results of the proposed scheme for linear and nonlinear models

The MSD convergence for the proposed scheme in the MINLP and MILP models is investigated in Table4. In the MINLP model, BARON, BONMIN, DISOPT, KNITRO, and OQNLP solvers are used [63]. According to Table 4, among these solvers, DISOPT and OQNLP failed to find the optimal solution to the proposed scheme. Other algorithms also have different optimum points, meaning that the values of their objective functions are not the same. Therefore, MINLP solvers cannot achieve the unique optimal solution, which reduces the reliability. However, among the mentioned algorithms, it can be seen that BONMIN has a more desirable situation than other solvers because it has the shortest computational time (3027.11 s) and fewest convergence iterations for the MSD (12 iterations), and calculates the minimum objective function ($2621.73). Then, BONMIN, CBC, and CPLEX solvers are used to solve the proposed MILP model [42]. According to Table 4, it can be seen that these algorithms obtain a unique optimal solution; i.e., they reach a constant value for the objective function in the three algorithms, in the least number of MSD convergence iterations (8 iterations). Thus, these results have higher reliability. They also have a lower computational time than solvers in the MINLP model, where CPLEX has the shortest possible time (102.68 s). Therefore, the MILP model with CPLEX solver is suitable for solving the proposed problem in MSD format (Section 3.3), which has low computational time and low computational error based on Section 3.1.Table 4
Convergence results of the proposed problem based on MSD obtained by different solvers.
ModelSolverObjective function ($) in problem ofConvergence iteration forCalculation timesModel statusFirst layerSecond layerFirst layer problemSecond layer problemMSDMINLPBARON2702.651589.349284143273.43Locally optimalBONMIN2621.731537.617568123027.11Locally optimalDISOPT——————InfeasibleKNITRO2768.281611.5910396173347.38Locally optimalOQNLP——————InfeasibleMILPBONMIN2182.8861275.2274698129.57OptimalCBC2182.8861275.2264608113.83OptimalCPLEX2182.8861275.2255518102.68Optimal


## 4.2.2. Optimal Scheduling of MG Devices

The results of the HFSO in each MG based on the first layer model of CEMS are presented in this section to determine the optimal scheduling of MG sources and active loads. All these numerical results are plotted in Figure6. Accordingly, the RESs, i.e., photovoltaic and wind systems, inject their maximum active power into MGs due to their zero operation cost so that the power level of total RESs in the MG1 and MG2 is generally less than network load in simulation for all hours. However, it is close to the total active load at the periods of 6 : 00–13 : 00 and 4 : 00–13 : 00 for MG1 and MG2, respectively. RES active power is also close to/greater than the network load at the period of 4 : 00–6:00/6 : 00–16 : 00 for MG3. Moreover, the NRESs, i.e., diesel generators existing in MG1 and MG2, are generally switched off from 1 : 00 to 6 : 00, because the fuel price of NRESs is greater than the energy price in this period. But, in the other hours, it injects high active power into MGs due to the small fuel price considering energy price and large capacity in comparison with other flexible sources. Finally, the active loads, i.e., ESS and DRP, are charged/discharged in the low/high energy price hours or period or 1 : 00–7:00/18 : 00–22 : 00, to obtain minimum energy cost for MGs. However, they are discharged at some hours of the period of medium energy price or 8 : 00–17 : 00 and 23 : 00–00 : 00, to achieve the maximum security energy or benefit. These loads are charged at other times of this period to provide the required energy for their discharge operation. As the last point, the MG devices operate according to Figure 6 to achieve the minimum energy cost received from the upstream network or SDN and to obtain maximum security and flexibility for MG by satisfying the MG operation limits such as voltage, capacity of MG lines, and power factor limitations.Figure 6
Daily expected total active power curve of different devices in (a) MG1, (b) MG2, and (c) MG3.
(a)(b)(c)


## 4.2.3. Optimal Operation of the SDN

The daily expected active and reactive power of MGs absorbed/injected from/into the SDN according to the first layer CEMS-based HFSO are shown in Figure7. Based on Figure 7(a), MG1 receives the active power from the SDN at the period 1 : 00–7:00 due to low energy price in this interval, but it injects 3 MW active power into the SDN in other hours, arising from high injection power of MG local sources, i.e., NRES, RES, ESS, and DRP, according to Figure 6(a). However, the power generation capacity of these sources in MG2 and MG3 is lower than that in MG1. Hence, these MGs cannot inject the active power into the SDN at more simulation times as seen in Figure 7(a). Moreover, the ESS and NRES can provide the reactive load of MG1. Thus, MG1 does not receive/deliver reactive power from/to the SDN to obtain unit power factor and suitable regulation of bus voltages based on Figure 7(b). However, this condition cannot happen in MG2 and MG3, due to the low reactive power generation capacity of ESS and NRES in these MGs. Therefore, they absorb the reactive power from the SDN regarding Figure 7(b).Figure 7
Daily expected power curve of different MGs in the SDN: (a) active power; (b) reactive power.
(a)(b)In addition, the optimal operation of the SDN devices based on HFSO matching the second layer CEMS is illustrated in Figure8, where the same result as Figure 6 can be observed to obtain the minimum cost of energy received from the upstream grid or the sub-transmission network and high flexibility and security to the SDN subject to SDN operation limits, (8)–(11). Figure 9 shows the daily scheduling active and reactive power of the SDN station regarding the proposed strategy.Figure 8
Daily expected total active power curve of different devices in the SDN.Figure 9
Daily expected active/reactive power curve of the SDN station.As can be seen in Figure9, the SDN absorbs the active power from the sub-transmission network from 1 : 00 to 7 : 00 due to low energy prices at these hours. However, it has been able to inject active power of 1 to 2.5 MW into the upstream network at other hours as MGs and SDN devices deliver high active power to the SDN in the period of 8 : 00–00 : 00 according to Figures 7(a) and 8. Moreover, since the reactive power generation capacity of ESS and NRES in the SDN is low, this network cannot obtain a unity power factor and generally receives reactive power from the upstream grid based on Figure 9.


## 4.2.4. Capabilities of HFSO-based two-layer CEMS

This section considers the following cases to investigate the capability of the proposed scheme:(i)
Case I: power flow analysis.(ii)
Case II: proposed scheme considering SIP and FIP = 0.(iii)
Case III: proposed scheme considering SIP = 0 and FIP = 10 $/MWh.(iv)
Case IV: proposed scheme considering SIP = 10 $/MWh and FIP = 0.(v)
Case V: proposed scheme considering SIP and FIP = 10 $/MWh.Tables5 and 6 express, respectively, the technical and economic capabilities of the proposed strategy. Based on Table 5, there is high energy loss in power flow studies, in which the presence of sources, ESSs, and DRP is not taken into account in MGs and SDN. According to Figures 5 and 7, the total energy demand of all loads in MG1–3 and SDN is 126.4364 MWh, 121.8664 MWh, 83.7832 MWh, and 56.6 MWh, respectively. Therefore, energy loss in the networks in Case I is roughly 5.35%, 4.83%, 2.15%, and 4.85% of their consumption energy. Nevertheless, by proper management of sources, ESSs, and DRP in other study cases (II–V), energy loss reduces compared to Case I. In the proposed scheme (Case V), energy loss in MG1–3 and SDN has been reduced by about 32.9% ((6.76–4.533)/6.76), 31%, 36%, and 32% compared to Case I. Such conditions are also true for the maximum voltage drop of the network. The maximum voltage drop in Case I for the mentioned MGs is higher than the permissible range of 0.1 p.u. (1–0.9). However, its value in the proposed scheme (Case V) has been enhanced by about 22.5%, 63.2%, 80.9%, and 49.5% compared to Case I for MG1–3 and SDN. Concerning voltage security, the mentioned networks in Case I have the sum value of WSI close to their minimum value of 19.2 (24 × 0.8). Yet the proposed scheme succeeds to improve this index by about 14.66%, 5.18%, 16.11%, and 19.37% for MG1–3 and SDN. It is worth noting that Case V in terms of different indices enjoys a lower percentage of improvement compared to Cases II-IV. For instance, in terms of energy loss, Case II provides more suitable results than Case V because, in this case, minimization of the expected cost of energy purchase from the upstream network and the expected operation cost of sources is only considered. To reduce energy cost, energy consumption and energy loss need to be small in the network as per (1). Hence, in Case II, less energy loss resulted compared to Case V in which the economic status of the operation, flexibility, and network security indices is taken into account simultaneously. Concerning voltage security and maximum voltage drop, Case IV provides more desirable results than Case V because objective function (1) in this case study only considers the improvement of operation and network security statuses. Consequently, network operators encourage sources, ESSs, and DRP to inject more energy into MGs and SDN. As a result, by supplying more security energy (SE) than other cases, this case succeeds to obtain a more suitable WSI and voltage drop than Case V and even other case studies. However, it should be noted that these circumstances are reached for an increase in the maximum overvoltage of the networks, but their values are smaller than the permissible limit of 0.1 p.u. (1.1−1). Concerning flexibility, Case III is superior to Case V because it takes into account the improvement of operation and flexibility indices in the objective function (1). By supplying more flexibility energy (FE), Case III provides more flexibility than other case studies. Furthermore, Cases II-V obtain suitable power factors for distribution substations of different networks so that the distribution substation in MG1 has always a power factor of unity. This parameter in Case I is smaller than the allowable value of 0.9 p.u. In the following, as seen in Table 5, MG1, MG2, and SDN have achieved suitable WSI by obtaining higher security energy (SE). However, MG3 needs a low SE to reach the suitable WSI. This fact is due to the high voltage drop in MG1, MG2, and SDN compared to MG3. Hence, MG3 elements inject low SE into MG3 compared to other networks. Finally, it should be noted that the active power generation capacity of flexible sources, i.e., NRES, ESS, and DRP, in MG1 is higher than that of other grids based on Section 4.1. Therefore, the FE in this network is greater than that of MG1, MG2, and SDN.Table 5
The values of different technical indices in MGs and SDN.
MG1IndexCase ICase IICase IIICase IVCase VOperationDaily energy loss (MWh)6.764.434.644.614.533Maximum voltage drop/overvoltage (p.u.)0.102/00.071/0.0450.080/0.0430.070/0.050.079/0.046Minimum power factor of station0.861111SecurityPoor voltage bus1413131313Sum of WSI19.3322.02121.98522.21222.164SE (MWh)—155.26154.11168.47164.013FlexibilityFE (MWh)—54.28165.84553.49663.976MG2IndexCase ICase IICase IIICase IVCase VOperationDaily energy loss (MWh)5.893.9744.1124.1064.067Maximum voltage drop/overvoltage (p.u.)0.143/00.051/0.0360.055/0.0350.050/0.0400.054/0.038Minimum power factor of station0.880.900.900.900.90SecurityPoor voltage bus1212121212Sum of WSI19.2720.09120.03520.31420.268SE (MWh)—117.562116.792126.113123.849FlexibilityFE (MWh)—25.19233.45824.27730.953MG3IndexCase ICase IICase IIICase IVCase VOperationDaily energy loss (MWh)1.8061.0581.2371.2181.156Maximum voltage drop/overvoltage (p.u.)0.136/00.022/0.0150.028/0.0130.021/0.0170.026/0.016Minimum power factor of station0.870.900.900.900.90SecurityPoor voltage bus1414141414Sum of WSI19.4822.17522.05722.84622.619SE (MWh)—67.31566.10478.97475.94FlexibilityFE (MWh)—16.3321.7215.5119.60SDNIndexCase ICase IICase IIICase IVCase VOperationDaily energy loss (MWh)2.751.731.9741.9521.87Maximum voltage drop/overvoltage (p.u.)0.097/00.046/0.0390.051/0.0380.045/0.0420.049/0.041Minimum power factor of station0.830.930.920.920.92SecurityPoor voltage bus1831313131Sum of WSI19.6722.9822.8323.7123.48SE (MWh)—88.37186.92299.21396.975FlexibilityFE (MWh)—14.54518.38713.79416.311Table 6
The values of different economic indices in MGs and SDN.
CaseIIIIIIIVVMG1Energy cost ($)2836.11−1145.37−1094.22−1109.23−1090.732NRES operation cost ($)02638.122678.432704.522749.906Security benefit ($)—001684.71640.13Flexibility benefit ($)—0658.450639.76Profit ($)−2836.11−1492.75−925.7689.41620.716MG2Energy cost ($)1022.43−213.18−200.33−204.27−198.324NRES operation cost ($)0982.341001.121015.461032.385Security benefit ($)—001261.131238.49Flexibility benefit ($)—0334.580309.53Profit ($)−1022.43−769.16−466.21449.94713.96MG3Energy cost ($)473.6399.45105.08103.96106.79NRES operation cost ($)00000Security benefit ($)—00789.74759.40Flexibility benefit ($)—0217.20195.60Profit ($)−473.63−99.45112.12685.78848.21SDNEnergy cost ($)1347.41−1085.29−1030.35−1045.74−1022.362NRES operation cost ($)0812.31833.64856.27880Security benefit ($)—00992.13969.75Flexibility benefit ($)—0183.870163.11Profit ($)−1347.41273.29380.581181.61275.22Table6 presents the economic indices of the proposed strategy that makes significant profit to all grids. The proposed strategy’s profit is equal to the difference between the sum of flexibility and security benefits and the sum of NRES operation and network energy costs. In Table 6, the cost of energy for this network is negative in Cases II-V. This means that the network is selling energy to upstream network in most of operating hours, as shown in Figure 7. Therefore, it earns income from the upstream network by selling energy. In Table 6, the energy cost is generally negative for this network in Cases II-V, meaning the revenue in this condition. Moreover, due to flexibility, security, and energy management, the sum of the revenue is greater than the NRES operation cost for all networks in the proposed scheme (Case V). It is worth noting that, by considering the improvement of operation, flexibility, and voltage security indices of the network, the proposed scheme (Case V) makes a higher profit for the mentioned networks compared to Cases II-IV. This is in accordance with the advantage of the suggested scheme compared to the literature listed in Table 1 or Section 1.


## 5. Conclusion

In this paper, the HFSO method based on two-layer CEMS has been introduced to obtain optimal flexibility and security in the distribution network in the presence of multi-microgrids. The first layer of CEMS is based on the coordination between MG operators and DG, ESSs, and DRP operators. It corresponds to the optimal operation of MGs, including an objective function that minimizes the difference between the sum of the cost of energy received by the SDN and the operating cost of NRESs, and the sum of flexibility and security benefits. Moreover, it is subject to AC power flow equations, system operation limits, network flexibility and security model, and renewable and flexible source constraints. The second layer of CEMS is based on coordination between operators of the SDN and MGs, and it includes the operation model of the SDN in the presence of the MMG. Using the same model as in the first layer for the second layer of CEMS achieves the HFSO in the SDN. It should be noted that the proposed strategy includes a MINLP method converted to a MILP model based on the Taylor series and polygon approaches to access global optimal solutions with low calculation time and error. Then, according to a hybrid RWM-Kantorovich method, the SBSP models the uncertainties of load, energy price, and maximum active power of RESs. Finally, based on the obtained numerical results, it was observed that the solvers of the MILP model for the proposed scheme compared with the MINLP model of the problem can obtain the unique optimal solution in the shortest possible time, while solvers of MINLP reach a variety of solutions. Furthermore, compared to power flow studies, the proposed scheme has been able to improve energy loss, voltage profile, and voltage security by more than 30%, 22%, and 5%, respectively, by achieving optimal power scheduling for sources, ESSs, and DRP. Additionally, considering an economic model for flexibility, operation, and security indices, the proposed scheme provides a suitable approach by achieving suitable profit for MGs and SDN to encourage these networks to improve flexibility, operation, and security indices.Operation in this article was for normal conditions. However, as the aforementioned energy management system can be implemented if the necessary telecommunication platform is established, cyberattacks are possible in this situation. Therefore, strengthening the network against cyberattacks is suggested as future work. In addition, the effect of cyberattacks on the proposed scheme is considered for future work.



---

*Source: 3378538-2022-10-07.xml*