# IGDT-Based Robust Operation of Thermal and Electricity Energy-Based Microgrid with Distributed Sources, Storages, and Responsive Loads

Authors: Mehdi Sarlak; Abouzar Samimi; Mehdi Nikzad; Amir Hossein Salemi
Journal: International Transactions on Electrical Energy Systems
 (2022)

---

## Abstract

In this paper, the optimal operation of microgrids (MGs) with thermal blocks, distributed generations (DGs), storage systems, and responsive loads is presented to achieve optimal scheduling of active, reactive, and thermal power of the mentioned elements in the day-ahead (DA) reactive power and energy market environment. The thermal block has a combined heat and power (CHP) system, a boiler, and thermally responsive loads. This scheme minimizes the difference between the total operating costs of the MG and power sources and the total revenue gained from the sale of energy and reactive power of the mentioned elements in the markets located in the MG. It is constrained by the AC power flow equations, network operation constraints, and the operating model of these elements. Furthermore, this scheme is subject to the uncertainties of energy price, load, and renewable power. In this paper, to access the optimal resistant solution against the maximum prediction error associated with the mentioned uncertainties, a robust model based on information gap decision theory (IGDT) is used. Finally, by implementing the proposed scheme on a 119-bus radial MG, the obtained numerical results confirm the ability of the scheme to simultaneously improve the economic and operational situation of the MG. The proposed scheme succeeded in improving energy cost, energy loss, voltage drop, and power factor of the distribution substation by roughly 101%, 44%, 41%, and 16% compared to power flow studies, even in the worst-case scenario caused by uncertainties.


---

## Body

## 1. Introduction

### 1.1. Motivation

In order to reduce environmental pollution due to the uncontrolled consumption of fossil fuels, the use of environmentally friendly resources, storage devices, and energy management programs at the energy consumption site has been considered [1–4]. In addition, the presence of these elements at consumption points improves the technical and economic indices of the network [5]. For example, renewable energy sources (RESs) provide clean and low-cost energy in the network due to their low level of pollution emissions and low operating costs [6]. Energy storage systems (ESSs) and demand response programs (DRPs) with energy management are able to peak-shave network load profiles, resulting in improved operation indices such as reduced energy losses and voltage drops, and reducing the operating cost of the network [7]. Therefore, it is predicted that the number of these elements on the consumption side will increase. The utilization of aggregating units such as microgrids (MGs) for power sources, storage devices, and responsive loads was considered to establish easier and more desirable control and management of these elements [8]. In this scheme, the operator of the mentioned elements is in bilateral coordination with the MG operator (MGO). Therefore, MGOs with superiority over downstream data such as energy demand and power sources capacity for energy generation and the operator of the mentioned elements with superiority over upstream data such as energy price can optimize the situation for the technical and economic indices of the MG along with achieving optimal operation scheduling of ESSs and responsive loads [8]. Note, however, that achieving the above capabilities requires the establishment of an appropriate power management system (PMS) in the MG and the provision of an optimal formulation of network operation.


### 1.2. Literature Review

Various studies have been carried out in the field of optimal operation of microgrids. Literature [9] provides a method to evaluate the echelon utilization batteries (EUBs) in MGs from an economic point of view, where the MG participates in demand-side management and the EUBs are combined with air-conditioner load. By investigating the optimal control of these loads, a new control approach is introduced. Reference [10] develops a smart DRP to be used in MGs with generation units, photovoltaics (PVs), and battery ESSs (BESS). By doing this, the day-ahead (DA) generation planning and operation of BESS are optimized. Moreover, the paper presents a model to predict the output power of PVs based on the available data. The uncertain output of RESs is addressed in [11], where the paper provides a solution to straightforwardly connect RESs to MGs. A power management system (PMS) to manage the demand side is using electric springs. In the following, the operation of grid-connected MGs is improved from a flexibility point of view. An electric spring is a DC/AC converter and can control its flowing active and reactive power simultaneously thanks to using an insulated-gate bipolar transistor (IGBT) in its structure. Moreover, with connection to the grid from its AC side, it uses a controllable load to play a role in regulating the grid voltage [11]. The optimization problem includes an objective function to minimize several parameters such as operating cost and voltage variations and maximize the flexibility of the system. The application of an energy management system (EMS) in ships is introduced in [12], where the fuel consumption of the ship and, as a result, pollutant emission is reduced. Instead of an AC system, the ship employs a DC MG connected to an ESS and generators that work based on RESs. Electric springs and electric vehicles (EVs) parking lots are integrated to conquer the challenges imposed by the intermittent nature of RESs and other uncertainties for constructing a secure modern MG [13]. This can be formulated as a hybrid stochastic/robust optimization (HSRO) problem, in which uncertainties such as the amount of demand, price of energy, the output power of RESs, and the outage of MG devices are modeled by adopting stochastic programming that works based on unscented transformation (UT). Moreover, the paper adopts a robust model so that the uncertainties of EV parking lots are modeled and the operational situation of the MG is enhanced. The UT method is one of the approximation-based methods that are suitable for modeling uncertainties due to its ability in nonlinear correlated/uncorrelated transitions [13], acceptable probability distribution function (PDF) estimation, and simple coding. UT includes the fewest number of scenarios, so that it includes 2n + 1 (2n) scenarios if the number of uncertainty parameters is n [13]. In an attempt to enhance several indices including operation, reliability, pollutant emission, and economics at the same time. Reference [14] applies an EMS to unbalanced MGs containing active and reactive power sources and active loads. The approach is mathematically expressed in the form of an optimization problem with four objectives subject to optimal power flow (OPF) equations, limitations of the unbalanced MGs, reliability bounds of the system, and some other parameters. Also, uncertainties are modeled by stochastic programming. Literature [15] focuses on the work in [14] but considers balanced MGs. The DA and real-time (RT) energy markets establish the basis for the operation of a grid-connected MG, where the MG contains many RESs and EVs [16]. A multilayer EMS is also employed to manage the MG. The method divided the MGs into two categories: single MGs and a separate dominant MG that is connected to the distribution system and supervises other MGs. The EMS consists of two layers, and the problem has two stages. The details of this can be found in [16]. A new type of optimal energy management, named two-stage stochastic p-robust is used for an MG with PV systems, wind turbines, and microturbine (MT) [17]. To this end, a combined DRP is adopted, where incentive-based and price-based demand responses are combined in an attempt to decrease the amount of peak consumption while guarantying enhanced reliability. Reference [18] proposes an adaptive robust optimization to introduce a new operation approach to include several uncertainties such as RESs, demand, price of energy, and arrival and departure time of EVs in the operation model of MGs. The problem has been established in two parts. One part of the problem addresses the unit commitment (UC) situation of distributed generations (DGs). Also, the other part deals with worst-case uncertainty parameters.Reference [19] presents a model for the unit commitment of renewable energy hubs (EHs) that contain various types of demands and energy generation units. The information gap decision theory (IGDT) is used for the day-ahead scheduling of Ehs while taking into account diverse uncertainty parameters associated with electricity, heat, cooling, PV, and wind power, together with energy price. Storage devices are modeled in detail. Moreover, the impacts of risk and huge cost variations on the operating cost of EHs and the planning of EHs’ elements are thoroughly analyzed. Reference [20] focuses on hybrid combined cooling, heat, and power (CCHP) systems. The paper investigates two methods of utilizing solar energy. These methods transform solar energy into thermal and electrical energy using solar thermal and photovoltaic collectors. Also, the non-dominated sorting genetic algorithm II (NSGA-II) is used to find solutions to the problem while taking into account economic, energy, and environmental aspects. In another study, a robust dispatch model is incorporated into integrated electricity and heat networks (IEHNs) considering demand response and energy price [21]. The heating network is first modeled linearly, and then a computationally tractable PBIDR model is proposed. Reference [22] presents a robust scheduling model for MGs that operate based on CHP. The model is constructed with the help of IGDT. Traditional generation units are used in the MG, although recent technologies such as PV systems, wind turbines, batteries, and similar elements are also adopted as distributed generation sources. Reference [23] discussed a type of energy management strategy that uses the risk-aware IGDT. The strategy targets CCHP MGs that contain battery charging terminals. The elements of the MG are scheduled according to the decision-maker policy. A collaborative transaction between a residential CCHP plant and a demand aggregator was discussed in [24]. The paper investigated the changing cost of energy in the CCHP and adopted it to suggest a new strategy of energy pricing, according to which the price of electrical energy is fixed and that of heating and cooling energy is changing. The changes rely on the proportion of electrical, heating, and cooling demand in the system. A modeling of uncertain parameters such as building load, weather conditions, and human behavior that impact the comfort of the customers were presented in [25]. The paper introduced an adaptive, robust economic dispatching that consists of two stages. Along with this model of dispatching, a robust thermal comfort management strategy was also suggested so that various variables are taken into account to address uncertain parameters. In another study, the authors [26] discuss an adaptive stochastic method and attempt to reach optimal operation decisions, in which the benefits of adopting specific features of energy types were considered. The paper used a volt/var control to provide an optimization of active power flow and reactive power flow so that voltage security was ensured. It also presented modeling of battery degradation and thermal network to reach effective operation. Besides that, the authors employed a risk evaluation technique based on conditional value-at-risk so that reasonable solutions are obtained. Reference [27] introduced an optimal model of operation of an integrated energy system, in which the thermal inertia of a heating network and buildings were considered to promote the amount of wind power absorbed by the turbines. Reference [28] discussed another operation model of a microgrid that contains various types of energies so that generating units and flexible demand are properly dispatched. The proposed method attempted to model coupling limitations of electrical and heating networks, the dynamic features of the heating network, and power flow limits in a comprehensive way. Table 1 reports the majority of the studies presented within the last decade.Table 1
Taxonomy of recent research works.
Ref.Model of thermal energy in MGMarket modelUncertainty modelConsidering maximum forecasting error of uncertaintiesRiskmodel[9]✗EnergyStochastic✗✗[10]✗✗Stochastic✗✗[11]✗✗Stochastic✗✗[12]✗✗Stochastic✗✗[13]✗✗Stochastic-robust✗✗[14]✗✗Stochastic✗✗[15]✗✗Stochastic✗✗[16]✗EnergyStochastic✗✗[17]✗✗Stochastic✗✗[18]✗EnergyRobust✗✗[19]✗✗Robust✓✓[20]✗✗Stochastic✗✗[21]✓✗Robust✗✗[22]✗EnergyRobust✓✓[23]✗✗Robust✓✓[24]✓✗—✗✗[25]✗✗Robust✗✗[26]✗✗Stochastic-robust✗✓[27]✓✗—✗✗[28]✓✗—✗✗Proposed model✓Energy and reactive powerRobust✓✓


### 1.3. Research Gaps

According to the literature and Table1, there are the following major research gaps in the field of MG operation:(i)
In most studies, the electrical section model of the combined heat and power (CHP) has been considered in MG operation and planning problems. Moreover, in most research related to MG energy management, the electric consumer model has been proposed. Nonetheless, on the consumption side, the exploitation of electrical and thermal energy is evident. The CHP also generates both electrical and thermal energies at its output. Therefore, to better model MG operation, it is necessary to formulate the thermal model of sources and consumers. This has rarely been discussed such as in [19–24, 27, 28], but note that references [19, 20, 22, 23] present only a hybrid system of sources, ESS, and heat storage, while neglecting their connection to the energy or electrical network.(ii)
Note that by managing the scheduling of power sources, storage devices, and responsive loads, a suitable financial benefit can be obtained for them from the energy and ancillary services markets. However, in most pieces of research, the MG operation model has been considered, the aim of which is to improve the technical indices of the network, such as reducing power losses and voltage drop, along with achieving the lowest operating cost of MG and the mentioned elements. In a few references such as [9, 16, 18, 22], their participation in the energy market has been addressed. Also, these resources, with the ability to control reactive power, can play a role in ancillary services such as reactive power compensation. Thus, they can also benefit from the ancillary services market, which has rarely been discussed in the literature.(iii)
In most research such as [9–12, 14–17, 20], stochastic modeling has been considered for uncertainties. In this method, it is necessary to identify the probability distribution function (PDF) of these parameters. Since the accurate determination of PDF requires data collection over a long-term study period such as one year, these calculations will be very time-demanding. In addition, to achieve a reliable optimal solution in this method, a significant number of scenarios are needed. Hence, it is expected that the stochastic problem model for the proposed scheme has a high computational time. However, it should be noted that in operation problems, the low computational time is of great importance. As the execution step in this type of problem is small and it is assumed to take up to 15 minutes in real-time operation, stochastic modeling of uncertainties may not meet the above conditions.(iv)
In some research such as [13, 18, 21], robust programming has been used to model uncertainties. Adaptive robust optimization (ARO) and bounded uncertainty-based robust optimization (BURO) are generally used in the robust scheduling of the MG operation problem. However, these methods provide a solution robust against the predicted error defined for uncertainties. In other words, they are not able to obtain a robust solution and the maximum possible prediction error simultaneously. Furthermore, to achieve a reliable solution in the presence of uncertainties, it is necessary to consider the risk model, which cannot be evaluated using the mentioned methods. To address these limitations, robust programming based on IGDT is appropriate, which has been considered in fewer studies such as [13, 22, 23] related to MG operation problems.


### 1.4. Solutions and Contributions

In this paper, to deal with the first and second research gaps, power management in an MG with a thermal block, ESSs, DRPs, and DGs is proposed. This scheme is based on the participation of these elements in the DA energy and reactive power markets. In this paper, it is assumed that DGs, i.e., CHP and boilers, can provide the amount of energy consumed by heat consumers who also participate in DRP. Hence, these elements are placed in a complex unit called a “thermal block,” which is located in the buses of the MG where the heat load is present. In the following, the proposed scheme is expressed in the form of an optimization problem. In its deterministic model, its objective function is to minimize the difference between the total operating costs of DGs and MGs and the total revenue of power sources, storage devices, and responsive loads in these markets. It is also bound by the AC optimal power flow (AC-OPF) equations, and the operating model of the thermal block, DGs, ESSs, and DRPs. This scheme has uncertainties in energy price, renewable power, and load consumption. Therefore, to address the third and fourth research gaps, they are modeled based on IGDT-based robust scheduling. Also, to consider the adverse effects of uncertainties on the proposed scheme and the risk model, the IGDT model in this paper is based on a risk-averse strategy. This method can obtain the optimal robust solution in the worst-case scenario resulting from the maximum prediction error of the uncertainty parameters. Finally, the contributions of the mentioned scheme can be summarized as follows:(i)
Presenting an MG operation model with thermal blocks, power sources, storage devices, and responsive loads to achieve optimal operation and economic status for the network and the elements.(ii)
Achieving financial benefits for power sources, storage devices, and responsive loads from the energy market and ancillary services such as reactive power to encourage owners of these elements in the optimal operation of the MG.(iii)
Achieving the optimal solution robust against the worst-case scenario resulting from the maximum prediction error of uncertainties of load, renewable power, and energy prices using IGDT modeling of the problem based on a risk-averse strategy.


### 1.5. Paper Organization

The paper is organized as follows: In Section2, the final model of the proposed scheme is stated. Then, in Section 3, uncertainty modeling is presented. The numerical results obtained from different case studies are evaluated in Section 4. Finally, the conclusions are summarized in Section 5.


## 1.1. Motivation

In order to reduce environmental pollution due to the uncontrolled consumption of fossil fuels, the use of environmentally friendly resources, storage devices, and energy management programs at the energy consumption site has been considered [1–4]. In addition, the presence of these elements at consumption points improves the technical and economic indices of the network [5]. For example, renewable energy sources (RESs) provide clean and low-cost energy in the network due to their low level of pollution emissions and low operating costs [6]. Energy storage systems (ESSs) and demand response programs (DRPs) with energy management are able to peak-shave network load profiles, resulting in improved operation indices such as reduced energy losses and voltage drops, and reducing the operating cost of the network [7]. Therefore, it is predicted that the number of these elements on the consumption side will increase. The utilization of aggregating units such as microgrids (MGs) for power sources, storage devices, and responsive loads was considered to establish easier and more desirable control and management of these elements [8]. In this scheme, the operator of the mentioned elements is in bilateral coordination with the MG operator (MGO). Therefore, MGOs with superiority over downstream data such as energy demand and power sources capacity for energy generation and the operator of the mentioned elements with superiority over upstream data such as energy price can optimize the situation for the technical and economic indices of the MG along with achieving optimal operation scheduling of ESSs and responsive loads [8]. Note, however, that achieving the above capabilities requires the establishment of an appropriate power management system (PMS) in the MG and the provision of an optimal formulation of network operation.


## 1.2. Literature Review

Various studies have been carried out in the field of optimal operation of microgrids. Literature [9] provides a method to evaluate the echelon utilization batteries (EUBs) in MGs from an economic point of view, where the MG participates in demand-side management and the EUBs are combined with air-conditioner load. By investigating the optimal control of these loads, a new control approach is introduced. Reference [10] develops a smart DRP to be used in MGs with generation units, photovoltaics (PVs), and battery ESSs (BESS). By doing this, the day-ahead (DA) generation planning and operation of BESS are optimized. Moreover, the paper presents a model to predict the output power of PVs based on the available data. The uncertain output of RESs is addressed in [11], where the paper provides a solution to straightforwardly connect RESs to MGs. A power management system (PMS) to manage the demand side is using electric springs. In the following, the operation of grid-connected MGs is improved from a flexibility point of view. An electric spring is a DC/AC converter and can control its flowing active and reactive power simultaneously thanks to using an insulated-gate bipolar transistor (IGBT) in its structure. Moreover, with connection to the grid from its AC side, it uses a controllable load to play a role in regulating the grid voltage [11]. The optimization problem includes an objective function to minimize several parameters such as operating cost and voltage variations and maximize the flexibility of the system. The application of an energy management system (EMS) in ships is introduced in [12], where the fuel consumption of the ship and, as a result, pollutant emission is reduced. Instead of an AC system, the ship employs a DC MG connected to an ESS and generators that work based on RESs. Electric springs and electric vehicles (EVs) parking lots are integrated to conquer the challenges imposed by the intermittent nature of RESs and other uncertainties for constructing a secure modern MG [13]. This can be formulated as a hybrid stochastic/robust optimization (HSRO) problem, in which uncertainties such as the amount of demand, price of energy, the output power of RESs, and the outage of MG devices are modeled by adopting stochastic programming that works based on unscented transformation (UT). Moreover, the paper adopts a robust model so that the uncertainties of EV parking lots are modeled and the operational situation of the MG is enhanced. The UT method is one of the approximation-based methods that are suitable for modeling uncertainties due to its ability in nonlinear correlated/uncorrelated transitions [13], acceptable probability distribution function (PDF) estimation, and simple coding. UT includes the fewest number of scenarios, so that it includes 2n + 1 (2n) scenarios if the number of uncertainty parameters is n [13]. In an attempt to enhance several indices including operation, reliability, pollutant emission, and economics at the same time. Reference [14] applies an EMS to unbalanced MGs containing active and reactive power sources and active loads. The approach is mathematically expressed in the form of an optimization problem with four objectives subject to optimal power flow (OPF) equations, limitations of the unbalanced MGs, reliability bounds of the system, and some other parameters. Also, uncertainties are modeled by stochastic programming. Literature [15] focuses on the work in [14] but considers balanced MGs. The DA and real-time (RT) energy markets establish the basis for the operation of a grid-connected MG, where the MG contains many RESs and EVs [16]. A multilayer EMS is also employed to manage the MG. The method divided the MGs into two categories: single MGs and a separate dominant MG that is connected to the distribution system and supervises other MGs. The EMS consists of two layers, and the problem has two stages. The details of this can be found in [16]. A new type of optimal energy management, named two-stage stochastic p-robust is used for an MG with PV systems, wind turbines, and microturbine (MT) [17]. To this end, a combined DRP is adopted, where incentive-based and price-based demand responses are combined in an attempt to decrease the amount of peak consumption while guarantying enhanced reliability. Reference [18] proposes an adaptive robust optimization to introduce a new operation approach to include several uncertainties such as RESs, demand, price of energy, and arrival and departure time of EVs in the operation model of MGs. The problem has been established in two parts. One part of the problem addresses the unit commitment (UC) situation of distributed generations (DGs). Also, the other part deals with worst-case uncertainty parameters.Reference [19] presents a model for the unit commitment of renewable energy hubs (EHs) that contain various types of demands and energy generation units. The information gap decision theory (IGDT) is used for the day-ahead scheduling of Ehs while taking into account diverse uncertainty parameters associated with electricity, heat, cooling, PV, and wind power, together with energy price. Storage devices are modeled in detail. Moreover, the impacts of risk and huge cost variations on the operating cost of EHs and the planning of EHs’ elements are thoroughly analyzed. Reference [20] focuses on hybrid combined cooling, heat, and power (CCHP) systems. The paper investigates two methods of utilizing solar energy. These methods transform solar energy into thermal and electrical energy using solar thermal and photovoltaic collectors. Also, the non-dominated sorting genetic algorithm II (NSGA-II) is used to find solutions to the problem while taking into account economic, energy, and environmental aspects. In another study, a robust dispatch model is incorporated into integrated electricity and heat networks (IEHNs) considering demand response and energy price [21]. The heating network is first modeled linearly, and then a computationally tractable PBIDR model is proposed. Reference [22] presents a robust scheduling model for MGs that operate based on CHP. The model is constructed with the help of IGDT. Traditional generation units are used in the MG, although recent technologies such as PV systems, wind turbines, batteries, and similar elements are also adopted as distributed generation sources. Reference [23] discussed a type of energy management strategy that uses the risk-aware IGDT. The strategy targets CCHP MGs that contain battery charging terminals. The elements of the MG are scheduled according to the decision-maker policy. A collaborative transaction between a residential CCHP plant and a demand aggregator was discussed in [24]. The paper investigated the changing cost of energy in the CCHP and adopted it to suggest a new strategy of energy pricing, according to which the price of electrical energy is fixed and that of heating and cooling energy is changing. The changes rely on the proportion of electrical, heating, and cooling demand in the system. A modeling of uncertain parameters such as building load, weather conditions, and human behavior that impact the comfort of the customers were presented in [25]. The paper introduced an adaptive, robust economic dispatching that consists of two stages. Along with this model of dispatching, a robust thermal comfort management strategy was also suggested so that various variables are taken into account to address uncertain parameters. In another study, the authors [26] discuss an adaptive stochastic method and attempt to reach optimal operation decisions, in which the benefits of adopting specific features of energy types were considered. The paper used a volt/var control to provide an optimization of active power flow and reactive power flow so that voltage security was ensured. It also presented modeling of battery degradation and thermal network to reach effective operation. Besides that, the authors employed a risk evaluation technique based on conditional value-at-risk so that reasonable solutions are obtained. Reference [27] introduced an optimal model of operation of an integrated energy system, in which the thermal inertia of a heating network and buildings were considered to promote the amount of wind power absorbed by the turbines. Reference [28] discussed another operation model of a microgrid that contains various types of energies so that generating units and flexible demand are properly dispatched. The proposed method attempted to model coupling limitations of electrical and heating networks, the dynamic features of the heating network, and power flow limits in a comprehensive way. Table 1 reports the majority of the studies presented within the last decade.Table 1
Taxonomy of recent research works.
Ref.Model of thermal energy in MGMarket modelUncertainty modelConsidering maximum forecasting error of uncertaintiesRiskmodel[9]✗EnergyStochastic✗✗[10]✗✗Stochastic✗✗[11]✗✗Stochastic✗✗[12]✗✗Stochastic✗✗[13]✗✗Stochastic-robust✗✗[14]✗✗Stochastic✗✗[15]✗✗Stochastic✗✗[16]✗EnergyStochastic✗✗[17]✗✗Stochastic✗✗[18]✗EnergyRobust✗✗[19]✗✗Robust✓✓[20]✗✗Stochastic✗✗[21]✓✗Robust✗✗[22]✗EnergyRobust✓✓[23]✗✗Robust✓✓[24]✓✗—✗✗[25]✗✗Robust✗✗[26]✗✗Stochastic-robust✗✓[27]✓✗—✗✗[28]✓✗—✗✗Proposed model✓Energy and reactive powerRobust✓✓


## 1.3. Research Gaps

According to the literature and Table1, there are the following major research gaps in the field of MG operation:(i)
In most studies, the electrical section model of the combined heat and power (CHP) has been considered in MG operation and planning problems. Moreover, in most research related to MG energy management, the electric consumer model has been proposed. Nonetheless, on the consumption side, the exploitation of electrical and thermal energy is evident. The CHP also generates both electrical and thermal energies at its output. Therefore, to better model MG operation, it is necessary to formulate the thermal model of sources and consumers. This has rarely been discussed such as in [19–24, 27, 28], but note that references [19, 20, 22, 23] present only a hybrid system of sources, ESS, and heat storage, while neglecting their connection to the energy or electrical network.(ii)
Note that by managing the scheduling of power sources, storage devices, and responsive loads, a suitable financial benefit can be obtained for them from the energy and ancillary services markets. However, in most pieces of research, the MG operation model has been considered, the aim of which is to improve the technical indices of the network, such as reducing power losses and voltage drop, along with achieving the lowest operating cost of MG and the mentioned elements. In a few references such as [9, 16, 18, 22], their participation in the energy market has been addressed. Also, these resources, with the ability to control reactive power, can play a role in ancillary services such as reactive power compensation. Thus, they can also benefit from the ancillary services market, which has rarely been discussed in the literature.(iii)
In most research such as [9–12, 14–17, 20], stochastic modeling has been considered for uncertainties. In this method, it is necessary to identify the probability distribution function (PDF) of these parameters. Since the accurate determination of PDF requires data collection over a long-term study period such as one year, these calculations will be very time-demanding. In addition, to achieve a reliable optimal solution in this method, a significant number of scenarios are needed. Hence, it is expected that the stochastic problem model for the proposed scheme has a high computational time. However, it should be noted that in operation problems, the low computational time is of great importance. As the execution step in this type of problem is small and it is assumed to take up to 15 minutes in real-time operation, stochastic modeling of uncertainties may not meet the above conditions.(iv)
In some research such as [13, 18, 21], robust programming has been used to model uncertainties. Adaptive robust optimization (ARO) and bounded uncertainty-based robust optimization (BURO) are generally used in the robust scheduling of the MG operation problem. However, these methods provide a solution robust against the predicted error defined for uncertainties. In other words, they are not able to obtain a robust solution and the maximum possible prediction error simultaneously. Furthermore, to achieve a reliable solution in the presence of uncertainties, it is necessary to consider the risk model, which cannot be evaluated using the mentioned methods. To address these limitations, robust programming based on IGDT is appropriate, which has been considered in fewer studies such as [13, 22, 23] related to MG operation problems.


## 1.4. Solutions and Contributions

In this paper, to deal with the first and second research gaps, power management in an MG with a thermal block, ESSs, DRPs, and DGs is proposed. This scheme is based on the participation of these elements in the DA energy and reactive power markets. In this paper, it is assumed that DGs, i.e., CHP and boilers, can provide the amount of energy consumed by heat consumers who also participate in DRP. Hence, these elements are placed in a complex unit called a “thermal block,” which is located in the buses of the MG where the heat load is present. In the following, the proposed scheme is expressed in the form of an optimization problem. In its deterministic model, its objective function is to minimize the difference between the total operating costs of DGs and MGs and the total revenue of power sources, storage devices, and responsive loads in these markets. It is also bound by the AC optimal power flow (AC-OPF) equations, and the operating model of the thermal block, DGs, ESSs, and DRPs. This scheme has uncertainties in energy price, renewable power, and load consumption. Therefore, to address the third and fourth research gaps, they are modeled based on IGDT-based robust scheduling. Also, to consider the adverse effects of uncertainties on the proposed scheme and the risk model, the IGDT model in this paper is based on a risk-averse strategy. This method can obtain the optimal robust solution in the worst-case scenario resulting from the maximum prediction error of the uncertainty parameters. Finally, the contributions of the mentioned scheme can be summarized as follows:(i)
Presenting an MG operation model with thermal blocks, power sources, storage devices, and responsive loads to achieve optimal operation and economic status for the network and the elements.(ii)
Achieving financial benefits for power sources, storage devices, and responsive loads from the energy market and ancillary services such as reactive power to encourage owners of these elements in the optimal operation of the MG.(iii)
Achieving the optimal solution robust against the worst-case scenario resulting from the maximum prediction error of uncertainties of load, renewable power, and energy prices using IGDT modeling of the problem based on a risk-averse strategy.


## 1.5. Paper Organization

The paper is organized as follows: In Section2, the final model of the proposed scheme is stated. Then, in Section 3, uncertainty modeling is presented. The numerical results obtained from different case studies are evaluated in Section 4. Finally, the conclusions are summarized in Section 5.


## 2. A Deterministic Model of the Proposed Scheme

In this section, the mathematical operation model of MG with power sources, storage devices, and responsive loads in accordance with the DA energy and reactive power markets, as shown in Figure1 is presented in the framework of deterministic optimization. This scheme is expressed as an optimization formulation. This formulation includes objective functions [29–33] and constraints [34–37]. This scheme is an attempt to minimize the difference between the operating costs of nonrenewable DGs (NRDGs) and MG with the revenue from the sale of energy/reactive power in the mentioned markets by the mentioned elements. It is also bound by the AC-OPF equations, the operating model of NRDGs, ESSs, DRPs, and thermal blocks. It is assumed that in the buses of MG where there is a heat consumer, the heat block is also located in this bus. It consists of two NRDGs of CHP and boiler types, and DRP is also considered for the thermal load. Note that, the optimization process based on the proposed scheme can be implemented on the network if there is smart grid technology in this network [38]. The smart grid needs telecommunications [39–42] and smart devices [43–46].Figure 1
Framework of the proposed MG operation.The deterministic model of the proposed scheme can be given as follows:(1)minCost=∑h∈OHλhPrefG+∑b∈BusαbMT+βbMTPb,hMT+χbMTPb,hMT2+αbCHP+βbCHPPb,hCHP+χbCHPPb,hCHP2+αbBO+βbBOHb,hBO+χbBOHb,hBO2−ρhPb,hMT+Pb,hCHP+Pb,hRES+Pb,hDR+Pb,hDCH−Pb,hCH+KQQb,hMT+Qb,hCHP+Qb,hRES+Qb,hESS−γhHb,hCHP+Hb,hBO+Hb,hDR.Subject to:(2)Pb,hG+Pb,hMT+Pb,hRES+Pb,hCHP+Pb,hDCH−Pb,hCH+Pb,hDR−∑k∈BusAb,kLPb,k,hL=Pb,hD∀b,h,(3)Qb,hG+Qb,hMT+Qb,hRES+Qb,hCHP+Qb,hESS−∑k∈BusAb,kLQb,k,hL=Qb,hD∀b,h,(4)Pb,k,hL=Gb,kVb,h2−Vb,hVk,hGb,kcosδb,h−δk,h+Bb,ksinδb,h−δk,h∀b,k,h,(5)Qb,k,hL=−Bb,kVb,h2+Vb,hVk,hGb,kcosδb,h−δk,h−Bb,ksinδb,h−δk,h∀b,k,h,(6)δb,h=0∀b=ref,h,(7)Pb,k,hL2+Qb,k,hL2≤Sb,kLu∀b,k,h,(8)Pb,hG2+Qb,hG2≤SbGu∀b=ref,h,(9)Pb,hGPb,hG2+Qb,hG2≥PFl∀b=ref,h,(10)Vl≤Vb,h≤Vu∀b,h,(11)PbMTl≤Pb,hMT≤PbMTu∀b,h,(12)QbMTl≤Qb,hMT≤QbMTu∀b,h,(13)Pb,hRES=Pb,hRESu∀b,h,(14)Pb,hRES2+Qb,hRES2≤SbRESu∀b,h,(15)0≤Pb,hCH≤PvCHuxb,h∀b,h,(16)0≤Pb,hDCH≤PbDCHu1−xb,h∀b,h,(17)Eb,h+1=Eb,h+ηbCHPb,hCH−1ηbDCHPb,hDCH∀b,h,(18)Eb,h=Ebini∀b,h=1,(19)Ebl≤Eb,h≤Ebu∀b,h,(20)Pb,hDCH−Pb,hCH2+Qb,hESS2≤SbESSu∀b,h,(21)Hb,hCHP+Hb,hBO+Hb,hDR=Hb,hD∀b,h,(22)Hb,hCHP=Pb,hCHP1−ηbT−ηbLηbT.ηbH∀b,h,(23)PbCHPl≤Pb,hCHP≤PbCHPu∀b,h,(24)QbCHPl≤Qb,hCHP≤QbCHPu∀b,h,(25)HbCHPl≤Hb,hCHP≤HbCHPu∀b,h,(26)HbBOl≤Hb,hBO≤HbBOu∀b,h,(27)−ξbPb,hDorHb,hD≤Pb,hDRorHb,hDR≤ξbPb,hDorHb,hD∀b,h,(28)∑h∈OHPb,hDRorHb,hDR=0∀b.


### 2.1. Objective Function

Equation (1) formulates the objective function of the proposed scheme, which is equal to minimizing the difference between the operating costs of MGs and NRDGs [47] and the revenue of power sources, storage devices, and responsive loads in the energy and reactive power markets. In the first part of (1), the operating cost of MG, which represents the cost of purchasing energy from the upstream network, is modeled [19]. In the second part of this equation, the first to third terms refer to the operating (fuel) costs of NRDGs of MT, CHP, and boiler types, respectively [47]. The fourth to sixth terms relate to the revenue gained by selling electrical energy by power sources, storage devices, and responsive loads in the DA electricity market, respectively; the sales revenue of reactive power by sources and storage devices in the DA electrical reactive market; and revenue from selling energy by thermal block (CHP, boiler, and thermal DRP) in the DA thermal energy market. It is noteworthy that in this section it is assumed that these markets are established as retail markets in the MG environment, and also MG itself can purchase from the wholesale market in the upstream network. Therefore, the energy purchase price of MG (λ) will be different from the energy sale price in the retail market (ρ and γ). Also, in this section, it is assumed that the price of reactive power is a multiplication of the price of electrical energy, i.e., KQ × ρ. In addition, based on [15, 16], reactive power sources generally inject/receive reactive power into/from the network in exchange for MG demand; hence, the absolute value of reactive power of these resources is adopted in (1).


### 2.2. MG Operation Constraints

These constraints are modeled as AC-OPF in equations (2)–(10) [48–51]. Equations (2)–(6) are proportional to AC-PF [48, 49], which represents the balance of active and reactive power in MG buses, the calculation of active and reactive power passing through distribution lines, and the voltage angle of the slack bus. Limitations of MG operation are also presented in equations (7)–(10) [50, 51]. In equations (7) and (8), the capacity limit of the distribution line/distribution substation is considered. Constraint (9) also refers to the limitation of the power factor of the distribution substation, according to which this substation must always have a power factor greater than PFl (generally has a value of 0.9) in all operating hours. The voltage magnitude limit of MG buses is also formulated in equations (10). If it operates in islanded mode, then MG is disconnected from the upstream network, so the variables PG and QG in all buses will be zero. The proposed scheme can be used for different types of distribution networks in accordance with the problem model given in equations (1)–(28). This is because the optimal power flow model, (2)–(10), can be applied to all distribution networks. Index ref should be defined as the buses with the distribution substation just to assure this point.


### 2.3. The Operation Model of DGs

Constraints (11)–(14) include the operation models of MT (11) and (12), and RESs (13) and (14) [7, 14]. Constraints (11) and (12) refer to the MT capability curve model, which indicates the limit of active power generation and controllable reactive power of this type of NRDG, respectively. The amount of active power generation by RESs is presented in equation (13) [7]. Constraint (14) also considers the limit of apparent power transferable by these sources. It should be noted that photovoltaic (PV)-type RES is connected to the grid via a DC/AC converter [52], which considering the presence of an IGBT bridge in this converter can control the active and reactive power of PV at the same time [53]. A wind system (WS)-type RES is also connected to the network via an AC/AC converter, and it also has an induction generator [54]. Therefore, this source is also able to control active and reactive power at the same time. Considering these cases, constraint (14) can be modeled for RESs.


### 2.4. The Operation Model of Electrical Energy Storage Devices

The formulation of this section is stated in constraints (15)–(20) [55], which represent respectively the limitation of charge and discharge rate, stored energy, and initial energy in the first hour of operation, the limit of storable energy, and the charger capacity limit of the storage. In constraints (15) and (16), the binary variable x is used to prevent simultaneous storage operations in charge and discharge modes. This paper also assumes that storage devices are either connected to the grid via a charger with an IGBT bridge (such as a battery) [53] or have a synchronous generator (such as a compressed-air energy storage unit, CAES) [56]. Therefore, the storage device will be able to control its active and reactive power simultaneously, and in these conditions, constraint (20) can be modeled for this element.


### 2.5. The Operation Model of the Thermal Block

The thermal block includes a CHP, boiler, and DRP-capable heat consumers. Its operation formulation is presented in constraints (21)–(28). In equation (21), the thermal power balance in this block is modeled. In equations (22)–(25), the CHP model is considered [5]. Thermal power generation by CHP is calculated based on equation (22). Constraints (23) and (24) formulate the CHP capability curve in the electrical section, which refers to the limitations of active and reactive power controlled by CHP, respectively. The limit of thermal power that can be produced by CHP is also considered in equation (25). There is no restriction on the type of CHP in this article because the formulations (22)–(25) can be used for different types of CHP [5]. There is a constraint similar to equation (25) for the boiler that is modeled in equation (26) [5]. Finally, the DRP operation model for the thermal (electrical) consumer is described in constraints (27) and (28) [7]. This type of DRP fits the incentive model so that consumers in this DRP reduce or increase their energy consumption according to the energy price signal. So it is named price-based DRP. In other words, in DRP, to minimize equation (1), consumers reduce (increase) their energy consumption when the energy price is high (low). Therefore, constraint (27) is proportional to the range of changes in consumer power in DRP, and based on this equation, they participate in DRP with a maximum participation rate of ξ. In addition, it is assumed that in this DRP, consumers can receive all the energy they need from the retail energy market up to the operation horizon. Therefore, to meet these conditions, constraint (28) is added to the formulation of the proposed scheme.


## 2.1. Objective Function

Equation (1) formulates the objective function of the proposed scheme, which is equal to minimizing the difference between the operating costs of MGs and NRDGs [47] and the revenue of power sources, storage devices, and responsive loads in the energy and reactive power markets. In the first part of (1), the operating cost of MG, which represents the cost of purchasing energy from the upstream network, is modeled [19]. In the second part of this equation, the first to third terms refer to the operating (fuel) costs of NRDGs of MT, CHP, and boiler types, respectively [47]. The fourth to sixth terms relate to the revenue gained by selling electrical energy by power sources, storage devices, and responsive loads in the DA electricity market, respectively; the sales revenue of reactive power by sources and storage devices in the DA electrical reactive market; and revenue from selling energy by thermal block (CHP, boiler, and thermal DRP) in the DA thermal energy market. It is noteworthy that in this section it is assumed that these markets are established as retail markets in the MG environment, and also MG itself can purchase from the wholesale market in the upstream network. Therefore, the energy purchase price of MG (λ) will be different from the energy sale price in the retail market (ρ and γ). Also, in this section, it is assumed that the price of reactive power is a multiplication of the price of electrical energy, i.e., KQ × ρ. In addition, based on [15, 16], reactive power sources generally inject/receive reactive power into/from the network in exchange for MG demand; hence, the absolute value of reactive power of these resources is adopted in (1).


## 2.2. MG Operation Constraints

These constraints are modeled as AC-OPF in equations (2)–(10) [48–51]. Equations (2)–(6) are proportional to AC-PF [48, 49], which represents the balance of active and reactive power in MG buses, the calculation of active and reactive power passing through distribution lines, and the voltage angle of the slack bus. Limitations of MG operation are also presented in equations (7)–(10) [50, 51]. In equations (7) and (8), the capacity limit of the distribution line/distribution substation is considered. Constraint (9) also refers to the limitation of the power factor of the distribution substation, according to which this substation must always have a power factor greater than PFl (generally has a value of 0.9) in all operating hours. The voltage magnitude limit of MG buses is also formulated in equations (10). If it operates in islanded mode, then MG is disconnected from the upstream network, so the variables PG and QG in all buses will be zero. The proposed scheme can be used for different types of distribution networks in accordance with the problem model given in equations (1)–(28). This is because the optimal power flow model, (2)–(10), can be applied to all distribution networks. Index ref should be defined as the buses with the distribution substation just to assure this point.


## 2.3. The Operation Model of DGs

Constraints (11)–(14) include the operation models of MT (11) and (12), and RESs (13) and (14) [7, 14]. Constraints (11) and (12) refer to the MT capability curve model, which indicates the limit of active power generation and controllable reactive power of this type of NRDG, respectively. The amount of active power generation by RESs is presented in equation (13) [7]. Constraint (14) also considers the limit of apparent power transferable by these sources. It should be noted that photovoltaic (PV)-type RES is connected to the grid via a DC/AC converter [52], which considering the presence of an IGBT bridge in this converter can control the active and reactive power of PV at the same time [53]. A wind system (WS)-type RES is also connected to the network via an AC/AC converter, and it also has an induction generator [54]. Therefore, this source is also able to control active and reactive power at the same time. Considering these cases, constraint (14) can be modeled for RESs.


## 2.4. The Operation Model of Electrical Energy Storage Devices

The formulation of this section is stated in constraints (15)–(20) [55], which represent respectively the limitation of charge and discharge rate, stored energy, and initial energy in the first hour of operation, the limit of storable energy, and the charger capacity limit of the storage. In constraints (15) and (16), the binary variable x is used to prevent simultaneous storage operations in charge and discharge modes. This paper also assumes that storage devices are either connected to the grid via a charger with an IGBT bridge (such as a battery) [53] or have a synchronous generator (such as a compressed-air energy storage unit, CAES) [56]. Therefore, the storage device will be able to control its active and reactive power simultaneously, and in these conditions, constraint (20) can be modeled for this element.


## 2.5. The Operation Model of the Thermal Block

The thermal block includes a CHP, boiler, and DRP-capable heat consumers. Its operation formulation is presented in constraints (21)–(28). In equation (21), the thermal power balance in this block is modeled. In equations (22)–(25), the CHP model is considered [5]. Thermal power generation by CHP is calculated based on equation (22). Constraints (23) and (24) formulate the CHP capability curve in the electrical section, which refers to the limitations of active and reactive power controlled by CHP, respectively. The limit of thermal power that can be produced by CHP is also considered in equation (25). There is no restriction on the type of CHP in this article because the formulations (22)–(25) can be used for different types of CHP [5]. There is a constraint similar to equation (25) for the boiler that is modeled in equation (26) [5]. Finally, the DRP operation model for the thermal (electrical) consumer is described in constraints (27) and (28) [7]. This type of DRP fits the incentive model so that consumers in this DRP reduce or increase their energy consumption according to the energy price signal. So it is named price-based DRP. In other words, in DRP, to minimize equation (1), consumers reduce (increase) their energy consumption when the energy price is high (low). Therefore, constraint (27) is proportional to the range of changes in consumer power in DRP, and based on this equation, they participate in DRP with a maximum participation rate of ξ. In addition, it is assumed that in this DRP, consumers can receive all the energy they need from the retail energy market up to the operation horizon. Therefore, to meet these conditions, constraint (28) is added to the formulation of the proposed scheme.


## 3. Uncertain Modeling

### 3.1. Uncertain Parameters

In the proposed scheme, equations (1)–(28), parameters such as energy price (λ, ρ, and γ), load (PD, QD, and HD), and maximum renewable active power (PRESu) are uncertain. Hence, stochastic [57, 58], probabilistic, or robust [59, 60] modeling is needed for the proposed design. Since the first two types of modeling require a significant number of scenarios and identification of the PDF for uncertainties [61], the computational time is high [9–12]. However, as the implementation step is low in operation problems (up to 15 minutes in real-time operation), the low computational time is of great importance in these problems [16]. Computational time depends on the solution space of the problem [62]. Following this, a robust model of the proposed scheme is presented [5, 17]. In this section, to achieve a robust optimal solution against the maximum prediction error caused by the mentioned uncertainties, a robust model based on IGDT [63] is employed. Also, to investigate the adverse effect of uncertainties in the worst-case scenario, IGDT in accordance with the risk-averse strategy is used for the proposed scheme [63]. In this model, the vector of uncertainty parameters (u¯) is presented as (29), which is the expected (predicted) values of uncertainties. However, the true value of these parameters is a variable in the IGDT model, which following these conditions, the vector of uncertainty variables (u˜) can be expressed as (30). In IGDT, uncertainty is not a parameter, but it is defined as a variable. Also, since this variable is to determine the worst-case scenario this variable is a decision variable in IGDT.(29)u¯=PDQDHDPRESuλργ,(30)u˜=P˜DQ˜DH˜DP˜RESuλ˜ρ˜γ˜.


### 3.2. IGDT-Based Uncertainty Management

An IGDT-based model is used in this study to address uncertainties of demand, price of energy, and output power of RESs. The approach is independent of PDFs. The problem modeling is described as follows:(31)f=minYFY,u˜,(32)LY,u˜≤0,(33)MY,u˜=0,(34)u˜∈Ξ,where, u˜ denotes input uncertainty variables, and Ξ represents the set of uncertainties. Y shows decision-making variables. The present paper used an IGDT based on an information gap model. Even though various models have been used for uncertainty parameters in the literature, the current research employs the envelope-bound model to provide the prior information regarding the uncertainties as given as follows [63]:(35)Ξ=Ξu¯,σ=u˜:u˜−u¯u¯≤σ,where, u˜ represents the uncertain variables matrix, and σ is the difference between the maximum deviation of uncertainty variables and its estimated or nominal value (u¯). The variable σ is known as the uncertainty level and its value is assumed constant for all uncertainties in this paper. Therefore, it has a scalar value. The deterministic model described by equations (31)–(34) is an optimization problem with an estimated or nominal value of uncertainty variables, i.e., equations (31)–(33) together with only u¯.This mathematical expression is called the based condition (BC), where its objective function is given byf=minYFY,u¯. The maximum radius of uncertainty (rc) is found as follows [63]:(36)rc=maxY,σσ.Constraints (32) and (33),(37)FY,u˜≤Δc,(38)Δc=fbY,u¯×1+π,(39)−σ≤u˜−u¯u¯≤σ,(40)0≤σ≤1,where Δc denotes the robustness value of the problem described in equations (31)–(34). The value objective function should be less than Δc. Note that Δc is an input parameter to the problem and is found using (38). In (38), fb is the value of the BC objective function and π∈0⋅1 shows the degree of robustness when fb deteriorates. Consequently, in equations (36)–(40), π is an input parameter and fb shows an input parameter found by equations (31)–(33) taking into account only u¯.In formulation (36)–(40), the uncertainty level, σ, and uncertainty u˜ are variables. So, this problem obtains the highest possible uncertainty level, and the values of uncertainty variables, (30), are determined. The values obtained for the uncertainty variables correspond to the worst-case scenario because, according to (37) and (38), with increasing the value of the uncertainty level, the solution space of the proposed problem decreases compared to the deterministic model; hence, the value of the objective function or cost in these conditions is more than the deterministic model. Thus, the scenario obtained from problem (36)–(40) represents the worst-case scenario.Note that, in the model presented by (36)–(40), the uncertainty level for all uncertainties is assumed to be constant. In other words, it is assumed that the uncertainty radius of all uncertainties is the same. However, the mentioned model has no limitation on considering different uncertainty radii for different uncertainties. Under such conditions, an index such as i is added to σ, which includes uncertainty variables of equation (30). Then, equation (36) is changed as rc=maxY,σ∑iσi, but generally, concerning the IGDT, the uncertainty level is considered fixed for all uncertainties. This can be found in [19].


### 3.3. Robust Formulation of the Proposed Scheme

As per subsection3.2, the IGDT-based robust model proposed for the proposed problem, described by equations (1)–(28), can be given as follows:(41)rc=maxσ.Subject to:Constraints (2)–(28) with u˜ instead of u¯,(42)∑h∈OHλ˜hPrefG+∑b∈BusαbMT+βbMTPb,hMT+χbMTPb,hMT2+αbCHP+βbCHPPb,hCHP+χbCHPPb,hCHP2+αbBO+βbBOHb,hBO+χbBOHb,hBO2−ρ˜hPb,hMT+Pb,hCHP+P˜b,hRES+Pb,hDR+Pb,hDCH−Pb,hCH+KQQb,hMT+Qb,hCHP+Qb,hRES+Qb,hESS−γ˜hHb,hCHP+Hb,hBO+Hb,hDRCostu˜≤Δc,(43)Δc=Costb×1+π.Constraints (39) and (40), where Costb is the value of the BC objective function found by equations (1)–(28). In (42), this adopts the uncertainties of the energy price variable and output power of RESs, namely. λ˜,ρ˜,γ˜,P˜RES, as given in (37). In the current study, the uncertainties are modeled using the IGDT method. The model requires one scenario named “worst-case scenario.” As a result, the IGDT method finds the maximum radius of the uncertainty; hence, the problem is feasible. In the end, this is mathematically expressed by equations (41)–(43).Finally, it should be noted that the uncertainty in the active power generation by RESs reduces the flexibility of the network [11, 13]. As a result, the balance between DA and real-time operations may not be met, and active power imbalances may occur in real-time scheduling for the network [16]. Therefore, there is a need to use flexibility resources such as DRP and storage in the network, which can compensate for fluctuations in RESs active power in real-time operation compared to DA operation and enhance the flexibility of the network [13]. To establish such conditions, from the point of view of mathematical modeling, the difference between the active power of the MG in the scenario corresponding to the deterministic model and the worst-case scenario can be minimized as much as possible. This is presented in (44), which should be added to the problem described in (41)–(43). If the flexibility tolerance (∆F) tends to zero, the flexibility will be 100%, and the imbalance between the mentioned operations in the network with RESs will be eliminated. Finally, the flowchart of the problem solution is shown in Figure 2.(44)Pb,hGu˜−Pb,hGu¯≤ΔF∀b,h.Figure 2
Flowchart of the problem solution.


## 3.1. Uncertain Parameters

In the proposed scheme, equations (1)–(28), parameters such as energy price (λ, ρ, and γ), load (PD, QD, and HD), and maximum renewable active power (PRESu) are uncertain. Hence, stochastic [57, 58], probabilistic, or robust [59, 60] modeling is needed for the proposed design. Since the first two types of modeling require a significant number of scenarios and identification of the PDF for uncertainties [61], the computational time is high [9–12]. However, as the implementation step is low in operation problems (up to 15 minutes in real-time operation), the low computational time is of great importance in these problems [16]. Computational time depends on the solution space of the problem [62]. Following this, a robust model of the proposed scheme is presented [5, 17]. In this section, to achieve a robust optimal solution against the maximum prediction error caused by the mentioned uncertainties, a robust model based on IGDT [63] is employed. Also, to investigate the adverse effect of uncertainties in the worst-case scenario, IGDT in accordance with the risk-averse strategy is used for the proposed scheme [63]. In this model, the vector of uncertainty parameters (u¯) is presented as (29), which is the expected (predicted) values of uncertainties. However, the true value of these parameters is a variable in the IGDT model, which following these conditions, the vector of uncertainty variables (u˜) can be expressed as (30). In IGDT, uncertainty is not a parameter, but it is defined as a variable. Also, since this variable is to determine the worst-case scenario this variable is a decision variable in IGDT.(29)u¯=PDQDHDPRESuλργ,(30)u˜=P˜DQ˜DH˜DP˜RESuλ˜ρ˜γ˜.


## 3.2. IGDT-Based Uncertainty Management

An IGDT-based model is used in this study to address uncertainties of demand, price of energy, and output power of RESs. The approach is independent of PDFs. The problem modeling is described as follows:(31)f=minYFY,u˜,(32)LY,u˜≤0,(33)MY,u˜=0,(34)u˜∈Ξ,where, u˜ denotes input uncertainty variables, and Ξ represents the set of uncertainties. Y shows decision-making variables. The present paper used an IGDT based on an information gap model. Even though various models have been used for uncertainty parameters in the literature, the current research employs the envelope-bound model to provide the prior information regarding the uncertainties as given as follows [63]:(35)Ξ=Ξu¯,σ=u˜:u˜−u¯u¯≤σ,where, u˜ represents the uncertain variables matrix, and σ is the difference between the maximum deviation of uncertainty variables and its estimated or nominal value (u¯). The variable σ is known as the uncertainty level and its value is assumed constant for all uncertainties in this paper. Therefore, it has a scalar value. The deterministic model described by equations (31)–(34) is an optimization problem with an estimated or nominal value of uncertainty variables, i.e., equations (31)–(33) together with only u¯.This mathematical expression is called the based condition (BC), where its objective function is given byf=minYFY,u¯. The maximum radius of uncertainty (rc) is found as follows [63]:(36)rc=maxY,σσ.Constraints (32) and (33),(37)FY,u˜≤Δc,(38)Δc=fbY,u¯×1+π,(39)−σ≤u˜−u¯u¯≤σ,(40)0≤σ≤1,where Δc denotes the robustness value of the problem described in equations (31)–(34). The value objective function should be less than Δc. Note that Δc is an input parameter to the problem and is found using (38). In (38), fb is the value of the BC objective function and π∈0⋅1 shows the degree of robustness when fb deteriorates. Consequently, in equations (36)–(40), π is an input parameter and fb shows an input parameter found by equations (31)–(33) taking into account only u¯.In formulation (36)–(40), the uncertainty level, σ, and uncertainty u˜ are variables. So, this problem obtains the highest possible uncertainty level, and the values of uncertainty variables, (30), are determined. The values obtained for the uncertainty variables correspond to the worst-case scenario because, according to (37) and (38), with increasing the value of the uncertainty level, the solution space of the proposed problem decreases compared to the deterministic model; hence, the value of the objective function or cost in these conditions is more than the deterministic model. Thus, the scenario obtained from problem (36)–(40) represents the worst-case scenario.Note that, in the model presented by (36)–(40), the uncertainty level for all uncertainties is assumed to be constant. In other words, it is assumed that the uncertainty radius of all uncertainties is the same. However, the mentioned model has no limitation on considering different uncertainty radii for different uncertainties. Under such conditions, an index such as i is added to σ, which includes uncertainty variables of equation (30). Then, equation (36) is changed as rc=maxY,σ∑iσi, but generally, concerning the IGDT, the uncertainty level is considered fixed for all uncertainties. This can be found in [19].


## 3.3. Robust Formulation of the Proposed Scheme

As per subsection3.2, the IGDT-based robust model proposed for the proposed problem, described by equations (1)–(28), can be given as follows:(41)rc=maxσ.Subject to:Constraints (2)–(28) with u˜ instead of u¯,(42)∑h∈OHλ˜hPrefG+∑b∈BusαbMT+βbMTPb,hMT+χbMTPb,hMT2+αbCHP+βbCHPPb,hCHP+χbCHPPb,hCHP2+αbBO+βbBOHb,hBO+χbBOHb,hBO2−ρ˜hPb,hMT+Pb,hCHP+P˜b,hRES+Pb,hDR+Pb,hDCH−Pb,hCH+KQQb,hMT+Qb,hCHP+Qb,hRES+Qb,hESS−γ˜hHb,hCHP+Hb,hBO+Hb,hDRCostu˜≤Δc,(43)Δc=Costb×1+π.Constraints (39) and (40), where Costb is the value of the BC objective function found by equations (1)–(28). In (42), this adopts the uncertainties of the energy price variable and output power of RESs, namely. λ˜,ρ˜,γ˜,P˜RES, as given in (37). In the current study, the uncertainties are modeled using the IGDT method. The model requires one scenario named “worst-case scenario.” As a result, the IGDT method finds the maximum radius of the uncertainty; hence, the problem is feasible. In the end, this is mathematically expressed by equations (41)–(43).Finally, it should be noted that the uncertainty in the active power generation by RESs reduces the flexibility of the network [11, 13]. As a result, the balance between DA and real-time operations may not be met, and active power imbalances may occur in real-time scheduling for the network [16]. Therefore, there is a need to use flexibility resources such as DRP and storage in the network, which can compensate for fluctuations in RESs active power in real-time operation compared to DA operation and enhance the flexibility of the network [13]. To establish such conditions, from the point of view of mathematical modeling, the difference between the active power of the MG in the scenario corresponding to the deterministic model and the worst-case scenario can be minimized as much as possible. This is presented in (44), which should be added to the problem described in (41)–(43). If the flexibility tolerance (∆F) tends to zero, the flexibility will be 100%, and the imbalance between the mentioned operations in the network with RESs will be eliminated. Finally, the flowchart of the problem solution is shown in Figure 2.(44)Pb,hGu˜−Pb,hGu¯≤ΔF∀b,h.Figure 2
Flowchart of the problem solution.


## 4. Numerical Results

### 4.1. Case Studies

Figure3 depicts a 119-bus test MG used to validate the efficacy of the suggested method [64]. The bass voltage and power of the test system are 11 kV and 10 MVA, respectively, and the voltage value can vary in the range of [0.9 1.05] p.u. [64]. The information about the peak electricity demand and the data on distribution lines and substations are available in [64]. Moreover, the values of hourly electricity demand are found by multiplying the peak demand and hourly electrical load factor curves [5]. The system consists of 7 thermal blocks located at places, as shown in Table 2. Each block has one CHP generation unit, one boiler, and a DRP. The maximum/minimum values of active, reactive, and thermal power of the CHP unit are 800/0 kW, 400/−400 kVAr, and 500/0 kW. The specified values for efficiencies ηT, ηH, ηL are assumed 0.4, 0.4, and 0.08 [5]. The maximum/minimum heat power output from the boiler is 500/0 kW. Coefficients α, β, and χ of the MT, CHP, and boiler are (100, 21, and 0.002), (93, 20, and 0.0017), and (112, 22, and 0.002). The battery-type ESS has an efficiency of 0.9, charge/discharge rate is 150 kW, charger capacity is 200 kVA, and its minimum/maximum permissible energy is 50/500 kWh. The location of the battery is given in Table 2. The block has a thermal load (heat demand) with a peak value of 750 kW. Reference [5] presents the hourly thermal load factor curve of this load. The parameter ξ in the DRP is assumed to be 0.3 [7].Figure 3
119-Bus MG [64].Table 2
Locations of DGs, ESSs, and thermal blocks in the 119-bus MG.
DG, ESS, and blockLocation (bus)PV10, 29, 64WS4, 11, 30, 65, 90, 104MT23, 37, 71, 108ESS4, 10, 30, 64, 104Thermal block6, 15, 51, 81, 97, 110, 115The location of the PV, WS, and MT are given in Table2. The maximum/minimum active and reactive power output of the MT is 600/0 kW and 250/−250 kVAr, and the capacity (maximum apparent power) of the WS and PV system is 600 kVA and 500 kVA, respectively. By multiplying the capacity of RES and the daily generation power rate curve, one can calculate the daily maximum active power output of the RESs [1]. The authors in Ref. [5] provide information about the hourly prices of electrical and heat energies. The parameter KQ is assumed 0.08 [65], and λ = 0.9 × ρ, and ΔF is assumed to be 0.05 p.u. to achieve high flexibility in the network. Note that in this paper choosing λ = 0.9 × ρ is proportional to the assumption that a retailer generally sells the goods purchased from a wholesaler to its retail customers at a higher price.


### 4.2. Results

In this section, the proposed problem is simulated in accordance with the data of the previous section in the MATLAB (GAMS [66–69]) software environment if evolutionary (mathematical) algorithms are used.


#### 4.2.1. Achieving a Suitable Solver

This section presents the convergence results of the proposed scheme for different values of robustness degree (RD orπ) obtained by mathematical solvers such as BARON, BONMIN, DISOPT, KNITRO, and OQNLP [66], and evolutionary algorithms such as teaching-learning-based optimization (TLBO) [70], differential equations (DE) [71], krill herd optimization (KHO) [72], gray wolf optimization (GWO) [73], and antlion optimization (ALO) [74]. To solve the proposed problem by mathematical algorithms, the problem model is coded in the GAMS optimization software, and each of the mentioned algorithms in this software has a toolbox [66]. To solve this problem using evolutionary algorithms, the proposed scheme along with the solution process based on this algorithm are coded in MATLAB software. First off, the algorithm (it can be any of the mentioned algorithms) generates N (population size) random values for decision-making variables including PMT, QMT, QRES, PCH, PDCH, x, QESS, PCHP, QCHP, PDR, HDR, σ, and u˜ in accordance with their operation limits given in equations (11) and (12), [0, SRESu], (15) and (16), {0, 1}, [0, SESSu], (23), (24), (27), [0, 1], and (40). Next, the values of dependent variables such as Cost, PG, QG, PL, QL, V, δ, PRES, E, HBO, and HCHP are calculated from equations (2)–(6), (13), (17), (18), (21), and (22). The values of PRES, E, HBO, and HCHP are determined from constraints (13), (17), (18), (21), and (22), while the values of other dependent variables are calculated from power flow constraints given in (2)–(6). In this paper, the forward-backward method is employed to solve the power flow problem [75]. In the following, the fitness function (FF) is calculated for N values of decision-making variables. To estimate the operation limits, (7)–(10), capacity limit of RESs, (14), the limit of storable energy in ESSs and their charger capacity, (19) and (20), thermal power limit of CHP and boiler, (25) and (26), the DRP limit, (27), and limit of changes of the objective function, (42), the FF in this paper will aim to maximize the difference between the main objective function, (41), and the sum of penalty functions (PeF), i.e., FF = max σ − PeF [14]. The PeF for limits a ≤ b and a = b is given as μ.max (0, a − b) and υ. (b − a), where μ ≥ 0 and υ ∈ (−∞, +∞) represents the Lagrangian multipliers, the values of which are determined like decision-making variables using the evolutionary algorithm [14]. In the following, the steps of updating random variables and Lagrangian multipliers using this algorithm based on their operating limits and the best value of FF in the previous step are implemented. The updating process differs depending on the type of evolutionary algorithm. This process for the TLBO, DE, KHO, GWO, and ALO algorithms can be found in [70–74], respectively. The current study assumes that the convergence point can be accessible after the maximum iterations (Imax). N and Imax are considered 80 and 4000 for evolutionary algorithms, and their other setting parameters were chosen from [70–74].The results of this section are presented in Table3. Based on this table, it can be seen that among the mathematical algorithms, BARON and KNITRO fail to achieve a feasible solution. Among BONMIN, DISOPT, and OQNLP, the BONMIN solver succeeds in obtaining the minimum cost and maximum rc compared to the other two solvers. It has achieved this point in the minimum convergence iterations (CI) and convergence time (CT) in comparison to other mathematical methods. In evolutionary algorithms, ALO has been able to find the most optimal solution (minimum cost and maximum rc) at the minimum CI and CT compared to TLBO, DE, KHO, and GWO. Furthermore, by comparing mathematical solution methods and evolutionary algorithms, according to Table 3, it is observed that ALO has more favorable conditions than BONMIN in terms of response speed and optimal point. This is because it has the minimum cost (maximum rc) compared to BONMIN. The CI of ALO is higher than BONMIN, but it has reached the optimal solution with a much lower CT than BONMIN. Therefore, in this paper, the ALO algorithm is a more suitable solver than the other solvers mentioned.Table 3
Convergence situation achieved by different solvers for different values of robustness degree.
π00.1SolverCost ($)rcCICT (min)Cost ($)rcCICT (min)BARONInfeasibleBONMIN−1892.1030739.7−113.50.10933842.8DISOPT−1611.4038451.4−78.20.10542255.2KNITROInfeasibleOQNLP−1537.5045767.2−75.70.10249772.5TLBO−1761.80242215.1−87.30.114267517.9DE−1756.30254915.6−84.50.111283118.4KHO−1812.50206112.3−99.70.124219714.1GWO−1838.60193811.5−104.20.127209513.2ALO−1976.3016729.2−139.80.134175610.1In Table4, the convergence status of the ALO-based problem using different uncertainty modeling methods for different values is presented. In this table, it can be seen that in the uncertainty modeling method based on scenario-based stochastic programming (SBSP) [15], a large number of scenarios are needed to find a reliable optimal solution. More than 80 scenarios show that the optimal solution has one solution (cost is constant), and the computational time and convergence iteration change. Therefore, it can be stated that 80 scenarios are needed to access a reliable solution in this method. The UT method [13] requires a total of 15 scenarios because, in (30), the number of uncertainty parameters is 7 and the total number of scenarios required in UT corresponding to 2n + 1 will be equal to 15 (2 × 7 + 1) scenarios. Therefore, with this number of scenarios, it has found an optimal solution such as SBSP in 80 scenarios, but its computational time is lower than SBSP. In the following, robust methods such as ARO [50] and IGDT have only one scenario. Of course, for π = 0, they consider deterministic modeling for uncertainties. Therefore, the value of the objective function in these methods is less than SBSP and UT, because the cost of modeling uncertainties is not obtained in these conditions. They also have one scenario, so their computational time is less than that of SBSP and UT. According to [50], the ARO formulation for a problem will complicate the problem. This makes the computation time longer than IGDT. Then, by increasing π, it is seen that the highest rc is obtained in IGDT. In ARO, the maximum radial uncertainty is obtained by the trial and error method and by changing the level of uncertainty. Hence, its computational time is longer than IGDT. Note that the details of SBSP, UT, and ARO are given in [13, 15], and [50], respectively. Based on equation (1), the term cost is equal to the difference between the sum of the operating cost of the network and nonrenewable sources and the sum of the revenues that power sources, storage, and responsive loads earn by selling electrical and thermal energy. With optimal management of sources, storage, and responsive loads (with their performance curves given in subsection 4.2.3), in the proposed scheme, it is inferred that the revenue of the mentioned elements overcomes their operating costs and that of the network (i.e., the revenue is greater than the operating cost). In these conditions, the value of cost in Tables 3 and 4 has become negative.Table 4
Convergence situation based on ALO in the different uncertainty modeling method for various values of robustness degree.
π00.1MethodCost ($)rcCICT (min)Cost ($)rcCICT (min)SBSP40 scenario−1669.8—263224.2There is no this condition60 scenarios−1663.3—286727.580 scenario−1657.1—307130.182 scenario−1657.1—307930.5UT (15 scenarios)−1657.6—243117.6ARO (1 scenario)−1976.3—182910.9−143.20.129216238.1 due to including the trial and error methodIGDT (1 scenario)−1976.3—16729.2−139.80.134175610.1


#### 4.2.2. Investigation of Uncertainty Parameters in the Worst-Case Scenario

Table5 reports the true value of uncertainty variables for different values of RD (π). According to this table, it is observed that increasing the degree of robustness to 0.4 increases the maximum radius of uncertainty (maximum forecast error or rc) to 0.345. Also, in the conditions of increasing RD with respect to 0.4, the value of rc is always constant. This indicates the fact that the implementation of the proposed scheme on the data in subsection 4.1 can cover the prediction error of 34.5% resulting from the uncertainty parameters including load, energy price, and renewable power. If the prediction error (radius of uncertainty) is more than 34.5%, the problem has no feasible solution. Moreover, it is observed that with increasing the maximum radius of uncertainty, the amount of load and price of energy purchased from the upstream network in the worst-case scenario have an upward trend, but the power generation by renewables and energy selling price have a downward trend in this scenario. So, in general, it can be stated that in the worst-case scenario, compared to the scenario corresponding to the deterministic model (π = 0), the amount of energy consumed by the network (energy purchase price) increases while the energy produced in the network (energy sales price) decreases.Table 5
True value of uncertainty parameters at the worst-case scenario for different values of robustness degree.
π00.10.20.30.40.5rc00.1340.2210.3020.3450.345∑b,hP˜b,hD (MW)429.9487.5524.9559.7578.2578.2∑b,hQ˜b,hD (MVAr)322.6365.8393.9420433.9433.9∑b,hH˜b,hD (MW)104.3118.3127.3135.8140.3140.3∑b,hP˜b,hRESu (MW)79.969.262.255.852.352.3∑hλ˜h ($/MWh)611.6693.5746.8796.3822.6822.6∑hρ˜h ($/MWh)550.4476.6428.8384.2360.5360.5∑hγ˜h ($/MWh)616533.5479.9430403.5403.5


#### 4.2.3. Investigation of the Operation and Economic Status of Distributed Energy Resources (DERs)

The daily curves of active and reactive power of sources, storage devices, and electrically responsive loads for different values of RD are shown in Figures4 and 5, respectively. Based on Figure 4(a), the data in subsections 4.1 and [1], corresponding to RESs, it can be seen that the active injection power of WSs and PVs at all operating hours is equal to the maximum active power they produce in accordance with the weather conditions. In addition, by increasing the RD, the amount of active power injection by RESs into the MG decreases in the worst-case scenario due to reducing P˜RESu in this condition based on Table 5. Figure 4(b) shows the daily power curves of CHPs and MTs. According to this figure, these DGs inject the maximum active power they can produce into the network during all operating hours. This is because (1) the energy purchase price (λ) is higher than the fuel price of DGs in some hours, hence the MG operator tends to utilize these DGs to supply energy consumption, and (2) these DGs sell energy to the MG at a price ρ, which is about 111% (100%/90%) of λ based on subsection 4.1. Therefore, to minimize the Cost function in equation (1), DGs are expected to inject as much energy as possible into the network, so they will always operate at the point corresponding to the maximum active output power. Even in the worst-case scenario (RD = 0.2), where λ (ρ) has increased (decreased) by about 22% compared to RD = 0 according to Table 5, the performance results of DGs are the same as deterministic studies (RD = 0). Regarding the performance of DRPs and ESSs in Figures 4(c) and 4(d), they perform charging operations during 1 : 00–16 : 00 and 23 : 00–00 : 00 so that during 1 : 00–7:00, high energy is received by DRP and ESSs from MG. They also operate in discharge mode from 17 : 00 to 22 : 00 and inject active power into the MG. The performance of these elements corresponds to the electricity price. According to [5], energy prices in the hours 1 : 00–7:00, 8 : 00–16 : 00, and 23 : 00–00 : 00, and 17 : 00–22 : 00 are equal to 17.6 $/MWh, 26.4 $/MWh, and 33 $/MWh. Hence, to minimize the Cost function, they will operate in charging (discharging) mode at times corresponding to the cheap (expensive) energy price. As another point, note that since the active power generation by RESs decreases in the worst-case scenario according to Figure 4(a), DRPs and ESSs inject a higher active power in this scenario during peak hours (17 : 00–22 : 00) into MG than in the case with RD = 0. They also increase their power consumption during other hours to be able to inject high energy into the network during other hours. Finally, note that CHP, MT, ESSs, and DRP are flexibility sources for RES. As RESs generate active power at all operating hours, it is expected that MG flexibility control is required at all hours. Therefore, flexible sources according to Figure 4 are not switched off at any time of operation.Figure 4
Daily active power curves of (a) RESs, (b) non-RESs, (c) DRPs, and (d) ESSs at different values of robustness degree.
(a)(b)(c)(d)Figure 5
Daily reactive power curve of (a) PVs, (b) WSs, (c) non-RESs, and (d) ESSs in the different values of robustness degree.
(a)(b)(c)(d)The daily reactive power curves of DGs and ESSs for different RD values are shown in Figure5. According to this figure, DGs inject the maximum active power into the network during most of the operating hours. For example, according to the data in subsection 4.1, the maximum capacity (apparent power) of PVs is 1.5 MVA (3 PV × 0.5 MVA for PV capacity), and PVs do not produce active power during 1 : 00–5:00 and 1900–00 : 00 according to Figure 4(a). Therefore, they inject reactive power of 1.5 MVAr into the network during these hours. However, in other hours, when their active power generation increases, their reactive power generation decreases. Moreover, with increasing RD, the active power generation by PVs decreases according to Figure 4(a), so they inject more reactive power into MG than in the case with RD = 0. These conditions hold for WSs according to Figure 5(b). But CHPs (MTs) injected the maximum reactive power that can be produced, i.e., 7 × 0.4 = 2.8 MVAr (4 × 0.25 = 1 MVAr) into the network at all hours and for all values of RD. In the following, comparing the load profile presented in [17] and Figure 5(d), it is observed that the daily reactive power curve of ESSs is the same as the daily reactive load profile. Since the amount of reactive power consumed by MG increases if RD is increased (Table 5), ESSs produce high reactive power. Note that this performance of sources and ESSs in the MG reactive power supply process corresponds to maximizing its revenue in the reactive power market based on equation (1).Figure6 depicts the thermal capacity balance curve of the sources and the response loads of the thermal block for different values of RD. According to Figure 6(a), thermal DRP reduces thermal energy consumption during peak thermal hours, i.e., 5 : 00–15 : 00, but at other times, it increases the energy consumption of thermal consumers. According to [5], the price of thermal energy in the hours 5 : 00–15 : 00 is 30 $/MWh and 22 $/MWh at other hours. Hence, the performance of thermal DRP in the figure corresponds to the minimization of the Cost function in (1). In addition, since the active power generation by CHP was constant at all operating hours according to Figure 3(b), the amount of heat generated by CHP, according to (22), will be constant for time changes and will have a daily curve as shown in Figure 6(b). Finally, according to (21), the boiler heat output will be equal to the difference between the thermal load and the total heat output of CHP and thermal DRP, so its daily profile will be as shown in Figure 6(b). Then, with increasing RD, the daily heat power curve of the boiler and DRP shifts upwards relative to the case RD = 0, because in the daily curve, the thermal load shifts upwards in the worst-case scenario according to Figure 6(a). Nonetheless, because the thermal power of CHP depends on its active power, and the active power of CHP, according to Figure 4(b), always has a constant value for different values of RD, the thermal power of CHP will not change with changes in RD. Finally, the economic indices curves of MG and DERs, including sources, storage devices, and responsive loads, are plotted in terms of RD in Figure 7. According to this figure, increasing the RD to 0.4 increases the operating costs of MGs and NRDGs, while the revenue of DERs decreases under these conditions. This leads to an increase in the net network cost (Cost function), i.e., the difference between the total costs of MG and NRDGs and the revenue of DERs in the energy and reactive power markets under the mentioned conditions. Nonetheless, the values of the economic indices mentioned for RD > 0.4 are the same as for RD = 0.4. All these results obtained in Figure 6 are proportional to the change in the status of energy consumption/generation and the purchase/sale price of energy in the worst-case scenario in Table 5 or subsection 4.2.2.Figure 6
Daily heat power curve of (a) block load without/with DRP and (b) sources in the different values of robustness degree.
(a)(b)Figure 7
Curve of economic indices in robustness degree.


#### 4.2.4. Evaluation of the Economic and Technical Status of MG

In this section, six case studies are analyzed to examine the capabilities of the proposed scheme, as follows:(i)
Case I: Power flow analysis (MG without sources, storage devices, DRP, and thermal block)(ii)
Case II: Proposed scheme considering only RESs(iii)
Case III: Proposed scheme considering only ESSs(iv)
Case IV: Proposed scheme considering only DRPs(v)
Case V: Proposed scheme considering NRDGs (MTs and thermal blocks)(vi)
Case VI: Proposed scheme considering all proposed DERsThe results of this section are shown in Table6, which reports the values of cost, energy loss (EL), minimum power factor (MPF) of the distribution substation, maximum voltage drop (MVD), and maximum overvoltage (MOV) in cases I to VI for different values of RD. According to this table, the presence of RESs alone in the proposed scheme (II) or the single presence of thermal blocks and MT in MG (V) leads to a significant reduction in energy costs, energy losses, and the voltage drop of MG compared to Case I. Also, MPF increases significantly; in Case V it has increased by more than 0.9. These conditions correspond to the occurrence of MOV around 0.02 p.u., which is less than its allowable limit, i.e., 0.05 (1.05–1). The performance of ESSs in Case III has a favorable effect on improving the status of EL, MPF, and MVD and a low impact on cost reduction without causing overvoltage compared to power flow studies, but the percentage of improvement of these parameters, in this case, is less than cases II and V because in cases II and V, high active and reactive power are simultaneously injected into the network by DGs based on Figures 4 and 5. But, in Case III, ESSs inject active power into the network only during peak hours and consume active power during other hours. The ESSs also failed to generate as much reactive power as the DGs; this is because of their low charger capacity and low number compared to the DGs in the 119-bus network based on subsection 4.1. In the case of DRP, it has been able to enhance cost, EL, and MVD in Case IV compared to Case I without causing overvoltage but has led to MPF degradation. Note that the percentage of improvement of cost, EL, and MVD in this case study is lower than in cases II and V but higher than in Case III. In addition, in Case IV, DRP only plays a role in controlling the active power and does not change the amount of reactive power. Therefore, according to Figure 4(c), during peak hours when the amount of active power consumption (active power of the distribution substation) is reduced by DRP with no change in reactive power status compared to Case I, based on equation (9), the power factor of the distribution substation will be reduced compared to Case I. Finally, the utilization of all DERs in the proposed scheme (VI) has enhanced all the proposed indices compared to Case I. In the worst-case scenario (RD = 0.1), the Cost function in Case VI is reduced by about 101% compared to Case I (12341-(-139.8))/12341). This value (percentage of improvement) for EL, MPF, and MVD in the worst-case scenario with a MOV of 0.023 p.u. is about 43.8%, 15.9%, and 40.7%, respectively. Note that according to Table 6 for an increase in RD, all the mentioned indices except for MOV and MPF increase compared to RD = 0, but MOV and MPF decrease. The reason is that, according to Table 5 in the worst-case scenario compared to the deterministic model, energy consumption (generation) increases (decreases). Finally, the cost curve in terms of flexibility tolerance (∆F) for different RD values is plotted in Figure 8. Based on this figure, it can be seen that with increasing ∆F, Cost decreases because, in these conditions, the operating cost of resources, storage devices, and responsive loads is reduced. For example, if ∆F increases, then the importance of network flexibility decreases. Therefore, electric DRPs and ESSs will be switched off at 8 : 00–16 : 00 and 23 : 00–00 : 00 and will not receive power from the network. This results in a reduction in their operating costs. As the RD increases, the curve shifts upwards, because in these conditions the energy consumption/generation and the purchase/sale price of energy increase/decrease compared to the case with RD = 0 according to Table 5.Table 6
Values of operation and economic indices in the different cases for various values of robustness degree.
RD00.1CaseIIIIIIIVVVIIIIIIIIVVVICost1119572571036286395376−1976.31234181211096492565376−139.8EL (MWh)7.2335.8926.9226.4625.1314.0487.9516.4457.6117.0235.5844.465MPF0.7990.8730.8350.7560.9080.9310.7990.8610.8300.7510.9020.926MVD (p.u)0.0920.0670.0840.0710.0640.0590.1080.0720.0910.0790.0700.064MOV (p.u)00.021000.0240.02600.018000.0220.023Figure 8
Curve of cost in ΔF in different values of RD.


## 4.1. Case Studies

Figure3 depicts a 119-bus test MG used to validate the efficacy of the suggested method [64]. The bass voltage and power of the test system are 11 kV and 10 MVA, respectively, and the voltage value can vary in the range of [0.9 1.05] p.u. [64]. The information about the peak electricity demand and the data on distribution lines and substations are available in [64]. Moreover, the values of hourly electricity demand are found by multiplying the peak demand and hourly electrical load factor curves [5]. The system consists of 7 thermal blocks located at places, as shown in Table 2. Each block has one CHP generation unit, one boiler, and a DRP. The maximum/minimum values of active, reactive, and thermal power of the CHP unit are 800/0 kW, 400/−400 kVAr, and 500/0 kW. The specified values for efficiencies ηT, ηH, ηL are assumed 0.4, 0.4, and 0.08 [5]. The maximum/minimum heat power output from the boiler is 500/0 kW. Coefficients α, β, and χ of the MT, CHP, and boiler are (100, 21, and 0.002), (93, 20, and 0.0017), and (112, 22, and 0.002). The battery-type ESS has an efficiency of 0.9, charge/discharge rate is 150 kW, charger capacity is 200 kVA, and its minimum/maximum permissible energy is 50/500 kWh. The location of the battery is given in Table 2. The block has a thermal load (heat demand) with a peak value of 750 kW. Reference [5] presents the hourly thermal load factor curve of this load. The parameter ξ in the DRP is assumed to be 0.3 [7].Figure 3
119-Bus MG [64].Table 2
Locations of DGs, ESSs, and thermal blocks in the 119-bus MG.
DG, ESS, and blockLocation (bus)PV10, 29, 64WS4, 11, 30, 65, 90, 104MT23, 37, 71, 108ESS4, 10, 30, 64, 104Thermal block6, 15, 51, 81, 97, 110, 115The location of the PV, WS, and MT are given in Table2. The maximum/minimum active and reactive power output of the MT is 600/0 kW and 250/−250 kVAr, and the capacity (maximum apparent power) of the WS and PV system is 600 kVA and 500 kVA, respectively. By multiplying the capacity of RES and the daily generation power rate curve, one can calculate the daily maximum active power output of the RESs [1]. The authors in Ref. [5] provide information about the hourly prices of electrical and heat energies. The parameter KQ is assumed 0.08 [65], and λ = 0.9 × ρ, and ΔF is assumed to be 0.05 p.u. to achieve high flexibility in the network. Note that in this paper choosing λ = 0.9 × ρ is proportional to the assumption that a retailer generally sells the goods purchased from a wholesaler to its retail customers at a higher price.


## 4.2. Results

In this section, the proposed problem is simulated in accordance with the data of the previous section in the MATLAB (GAMS [66–69]) software environment if evolutionary (mathematical) algorithms are used.


### 4.2.1. Achieving a Suitable Solver

This section presents the convergence results of the proposed scheme for different values of robustness degree (RD orπ) obtained by mathematical solvers such as BARON, BONMIN, DISOPT, KNITRO, and OQNLP [66], and evolutionary algorithms such as teaching-learning-based optimization (TLBO) [70], differential equations (DE) [71], krill herd optimization (KHO) [72], gray wolf optimization (GWO) [73], and antlion optimization (ALO) [74]. To solve the proposed problem by mathematical algorithms, the problem model is coded in the GAMS optimization software, and each of the mentioned algorithms in this software has a toolbox [66]. To solve this problem using evolutionary algorithms, the proposed scheme along with the solution process based on this algorithm are coded in MATLAB software. First off, the algorithm (it can be any of the mentioned algorithms) generates N (population size) random values for decision-making variables including PMT, QMT, QRES, PCH, PDCH, x, QESS, PCHP, QCHP, PDR, HDR, σ, and u˜ in accordance with their operation limits given in equations (11) and (12), [0, SRESu], (15) and (16), {0, 1}, [0, SESSu], (23), (24), (27), [0, 1], and (40). Next, the values of dependent variables such as Cost, PG, QG, PL, QL, V, δ, PRES, E, HBO, and HCHP are calculated from equations (2)–(6), (13), (17), (18), (21), and (22). The values of PRES, E, HBO, and HCHP are determined from constraints (13), (17), (18), (21), and (22), while the values of other dependent variables are calculated from power flow constraints given in (2)–(6). In this paper, the forward-backward method is employed to solve the power flow problem [75]. In the following, the fitness function (FF) is calculated for N values of decision-making variables. To estimate the operation limits, (7)–(10), capacity limit of RESs, (14), the limit of storable energy in ESSs and their charger capacity, (19) and (20), thermal power limit of CHP and boiler, (25) and (26), the DRP limit, (27), and limit of changes of the objective function, (42), the FF in this paper will aim to maximize the difference between the main objective function, (41), and the sum of penalty functions (PeF), i.e., FF = max σ − PeF [14]. The PeF for limits a ≤ b and a = b is given as μ.max (0, a − b) and υ. (b − a), where μ ≥ 0 and υ ∈ (−∞, +∞) represents the Lagrangian multipliers, the values of which are determined like decision-making variables using the evolutionary algorithm [14]. In the following, the steps of updating random variables and Lagrangian multipliers using this algorithm based on their operating limits and the best value of FF in the previous step are implemented. The updating process differs depending on the type of evolutionary algorithm. This process for the TLBO, DE, KHO, GWO, and ALO algorithms can be found in [70–74], respectively. The current study assumes that the convergence point can be accessible after the maximum iterations (Imax). N and Imax are considered 80 and 4000 for evolutionary algorithms, and their other setting parameters were chosen from [70–74].The results of this section are presented in Table3. Based on this table, it can be seen that among the mathematical algorithms, BARON and KNITRO fail to achieve a feasible solution. Among BONMIN, DISOPT, and OQNLP, the BONMIN solver succeeds in obtaining the minimum cost and maximum rc compared to the other two solvers. It has achieved this point in the minimum convergence iterations (CI) and convergence time (CT) in comparison to other mathematical methods. In evolutionary algorithms, ALO has been able to find the most optimal solution (minimum cost and maximum rc) at the minimum CI and CT compared to TLBO, DE, KHO, and GWO. Furthermore, by comparing mathematical solution methods and evolutionary algorithms, according to Table 3, it is observed that ALO has more favorable conditions than BONMIN in terms of response speed and optimal point. This is because it has the minimum cost (maximum rc) compared to BONMIN. The CI of ALO is higher than BONMIN, but it has reached the optimal solution with a much lower CT than BONMIN. Therefore, in this paper, the ALO algorithm is a more suitable solver than the other solvers mentioned.Table 3
Convergence situation achieved by different solvers for different values of robustness degree.
π00.1SolverCost ($)rcCICT (min)Cost ($)rcCICT (min)BARONInfeasibleBONMIN−1892.1030739.7−113.50.10933842.8DISOPT−1611.4038451.4−78.20.10542255.2KNITROInfeasibleOQNLP−1537.5045767.2−75.70.10249772.5TLBO−1761.80242215.1−87.30.114267517.9DE−1756.30254915.6−84.50.111283118.4KHO−1812.50206112.3−99.70.124219714.1GWO−1838.60193811.5−104.20.127209513.2ALO−1976.3016729.2−139.80.134175610.1In Table4, the convergence status of the ALO-based problem using different uncertainty modeling methods for different values is presented. In this table, it can be seen that in the uncertainty modeling method based on scenario-based stochastic programming (SBSP) [15], a large number of scenarios are needed to find a reliable optimal solution. More than 80 scenarios show that the optimal solution has one solution (cost is constant), and the computational time and convergence iteration change. Therefore, it can be stated that 80 scenarios are needed to access a reliable solution in this method. The UT method [13] requires a total of 15 scenarios because, in (30), the number of uncertainty parameters is 7 and the total number of scenarios required in UT corresponding to 2n + 1 will be equal to 15 (2 × 7 + 1) scenarios. Therefore, with this number of scenarios, it has found an optimal solution such as SBSP in 80 scenarios, but its computational time is lower than SBSP. In the following, robust methods such as ARO [50] and IGDT have only one scenario. Of course, for π = 0, they consider deterministic modeling for uncertainties. Therefore, the value of the objective function in these methods is less than SBSP and UT, because the cost of modeling uncertainties is not obtained in these conditions. They also have one scenario, so their computational time is less than that of SBSP and UT. According to [50], the ARO formulation for a problem will complicate the problem. This makes the computation time longer than IGDT. Then, by increasing π, it is seen that the highest rc is obtained in IGDT. In ARO, the maximum radial uncertainty is obtained by the trial and error method and by changing the level of uncertainty. Hence, its computational time is longer than IGDT. Note that the details of SBSP, UT, and ARO are given in [13, 15], and [50], respectively. Based on equation (1), the term cost is equal to the difference between the sum of the operating cost of the network and nonrenewable sources and the sum of the revenues that power sources, storage, and responsive loads earn by selling electrical and thermal energy. With optimal management of sources, storage, and responsive loads (with their performance curves given in subsection 4.2.3), in the proposed scheme, it is inferred that the revenue of the mentioned elements overcomes their operating costs and that of the network (i.e., the revenue is greater than the operating cost). In these conditions, the value of cost in Tables 3 and 4 has become negative.Table 4
Convergence situation based on ALO in the different uncertainty modeling method for various values of robustness degree.
π00.1MethodCost ($)rcCICT (min)Cost ($)rcCICT (min)SBSP40 scenario−1669.8—263224.2There is no this condition60 scenarios−1663.3—286727.580 scenario−1657.1—307130.182 scenario−1657.1—307930.5UT (15 scenarios)−1657.6—243117.6ARO (1 scenario)−1976.3—182910.9−143.20.129216238.1 due to including the trial and error methodIGDT (1 scenario)−1976.3—16729.2−139.80.134175610.1


### 4.2.2. Investigation of Uncertainty Parameters in the Worst-Case Scenario

Table5 reports the true value of uncertainty variables for different values of RD (π). According to this table, it is observed that increasing the degree of robustness to 0.4 increases the maximum radius of uncertainty (maximum forecast error or rc) to 0.345. Also, in the conditions of increasing RD with respect to 0.4, the value of rc is always constant. This indicates the fact that the implementation of the proposed scheme on the data in subsection 4.1 can cover the prediction error of 34.5% resulting from the uncertainty parameters including load, energy price, and renewable power. If the prediction error (radius of uncertainty) is more than 34.5%, the problem has no feasible solution. Moreover, it is observed that with increasing the maximum radius of uncertainty, the amount of load and price of energy purchased from the upstream network in the worst-case scenario have an upward trend, but the power generation by renewables and energy selling price have a downward trend in this scenario. So, in general, it can be stated that in the worst-case scenario, compared to the scenario corresponding to the deterministic model (π = 0), the amount of energy consumed by the network (energy purchase price) increases while the energy produced in the network (energy sales price) decreases.Table 5
True value of uncertainty parameters at the worst-case scenario for different values of robustness degree.
π00.10.20.30.40.5rc00.1340.2210.3020.3450.345∑b,hP˜b,hD (MW)429.9487.5524.9559.7578.2578.2∑b,hQ˜b,hD (MVAr)322.6365.8393.9420433.9433.9∑b,hH˜b,hD (MW)104.3118.3127.3135.8140.3140.3∑b,hP˜b,hRESu (MW)79.969.262.255.852.352.3∑hλ˜h ($/MWh)611.6693.5746.8796.3822.6822.6∑hρ˜h ($/MWh)550.4476.6428.8384.2360.5360.5∑hγ˜h ($/MWh)616533.5479.9430403.5403.5


### 4.2.3. Investigation of the Operation and Economic Status of Distributed Energy Resources (DERs)

The daily curves of active and reactive power of sources, storage devices, and electrically responsive loads for different values of RD are shown in Figures4 and 5, respectively. Based on Figure 4(a), the data in subsections 4.1 and [1], corresponding to RESs, it can be seen that the active injection power of WSs and PVs at all operating hours is equal to the maximum active power they produce in accordance with the weather conditions. In addition, by increasing the RD, the amount of active power injection by RESs into the MG decreases in the worst-case scenario due to reducing P˜RESu in this condition based on Table 5. Figure 4(b) shows the daily power curves of CHPs and MTs. According to this figure, these DGs inject the maximum active power they can produce into the network during all operating hours. This is because (1) the energy purchase price (λ) is higher than the fuel price of DGs in some hours, hence the MG operator tends to utilize these DGs to supply energy consumption, and (2) these DGs sell energy to the MG at a price ρ, which is about 111% (100%/90%) of λ based on subsection 4.1. Therefore, to minimize the Cost function in equation (1), DGs are expected to inject as much energy as possible into the network, so they will always operate at the point corresponding to the maximum active output power. Even in the worst-case scenario (RD = 0.2), where λ (ρ) has increased (decreased) by about 22% compared to RD = 0 according to Table 5, the performance results of DGs are the same as deterministic studies (RD = 0). Regarding the performance of DRPs and ESSs in Figures 4(c) and 4(d), they perform charging operations during 1 : 00–16 : 00 and 23 : 00–00 : 00 so that during 1 : 00–7:00, high energy is received by DRP and ESSs from MG. They also operate in discharge mode from 17 : 00 to 22 : 00 and inject active power into the MG. The performance of these elements corresponds to the electricity price. According to [5], energy prices in the hours 1 : 00–7:00, 8 : 00–16 : 00, and 23 : 00–00 : 00, and 17 : 00–22 : 00 are equal to 17.6 $/MWh, 26.4 $/MWh, and 33 $/MWh. Hence, to minimize the Cost function, they will operate in charging (discharging) mode at times corresponding to the cheap (expensive) energy price. As another point, note that since the active power generation by RESs decreases in the worst-case scenario according to Figure 4(a), DRPs and ESSs inject a higher active power in this scenario during peak hours (17 : 00–22 : 00) into MG than in the case with RD = 0. They also increase their power consumption during other hours to be able to inject high energy into the network during other hours. Finally, note that CHP, MT, ESSs, and DRP are flexibility sources for RES. As RESs generate active power at all operating hours, it is expected that MG flexibility control is required at all hours. Therefore, flexible sources according to Figure 4 are not switched off at any time of operation.Figure 4
Daily active power curves of (a) RESs, (b) non-RESs, (c) DRPs, and (d) ESSs at different values of robustness degree.
(a)(b)(c)(d)Figure 5
Daily reactive power curve of (a) PVs, (b) WSs, (c) non-RESs, and (d) ESSs in the different values of robustness degree.
(a)(b)(c)(d)The daily reactive power curves of DGs and ESSs for different RD values are shown in Figure5. According to this figure, DGs inject the maximum active power into the network during most of the operating hours. For example, according to the data in subsection 4.1, the maximum capacity (apparent power) of PVs is 1.5 MVA (3 PV × 0.5 MVA for PV capacity), and PVs do not produce active power during 1 : 00–5:00 and 1900–00 : 00 according to Figure 4(a). Therefore, they inject reactive power of 1.5 MVAr into the network during these hours. However, in other hours, when their active power generation increases, their reactive power generation decreases. Moreover, with increasing RD, the active power generation by PVs decreases according to Figure 4(a), so they inject more reactive power into MG than in the case with RD = 0. These conditions hold for WSs according to Figure 5(b). But CHPs (MTs) injected the maximum reactive power that can be produced, i.e., 7 × 0.4 = 2.8 MVAr (4 × 0.25 = 1 MVAr) into the network at all hours and for all values of RD. In the following, comparing the load profile presented in [17] and Figure 5(d), it is observed that the daily reactive power curve of ESSs is the same as the daily reactive load profile. Since the amount of reactive power consumed by MG increases if RD is increased (Table 5), ESSs produce high reactive power. Note that this performance of sources and ESSs in the MG reactive power supply process corresponds to maximizing its revenue in the reactive power market based on equation (1).Figure6 depicts the thermal capacity balance curve of the sources and the response loads of the thermal block for different values of RD. According to Figure 6(a), thermal DRP reduces thermal energy consumption during peak thermal hours, i.e., 5 : 00–15 : 00, but at other times, it increases the energy consumption of thermal consumers. According to [5], the price of thermal energy in the hours 5 : 00–15 : 00 is 30 $/MWh and 22 $/MWh at other hours. Hence, the performance of thermal DRP in the figure corresponds to the minimization of the Cost function in (1). In addition, since the active power generation by CHP was constant at all operating hours according to Figure 3(b), the amount of heat generated by CHP, according to (22), will be constant for time changes and will have a daily curve as shown in Figure 6(b). Finally, according to (21), the boiler heat output will be equal to the difference between the thermal load and the total heat output of CHP and thermal DRP, so its daily profile will be as shown in Figure 6(b). Then, with increasing RD, the daily heat power curve of the boiler and DRP shifts upwards relative to the case RD = 0, because in the daily curve, the thermal load shifts upwards in the worst-case scenario according to Figure 6(a). Nonetheless, because the thermal power of CHP depends on its active power, and the active power of CHP, according to Figure 4(b), always has a constant value for different values of RD, the thermal power of CHP will not change with changes in RD. Finally, the economic indices curves of MG and DERs, including sources, storage devices, and responsive loads, are plotted in terms of RD in Figure 7. According to this figure, increasing the RD to 0.4 increases the operating costs of MGs and NRDGs, while the revenue of DERs decreases under these conditions. This leads to an increase in the net network cost (Cost function), i.e., the difference between the total costs of MG and NRDGs and the revenue of DERs in the energy and reactive power markets under the mentioned conditions. Nonetheless, the values of the economic indices mentioned for RD > 0.4 are the same as for RD = 0.4. All these results obtained in Figure 6 are proportional to the change in the status of energy consumption/generation and the purchase/sale price of energy in the worst-case scenario in Table 5 or subsection 4.2.2.Figure 6
Daily heat power curve of (a) block load without/with DRP and (b) sources in the different values of robustness degree.
(a)(b)Figure 7
Curve of economic indices in robustness degree.


### 4.2.4. Evaluation of the Economic and Technical Status of MG

In this section, six case studies are analyzed to examine the capabilities of the proposed scheme, as follows:(i)
Case I: Power flow analysis (MG without sources, storage devices, DRP, and thermal block)(ii)
Case II: Proposed scheme considering only RESs(iii)
Case III: Proposed scheme considering only ESSs(iv)
Case IV: Proposed scheme considering only DRPs(v)
Case V: Proposed scheme considering NRDGs (MTs and thermal blocks)(vi)
Case VI: Proposed scheme considering all proposed DERsThe results of this section are shown in Table6, which reports the values of cost, energy loss (EL), minimum power factor (MPF) of the distribution substation, maximum voltage drop (MVD), and maximum overvoltage (MOV) in cases I to VI for different values of RD. According to this table, the presence of RESs alone in the proposed scheme (II) or the single presence of thermal blocks and MT in MG (V) leads to a significant reduction in energy costs, energy losses, and the voltage drop of MG compared to Case I. Also, MPF increases significantly; in Case V it has increased by more than 0.9. These conditions correspond to the occurrence of MOV around 0.02 p.u., which is less than its allowable limit, i.e., 0.05 (1.05–1). The performance of ESSs in Case III has a favorable effect on improving the status of EL, MPF, and MVD and a low impact on cost reduction without causing overvoltage compared to power flow studies, but the percentage of improvement of these parameters, in this case, is less than cases II and V because in cases II and V, high active and reactive power are simultaneously injected into the network by DGs based on Figures 4 and 5. But, in Case III, ESSs inject active power into the network only during peak hours and consume active power during other hours. The ESSs also failed to generate as much reactive power as the DGs; this is because of their low charger capacity and low number compared to the DGs in the 119-bus network based on subsection 4.1. In the case of DRP, it has been able to enhance cost, EL, and MVD in Case IV compared to Case I without causing overvoltage but has led to MPF degradation. Note that the percentage of improvement of cost, EL, and MVD in this case study is lower than in cases II and V but higher than in Case III. In addition, in Case IV, DRP only plays a role in controlling the active power and does not change the amount of reactive power. Therefore, according to Figure 4(c), during peak hours when the amount of active power consumption (active power of the distribution substation) is reduced by DRP with no change in reactive power status compared to Case I, based on equation (9), the power factor of the distribution substation will be reduced compared to Case I. Finally, the utilization of all DERs in the proposed scheme (VI) has enhanced all the proposed indices compared to Case I. In the worst-case scenario (RD = 0.1), the Cost function in Case VI is reduced by about 101% compared to Case I (12341-(-139.8))/12341). This value (percentage of improvement) for EL, MPF, and MVD in the worst-case scenario with a MOV of 0.023 p.u. is about 43.8%, 15.9%, and 40.7%, respectively. Note that according to Table 6 for an increase in RD, all the mentioned indices except for MOV and MPF increase compared to RD = 0, but MOV and MPF decrease. The reason is that, according to Table 5 in the worst-case scenario compared to the deterministic model, energy consumption (generation) increases (decreases). Finally, the cost curve in terms of flexibility tolerance (∆F) for different RD values is plotted in Figure 8. Based on this figure, it can be seen that with increasing ∆F, Cost decreases because, in these conditions, the operating cost of resources, storage devices, and responsive loads is reduced. For example, if ∆F increases, then the importance of network flexibility decreases. Therefore, electric DRPs and ESSs will be switched off at 8 : 00–16 : 00 and 23 : 00–00 : 00 and will not receive power from the network. This results in a reduction in their operating costs. As the RD increases, the curve shifts upwards, because in these conditions the energy consumption/generation and the purchase/sale price of energy increase/decrease compared to the case with RD = 0 according to Table 5.Table 6
Values of operation and economic indices in the different cases for various values of robustness degree.
RD00.1CaseIIIIIIIVVVIIIIIIIIVVVICost1119572571036286395376−1976.31234181211096492565376−139.8EL (MWh)7.2335.8926.9226.4625.1314.0487.9516.4457.6117.0235.5844.465MPF0.7990.8730.8350.7560.9080.9310.7990.8610.8300.7510.9020.926MVD (p.u)0.0920.0670.0840.0710.0640.0590.1080.0720.0910.0790.0700.064MOV (p.u)00.021000.0240.02600.018000.0220.023Figure 8
Curve of cost in ΔF in different values of RD.


## 4.2.1. Achieving a Suitable Solver

This section presents the convergence results of the proposed scheme for different values of robustness degree (RD orπ) obtained by mathematical solvers such as BARON, BONMIN, DISOPT, KNITRO, and OQNLP [66], and evolutionary algorithms such as teaching-learning-based optimization (TLBO) [70], differential equations (DE) [71], krill herd optimization (KHO) [72], gray wolf optimization (GWO) [73], and antlion optimization (ALO) [74]. To solve the proposed problem by mathematical algorithms, the problem model is coded in the GAMS optimization software, and each of the mentioned algorithms in this software has a toolbox [66]. To solve this problem using evolutionary algorithms, the proposed scheme along with the solution process based on this algorithm are coded in MATLAB software. First off, the algorithm (it can be any of the mentioned algorithms) generates N (population size) random values for decision-making variables including PMT, QMT, QRES, PCH, PDCH, x, QESS, PCHP, QCHP, PDR, HDR, σ, and u˜ in accordance with their operation limits given in equations (11) and (12), [0, SRESu], (15) and (16), {0, 1}, [0, SESSu], (23), (24), (27), [0, 1], and (40). Next, the values of dependent variables such as Cost, PG, QG, PL, QL, V, δ, PRES, E, HBO, and HCHP are calculated from equations (2)–(6), (13), (17), (18), (21), and (22). The values of PRES, E, HBO, and HCHP are determined from constraints (13), (17), (18), (21), and (22), while the values of other dependent variables are calculated from power flow constraints given in (2)–(6). In this paper, the forward-backward method is employed to solve the power flow problem [75]. In the following, the fitness function (FF) is calculated for N values of decision-making variables. To estimate the operation limits, (7)–(10), capacity limit of RESs, (14), the limit of storable energy in ESSs and their charger capacity, (19) and (20), thermal power limit of CHP and boiler, (25) and (26), the DRP limit, (27), and limit of changes of the objective function, (42), the FF in this paper will aim to maximize the difference between the main objective function, (41), and the sum of penalty functions (PeF), i.e., FF = max σ − PeF [14]. The PeF for limits a ≤ b and a = b is given as μ.max (0, a − b) and υ. (b − a), where μ ≥ 0 and υ ∈ (−∞, +∞) represents the Lagrangian multipliers, the values of which are determined like decision-making variables using the evolutionary algorithm [14]. In the following, the steps of updating random variables and Lagrangian multipliers using this algorithm based on their operating limits and the best value of FF in the previous step are implemented. The updating process differs depending on the type of evolutionary algorithm. This process for the TLBO, DE, KHO, GWO, and ALO algorithms can be found in [70–74], respectively. The current study assumes that the convergence point can be accessible after the maximum iterations (Imax). N and Imax are considered 80 and 4000 for evolutionary algorithms, and their other setting parameters were chosen from [70–74].The results of this section are presented in Table3. Based on this table, it can be seen that among the mathematical algorithms, BARON and KNITRO fail to achieve a feasible solution. Among BONMIN, DISOPT, and OQNLP, the BONMIN solver succeeds in obtaining the minimum cost and maximum rc compared to the other two solvers. It has achieved this point in the minimum convergence iterations (CI) and convergence time (CT) in comparison to other mathematical methods. In evolutionary algorithms, ALO has been able to find the most optimal solution (minimum cost and maximum rc) at the minimum CI and CT compared to TLBO, DE, KHO, and GWO. Furthermore, by comparing mathematical solution methods and evolutionary algorithms, according to Table 3, it is observed that ALO has more favorable conditions than BONMIN in terms of response speed and optimal point. This is because it has the minimum cost (maximum rc) compared to BONMIN. The CI of ALO is higher than BONMIN, but it has reached the optimal solution with a much lower CT than BONMIN. Therefore, in this paper, the ALO algorithm is a more suitable solver than the other solvers mentioned.Table 3
Convergence situation achieved by different solvers for different values of robustness degree.
π00.1SolverCost ($)rcCICT (min)Cost ($)rcCICT (min)BARONInfeasibleBONMIN−1892.1030739.7−113.50.10933842.8DISOPT−1611.4038451.4−78.20.10542255.2KNITROInfeasibleOQNLP−1537.5045767.2−75.70.10249772.5TLBO−1761.80242215.1−87.30.114267517.9DE−1756.30254915.6−84.50.111283118.4KHO−1812.50206112.3−99.70.124219714.1GWO−1838.60193811.5−104.20.127209513.2ALO−1976.3016729.2−139.80.134175610.1In Table4, the convergence status of the ALO-based problem using different uncertainty modeling methods for different values is presented. In this table, it can be seen that in the uncertainty modeling method based on scenario-based stochastic programming (SBSP) [15], a large number of scenarios are needed to find a reliable optimal solution. More than 80 scenarios show that the optimal solution has one solution (cost is constant), and the computational time and convergence iteration change. Therefore, it can be stated that 80 scenarios are needed to access a reliable solution in this method. The UT method [13] requires a total of 15 scenarios because, in (30), the number of uncertainty parameters is 7 and the total number of scenarios required in UT corresponding to 2n + 1 will be equal to 15 (2 × 7 + 1) scenarios. Therefore, with this number of scenarios, it has found an optimal solution such as SBSP in 80 scenarios, but its computational time is lower than SBSP. In the following, robust methods such as ARO [50] and IGDT have only one scenario. Of course, for π = 0, they consider deterministic modeling for uncertainties. Therefore, the value of the objective function in these methods is less than SBSP and UT, because the cost of modeling uncertainties is not obtained in these conditions. They also have one scenario, so their computational time is less than that of SBSP and UT. According to [50], the ARO formulation for a problem will complicate the problem. This makes the computation time longer than IGDT. Then, by increasing π, it is seen that the highest rc is obtained in IGDT. In ARO, the maximum radial uncertainty is obtained by the trial and error method and by changing the level of uncertainty. Hence, its computational time is longer than IGDT. Note that the details of SBSP, UT, and ARO are given in [13, 15], and [50], respectively. Based on equation (1), the term cost is equal to the difference between the sum of the operating cost of the network and nonrenewable sources and the sum of the revenues that power sources, storage, and responsive loads earn by selling electrical and thermal energy. With optimal management of sources, storage, and responsive loads (with their performance curves given in subsection 4.2.3), in the proposed scheme, it is inferred that the revenue of the mentioned elements overcomes their operating costs and that of the network (i.e., the revenue is greater than the operating cost). In these conditions, the value of cost in Tables 3 and 4 has become negative.Table 4
Convergence situation based on ALO in the different uncertainty modeling method for various values of robustness degree.
π00.1MethodCost ($)rcCICT (min)Cost ($)rcCICT (min)SBSP40 scenario−1669.8—263224.2There is no this condition60 scenarios−1663.3—286727.580 scenario−1657.1—307130.182 scenario−1657.1—307930.5UT (15 scenarios)−1657.6—243117.6ARO (1 scenario)−1976.3—182910.9−143.20.129216238.1 due to including the trial and error methodIGDT (1 scenario)−1976.3—16729.2−139.80.134175610.1


## 4.2.2. Investigation of Uncertainty Parameters in the Worst-Case Scenario

Table5 reports the true value of uncertainty variables for different values of RD (π). According to this table, it is observed that increasing the degree of robustness to 0.4 increases the maximum radius of uncertainty (maximum forecast error or rc) to 0.345. Also, in the conditions of increasing RD with respect to 0.4, the value of rc is always constant. This indicates the fact that the implementation of the proposed scheme on the data in subsection 4.1 can cover the prediction error of 34.5% resulting from the uncertainty parameters including load, energy price, and renewable power. If the prediction error (radius of uncertainty) is more than 34.5%, the problem has no feasible solution. Moreover, it is observed that with increasing the maximum radius of uncertainty, the amount of load and price of energy purchased from the upstream network in the worst-case scenario have an upward trend, but the power generation by renewables and energy selling price have a downward trend in this scenario. So, in general, it can be stated that in the worst-case scenario, compared to the scenario corresponding to the deterministic model (π = 0), the amount of energy consumed by the network (energy purchase price) increases while the energy produced in the network (energy sales price) decreases.Table 5
True value of uncertainty parameters at the worst-case scenario for different values of robustness degree.
π00.10.20.30.40.5rc00.1340.2210.3020.3450.345∑b,hP˜b,hD (MW)429.9487.5524.9559.7578.2578.2∑b,hQ˜b,hD (MVAr)322.6365.8393.9420433.9433.9∑b,hH˜b,hD (MW)104.3118.3127.3135.8140.3140.3∑b,hP˜b,hRESu (MW)79.969.262.255.852.352.3∑hλ˜h ($/MWh)611.6693.5746.8796.3822.6822.6∑hρ˜h ($/MWh)550.4476.6428.8384.2360.5360.5∑hγ˜h ($/MWh)616533.5479.9430403.5403.5


## 4.2.3. Investigation of the Operation and Economic Status of Distributed Energy Resources (DERs)

The daily curves of active and reactive power of sources, storage devices, and electrically responsive loads for different values of RD are shown in Figures4 and 5, respectively. Based on Figure 4(a), the data in subsections 4.1 and [1], corresponding to RESs, it can be seen that the active injection power of WSs and PVs at all operating hours is equal to the maximum active power they produce in accordance with the weather conditions. In addition, by increasing the RD, the amount of active power injection by RESs into the MG decreases in the worst-case scenario due to reducing P˜RESu in this condition based on Table 5. Figure 4(b) shows the daily power curves of CHPs and MTs. According to this figure, these DGs inject the maximum active power they can produce into the network during all operating hours. This is because (1) the energy purchase price (λ) is higher than the fuel price of DGs in some hours, hence the MG operator tends to utilize these DGs to supply energy consumption, and (2) these DGs sell energy to the MG at a price ρ, which is about 111% (100%/90%) of λ based on subsection 4.1. Therefore, to minimize the Cost function in equation (1), DGs are expected to inject as much energy as possible into the network, so they will always operate at the point corresponding to the maximum active output power. Even in the worst-case scenario (RD = 0.2), where λ (ρ) has increased (decreased) by about 22% compared to RD = 0 according to Table 5, the performance results of DGs are the same as deterministic studies (RD = 0). Regarding the performance of DRPs and ESSs in Figures 4(c) and 4(d), they perform charging operations during 1 : 00–16 : 00 and 23 : 00–00 : 00 so that during 1 : 00–7:00, high energy is received by DRP and ESSs from MG. They also operate in discharge mode from 17 : 00 to 22 : 00 and inject active power into the MG. The performance of these elements corresponds to the electricity price. According to [5], energy prices in the hours 1 : 00–7:00, 8 : 00–16 : 00, and 23 : 00–00 : 00, and 17 : 00–22 : 00 are equal to 17.6 $/MWh, 26.4 $/MWh, and 33 $/MWh. Hence, to minimize the Cost function, they will operate in charging (discharging) mode at times corresponding to the cheap (expensive) energy price. As another point, note that since the active power generation by RESs decreases in the worst-case scenario according to Figure 4(a), DRPs and ESSs inject a higher active power in this scenario during peak hours (17 : 00–22 : 00) into MG than in the case with RD = 0. They also increase their power consumption during other hours to be able to inject high energy into the network during other hours. Finally, note that CHP, MT, ESSs, and DRP are flexibility sources for RES. As RESs generate active power at all operating hours, it is expected that MG flexibility control is required at all hours. Therefore, flexible sources according to Figure 4 are not switched off at any time of operation.Figure 4
Daily active power curves of (a) RESs, (b) non-RESs, (c) DRPs, and (d) ESSs at different values of robustness degree.
(a)(b)(c)(d)Figure 5
Daily reactive power curve of (a) PVs, (b) WSs, (c) non-RESs, and (d) ESSs in the different values of robustness degree.
(a)(b)(c)(d)The daily reactive power curves of DGs and ESSs for different RD values are shown in Figure5. According to this figure, DGs inject the maximum active power into the network during most of the operating hours. For example, according to the data in subsection 4.1, the maximum capacity (apparent power) of PVs is 1.5 MVA (3 PV × 0.5 MVA for PV capacity), and PVs do not produce active power during 1 : 00–5:00 and 1900–00 : 00 according to Figure 4(a). Therefore, they inject reactive power of 1.5 MVAr into the network during these hours. However, in other hours, when their active power generation increases, their reactive power generation decreases. Moreover, with increasing RD, the active power generation by PVs decreases according to Figure 4(a), so they inject more reactive power into MG than in the case with RD = 0. These conditions hold for WSs according to Figure 5(b). But CHPs (MTs) injected the maximum reactive power that can be produced, i.e., 7 × 0.4 = 2.8 MVAr (4 × 0.25 = 1 MVAr) into the network at all hours and for all values of RD. In the following, comparing the load profile presented in [17] and Figure 5(d), it is observed that the daily reactive power curve of ESSs is the same as the daily reactive load profile. Since the amount of reactive power consumed by MG increases if RD is increased (Table 5), ESSs produce high reactive power. Note that this performance of sources and ESSs in the MG reactive power supply process corresponds to maximizing its revenue in the reactive power market based on equation (1).Figure6 depicts the thermal capacity balance curve of the sources and the response loads of the thermal block for different values of RD. According to Figure 6(a), thermal DRP reduces thermal energy consumption during peak thermal hours, i.e., 5 : 00–15 : 00, but at other times, it increases the energy consumption of thermal consumers. According to [5], the price of thermal energy in the hours 5 : 00–15 : 00 is 30 $/MWh and 22 $/MWh at other hours. Hence, the performance of thermal DRP in the figure corresponds to the minimization of the Cost function in (1). In addition, since the active power generation by CHP was constant at all operating hours according to Figure 3(b), the amount of heat generated by CHP, according to (22), will be constant for time changes and will have a daily curve as shown in Figure 6(b). Finally, according to (21), the boiler heat output will be equal to the difference between the thermal load and the total heat output of CHP and thermal DRP, so its daily profile will be as shown in Figure 6(b). Then, with increasing RD, the daily heat power curve of the boiler and DRP shifts upwards relative to the case RD = 0, because in the daily curve, the thermal load shifts upwards in the worst-case scenario according to Figure 6(a). Nonetheless, because the thermal power of CHP depends on its active power, and the active power of CHP, according to Figure 4(b), always has a constant value for different values of RD, the thermal power of CHP will not change with changes in RD. Finally, the economic indices curves of MG and DERs, including sources, storage devices, and responsive loads, are plotted in terms of RD in Figure 7. According to this figure, increasing the RD to 0.4 increases the operating costs of MGs and NRDGs, while the revenue of DERs decreases under these conditions. This leads to an increase in the net network cost (Cost function), i.e., the difference between the total costs of MG and NRDGs and the revenue of DERs in the energy and reactive power markets under the mentioned conditions. Nonetheless, the values of the economic indices mentioned for RD > 0.4 are the same as for RD = 0.4. All these results obtained in Figure 6 are proportional to the change in the status of energy consumption/generation and the purchase/sale price of energy in the worst-case scenario in Table 5 or subsection 4.2.2.Figure 6
Daily heat power curve of (a) block load without/with DRP and (b) sources in the different values of robustness degree.
(a)(b)Figure 7
Curve of economic indices in robustness degree.


## 4.2.4. Evaluation of the Economic and Technical Status of MG

In this section, six case studies are analyzed to examine the capabilities of the proposed scheme, as follows:(i)
Case I: Power flow analysis (MG without sources, storage devices, DRP, and thermal block)(ii)
Case II: Proposed scheme considering only RESs(iii)
Case III: Proposed scheme considering only ESSs(iv)
Case IV: Proposed scheme considering only DRPs(v)
Case V: Proposed scheme considering NRDGs (MTs and thermal blocks)(vi)
Case VI: Proposed scheme considering all proposed DERsThe results of this section are shown in Table6, which reports the values of cost, energy loss (EL), minimum power factor (MPF) of the distribution substation, maximum voltage drop (MVD), and maximum overvoltage (MOV) in cases I to VI for different values of RD. According to this table, the presence of RESs alone in the proposed scheme (II) or the single presence of thermal blocks and MT in MG (V) leads to a significant reduction in energy costs, energy losses, and the voltage drop of MG compared to Case I. Also, MPF increases significantly; in Case V it has increased by more than 0.9. These conditions correspond to the occurrence of MOV around 0.02 p.u., which is less than its allowable limit, i.e., 0.05 (1.05–1). The performance of ESSs in Case III has a favorable effect on improving the status of EL, MPF, and MVD and a low impact on cost reduction without causing overvoltage compared to power flow studies, but the percentage of improvement of these parameters, in this case, is less than cases II and V because in cases II and V, high active and reactive power are simultaneously injected into the network by DGs based on Figures 4 and 5. But, in Case III, ESSs inject active power into the network only during peak hours and consume active power during other hours. The ESSs also failed to generate as much reactive power as the DGs; this is because of their low charger capacity and low number compared to the DGs in the 119-bus network based on subsection 4.1. In the case of DRP, it has been able to enhance cost, EL, and MVD in Case IV compared to Case I without causing overvoltage but has led to MPF degradation. Note that the percentage of improvement of cost, EL, and MVD in this case study is lower than in cases II and V but higher than in Case III. In addition, in Case IV, DRP only plays a role in controlling the active power and does not change the amount of reactive power. Therefore, according to Figure 4(c), during peak hours when the amount of active power consumption (active power of the distribution substation) is reduced by DRP with no change in reactive power status compared to Case I, based on equation (9), the power factor of the distribution substation will be reduced compared to Case I. Finally, the utilization of all DERs in the proposed scheme (VI) has enhanced all the proposed indices compared to Case I. In the worst-case scenario (RD = 0.1), the Cost function in Case VI is reduced by about 101% compared to Case I (12341-(-139.8))/12341). This value (percentage of improvement) for EL, MPF, and MVD in the worst-case scenario with a MOV of 0.023 p.u. is about 43.8%, 15.9%, and 40.7%, respectively. Note that according to Table 6 for an increase in RD, all the mentioned indices except for MOV and MPF increase compared to RD = 0, but MOV and MPF decrease. The reason is that, according to Table 5 in the worst-case scenario compared to the deterministic model, energy consumption (generation) increases (decreases). Finally, the cost curve in terms of flexibility tolerance (∆F) for different RD values is plotted in Figure 8. Based on this figure, it can be seen that with increasing ∆F, Cost decreases because, in these conditions, the operating cost of resources, storage devices, and responsive loads is reduced. For example, if ∆F increases, then the importance of network flexibility decreases. Therefore, electric DRPs and ESSs will be switched off at 8 : 00–16 : 00 and 23 : 00–00 : 00 and will not receive power from the network. This results in a reduction in their operating costs. As the RD increases, the curve shifts upwards, because in these conditions the energy consumption/generation and the purchase/sale price of energy increase/decrease compared to the case with RD = 0 according to Table 5.Table 6
Values of operation and economic indices in the different cases for various values of robustness degree.
RD00.1CaseIIIIIIIVVVIIIIIIIIVVVICost1119572571036286395376−1976.31234181211096492565376−139.8EL (MWh)7.2335.8926.9226.4625.1314.0487.9516.4457.6117.0235.5844.465MPF0.7990.8730.8350.7560.9080.9310.7990.8610.8300.7510.9020.926MVD (p.u)0.0920.0670.0840.0710.0640.0590.1080.0720.0910.0790.0700.064MOV (p.u)00.021000.0240.02600.018000.0220.023Figure 8
Curve of cost in ΔF in different values of RD.


## 5. Conclusion

This paper presents the electrical and thermal power management in microgrids with thermal blocks and distributed energy sources in accordance with the conditions of the day-ahead energy and reactive power markets. The thermal block is equipped with a CHP, a boiler, and a DRP to supply its energy consumption. Subsequently, the scheme aims for minimizing the difference between the total operating cost of the MG and NRDGs and the revenue of DERs in the mentioned markets. It was also bound by the constraints of optimal AC power flow and the operating model of the thermal block and DERs. In this paper, robust IGDT-based programming based on the network flexibility model was adopted to model uncertainties of load, market price, and renewable power. In the end, based on the obtained numerical results, it was observed that the ALO algorithm has a good ability to achieve the optimal solution in the shortest computational time compared to mathematical solution methods and some other evolutionary algorithms. Also, by achieving the optimal performance curve in active, reactive, and thermal power of sources, storage devices, and responsive loads, the proposed scheme was able to enhance energy cost, energy losses, the minimum power factor of the distribution substation, and the maximum voltage drop of the network in the worst-case scenario compared to power flow studies by about 101%, 43.8%, 15.9%, and 40.7%, respectively, by creating a maximum overvoltage of 0.023 p.u. Furthermore, to achieve high flexibility in the network, it is necessary to increase the operating cost of flexibility resources, which increases the energy cost. It should be noted that the proposed design based on the IGDT model had a lower computing time than other uncertainty modelings due to its low volume and simplicity of the process. In addition, it can calculate the maximum radial uncertainty or the maximum prediction error in the worst-case scenario. Nonetheless, this issue is not established in stochastic modeling of uncertainties and other methods of a robust model. Furthermore, the increase in robustness degree in IGDT increases the cost of purchasing energy and the exploitation of renewable resources, but the income from the sale of energy sources and storage decreases. This is caused by the increase (or decrease) in energy demand (resource capacity in power generation) in the worst-case scenario.The proposed scheme corresponds to a balanced MG, but an imbalance between the three phases may arise due to the imbalance. Therefore, there is a need for unbalanced MG modeling in the proposed design, which will be explored in future work. The proposed scheme is also able to improve the reliability of the network. Since the sources, storage units, and thermal blocks are distributed in the MG consumption areas, they can feed a specific percentage of consumers in the event of a fault and reduce network outages. This depends on the modeling of the proposed design based on reliability, which is considered for future work. By distributing the mentioned elements in different areas of the MG, it is predicted that the proposed scheme can achieve a high safety margin for voltage. This can be achieved by modeling the voltage security index in the proposed scheme, which will also be proposed in future work.



---

*Source: 6002695-2022-12-13.xml*