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

,


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][2][3][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 profles, resulting in improved operation indices such as reduced energy losses and voltage drops, and reducing the operating cost of the network [7]. Terefore, it is predicted that the number of these elements on the consumption side will increase. Te 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). Terefore, 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.

Literature Review.
Various studies have been carried out in the feld 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. Te 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 fexibility point of view. An electric spring is a DC/AC converter and can control its fowing 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]. Te optimization problem includes an objective function to minimize several parameters such as operating cost and voltage variations and maximize the fexibility of the system. Te 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]. Tis 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. Te UTmethod 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. Te approach is mathematically expressed in the form of an optimization problem with four objectives subject to optimal power fow (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. Te 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. Te 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. Te EMS consists of two layers, and the problem has two stages. Te 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 pricebased 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. Te 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. Te 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. Te paper investigates two methods of utilizing solar energy. Tese 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 fnd 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]. Te heating network is frst 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. Te 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. Te strategy targets CCHP MGs that contain battery charging terminals. Te 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]. Te 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 fxed and that of heating and cooling energy is changing. Te 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]. Te 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 benefts of adopting specifc features of energy types were considered. Te paper used a volt/var control to provide an optimization of active power fow and reactive power fow so that voltage security was ensured. It also presented modeling of battery degradation and thermal network to reach efective 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 fexible demand are properly dispatched. Te proposed method attempted to model coupling limitations of electrical and heating networks, the dynamic features of the heating network, and power fow limits in a comprehensive way. Table 1 reports the majority of the studies presented within the last decade.

Research Gaps.
According to the literature and Table 1, there are the following major research gaps in the feld 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. Te CHP also generates both electrical and thermal energies at its output. Terefore, to better model MG operation, it is necessary to formulate the thermal model of sources and consumers. Tis 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 fnancial beneft 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. Tus, they can also beneft 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 signifcant 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 defned 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.
International Transactions on Electrical Energy Systems

Solutions and Contributions.
In this paper, to deal with the frst and second research gaps, power management in an MG with a thermal block, ESSs, DRPs, and DGs is proposed. Tis 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 diference 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 fow (AC-OPF) equations, and the operating model of the thermal block, DGs, ESSs, and DRPs. Tis scheme has uncertainties in energy price, renewable power, and load consumption. Terefore, to address the third and fourth research gaps, they are modeled based on IGDT-based robust scheduling. Also, to consider the adverse efects of uncertainties on the proposed scheme and the risk model, the IGDT model in this paper is based on a risk-averse strategy. Tis 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 fnancial benefts 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.

Paper
Organization. Te paper is organized as follows: In Section 2, the fnal model of the proposed scheme is stated. Ten, in Section 3, uncertainty modeling is presented. Te numerical results obtained from diferent case studies are evaluated in Section 4. Finally, the conclusions are summarized in Section 5.

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 Figure 1 is presented in the framework of deterministic optimization. Tis scheme is expressed as an optimization formulation. Tis formulation includes objective functions [29][30][31][32][33] and constraints [34][35][36][37]. Tis scheme is an attempt to minimize the diference 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]. Te smart grid needs telecommunications [39][40][41][42] and smart devices [43][44][45][46]. Te deterministic model of the proposed scheme can be given as follows: Subject to:  International Transactions on Electrical Energy Systems ������������� � 2.1. Objective Function. Equation (1) formulates the objective function of the proposed scheme, which is equal to minimizing the diference 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 frst 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 frst to third terms refer to the operating (fuel) costs of NRDGs of MT, CHP, and boiler types, respectively [47]. Te 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. Terefore, the energy purchase price of MG (λ) will be diferent from the energy sale price in the retail market (ρ and c). Also, in this section, it is assumed that the price of reactive power is a multiplication of the price of electrical energy, i.e., K Q × ρ. 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).

MG Operation Constraints.
Tese constraints are modeled as AC-OPF in equations (2)-(10) [48][49][50][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 PF l (generally has a value of 0.9) in all operating hours. Te 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 P G and Q G in all buses will be zero. Te proposed scheme can be used for diferent types of distribution networks in accordance with the problem model given in equations (1)- (28). Tis is because the optimal power fow model, (2)-(10), can be applied to all distribution networks. Index ref should be defned as the buses with the distribution substation just to assure this point. (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. Te 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]. Terefore, 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.

Te Operation Model of Electrical Energy Storage Devices.
Te 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 frst 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. Tis 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]. Terefore, 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.

Te Operation Model of the Termal Block.
Te 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]. Termal 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. Te limit of thermal power that can be produced by CHP is also considered in equation (25). Tere is no restriction on the type of CHP in this article because the formulations (22)-(25) can be used International Transactions on Electrical Energy Systems for diferent types of CHP [5]. Tere 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]. Tis type of DRP fts 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). Terefore, 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. Terefore, to meet these conditions, constraint (28) is added to the formulation of the proposed scheme.
Since the frst two types of modeling require a signifcant number of scenarios and identifcation of the PDF for uncertainties [61], the computational time is high [9][10][11][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 efect 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 defned as a variable. Also, since this variable is to determine the worst-case scenario this variable is a decision variable in IGDT.

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. Te approach is independent of PDFs. Te problem modeling is described as follows: where, u denotes input uncertainty variables, and Ξ represents the set of uncertainties. Y shows decision-making variables. Te 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]: where, u represents the uncertain variables matrix, and σ is the diference between the maximum deviation of uncertainty variables and its estimated or nominal value (u). Te variable σ is known as the uncertainty level and its value is assumed constant for all uncertainties in this paper. Terefore, it has a scalar value. Te 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. Tis mathematical expression is called the based condition (BC), where its objective function is given by f � min Y F(Y, u). Te maximum radius of uncertainty (r c ) is found as follows [63]: Constraints (32) and (33), where Δ c denotes the robustness value of the problem described in equations (31)- (34). Te 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), f b is the value of the BC objective function and π ∈ [0 · 1] shows the degree of robustness when f b deteriorates. Consequently, in equations (36)- (40), π is an input parameter and f b shows an input parameter found by equations (31)-(33) taking into account only u.

International Transactions on Electrical Energy Systems
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. Te 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. Tus, 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 diferent uncertainty radii for diferent uncertainties. Under such conditions, an index such as i is added to σ, which includes uncertainty variables of equation (30). Ten, equation (36) is changed as r c � max Y,σ i σ i , but generally, concerning the IGDT, the uncertainty level is considered fxed for all uncertainties. Tis can be found in [19].

Robust
Formulation of the Proposed Scheme. As per subsection 3.2, the IGDT-based robust model proposed for the proposed problem, described by equations (1)- (28), can be given as follows: Subject to: Constraints (39) and (40), where Cost b 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. λ, ρ, c, P RES , as given in (37).
In the current study, the uncertainties are modeled using the IGDT method. Te model requires one scenario named "worst-case scenario." As a result, the IGDT method fnds 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 fexibility 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]. Terefore, there is a need to use fexibility resources such as DRP and storage in the network, which can compensate for fuctuations in RESs active power in real-time operation compared to DA operation and enhance the fexibility of the network [13]. To establish such conditions, from the point of view of mathematical modeling, the diference 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. Tis is presented in (44), which should be added to the problem described in (41)- (43). If the fexibility tolerance (∆F) tends to zero, the fexibility will be 100%, and the imbalance between the mentioned operations in the network with RESs will be eliminated. Finally, the fowchart of the problem solution is shown in Figure 2. Figure 3 depicts a 119-bus test MG used to validate the efcacy of the suggested method [64]. Te 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]. Te 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]. Te 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. Te maximum/minimum values of active, reactive, and thermal power of the CHP unit are 800/ International Transactions on Electrical Energy Systems 9 0 kW, 400/−400 kVAr, and 500/0 kW. Te specifed values for efciencies η T , η H , η L are assumed 0.4, 0.4, and 0.08 [5]. Te maximum/minimum heat power output from the boiler is 500/0 kW. Coefcients α, β, and χ of the MT, CHP, and boiler are (100, 21, and 0.002), (93, 20, and 0.0017), and (112, 22, and 0.002). Te battery-type ESS has an efciency of 0.9, charge/discharge rate is 150 kW, charger capacity is 200 kVA, and its minimum/maximum permissible energy is 50/500 kWh. Te location of the battery is given in Table 2. Te 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. Te parameter ξ in the DRP is assumed to be 0.3 [7]. Te location of the PV, WS, and MT are given in Table 2. Te 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]. Te authors in Ref. [5] provide information about the hourly prices of electrical and heat energies. Te parameter K Q is assumed 0.08 [65], and λ � 0.9 × ρ, and ΔF is assumed to be 0.05 p.u. to achieve high fexibility 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.

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

Achieving a Suitable Solver.
Tis section presents the convergence results of the proposed scheme for diferent 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], diferential 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 of, the algorithm (it can be any of the mentioned algorithms) generates N (population size) random values for decisionmaking variables including P MT , Q MT , Q RES , P CH , P DCH , x, Q ESS , P CHP , Q CHP , P DR , H DR , σ, and u in accordance with their operation limits given in equations (11) and (12) (22), while the values of other dependent variables are calculated from power fow constraints given in (2)-(6). In this paper, the forward-backward method is employed to solve the power fow problem [75]. In the following, the ftness 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 diference between the main objective function, (41), and the sum of penalty functions (PeF), i.e., FF � max σ − PeF [14]. Te 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. Te updating process difers depending on the type of evolutionary algorithm. Tis process for the TLBO, DE, KHO, GWO, and ALO algorithms can be found in [70][71][72][73][74], respectively. Te current study assumes that the convergence point can be accessible after the maximum iterations (I max ). N and I max are considered 80 and 4000 for evolutionary algorithms, and their other setting parameters were chosen from [70][71][72][73][74].
Te results of this section are presented in Table 3. 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 fnd the most optimal solution (minimum cost and maximum r c ) 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. Tis is because it has the minimum cost (maximum r c ) compared to BONMIN. Te CI of ALO is higher than BONMIN, but it has reached the optimal solution with a much lower CT than BONMIN. Terefore, in this paper, the ALO algorithm is a more suitable solver than the other solvers mentioned. In Table 4, the convergence status of the ALO-based problem using diferent uncertainty modeling methods for diferent 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 fnd 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. Terefore, it can be stated that 80 scenarios are needed to access a reliable solution in this method. Te UT method [13] requires a total of 15 scenarios because, in (30), the number of uncertainty parameters is 7 and the total 19 Figure 3: 119-Bus MG [64]. number of scenarios required in UT corresponding to 2n + 1 will be equal to 15 (2 × 7 + 1) scenarios. Terefore, 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. Terefore, 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. Tey 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. Tis makes the computation time longer than IGDT. Ten, by increasing π, it is seen that the highest r c 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 diference 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 5 reports the true value of uncertainty variables for diferent 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 r c ) to 0.345. Also, in the conditions of increasing RD with respect to 0.4, the value of r c is always constant. Tis 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.

Investigation of the Operation and Economic Status of Distributed Energy Resources (DERs).
Te daily curves of active and reactive power of sources, storage devices, and electrically responsive loads for diferent values of RD are shown in Figures 4 and 5, respectively. Based on Figure 4(a), the data in subsections 4.1 and [1], corresponding to RESs, it   International Transactions on Electrical Energy Systems 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 fgure, these DGs inject the maximum active power they can produce into the network during all operating hours. Tis 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. Terefore, 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. Tey also operate in discharge mode from 17 : 00 to 22 : 00 and inject active power into the MG. Te 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. Tey 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 fexibility sources for RES. As RESs generate active power at all operating hours, it is expected that MG fexibility control is required at all hours. Terefore, fexible sources according to Figure 4 are not switched of at any time of operation.
Te daily reactive power curves of DGs and ESSs for diferent RD values are shown in Figure 5. According to this fgure, 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). Terefore, 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. Tese 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 profle 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 profle. 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). Figure 6 depicts the thermal capacity balance curve of the sources and the response loads of the thermal block for diferent 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 fgure 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 diference between the thermal load and the total heat output of CHP and thermal DRP, so its daily profle will be as shown in Figure 6(b). Ten, 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 diferent 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 fgure, increasing the RD to 0.4 increases the operating costs of MGs and NRDGs, while the revenue of DERs decreases under these conditions. Tis leads to an increase in the net network cost (Cost function), i.e., the diference 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  Te results of this section are shown in Table 6, 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 diferent 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 signifcant reduction in energy costs, energy losses, and the voltage drop of MG compared to Case I. Also, MPF increases signifcantly; in Case V it has increased by more than 0.9. Tese 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). Te performance of ESSs in Case III has a favorable efect on improving the status of EL, MPF, and MVD and a low impact on cost reduction without causing overvoltage compared to power fow 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. Te 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. Terefore, 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). Tis 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. Te 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 fexibility tolerance (∆F) for diferent RD values is plotted in Figure 8. Based on this fgure, 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 fexibility decreases. Terefore, electric DRPs and ESSs will be switched of at 8 : 00-16 : 00 and 23 : 00-00 : 00 and will not receive power from the network. Tis 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.

Conclusion
Tis 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. Te thermal block is equipped with a CHP, a boiler, and a DRP to supply its energy consumption. Subsequently, the scheme aims for minimizing the diference 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 fow and the operating model of the thermal block and DERs. In this paper, robust IGDTbased programming based on the network fexibility 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 fow 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 fexibility in the network, it is necessary to increase the operating cost of fexibility 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. Tis is caused by the increase (or decrease) in energy demand (resource capacity in power generation) in the worst-case scenario.
Te proposed scheme corresponds to a balanced MG, but an imbalance between the three phases may arise due to the imbalance. Terefore, there is a need for unbalanced MG modeling in the proposed design, which will be explored in future work. Te 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 specifc percentage of consumers in the event of a fault and reduce network outages. Tis depends on the modeling of the proposed design based on reliability, which is considered for future work. By distributing the mentioned elements in diferent areas of the MG, it is predicted that the proposed scheme can achieve a high safety margin for voltage. Tis can be achieved by modeling the voltage security index in the proposed scheme, which will also be proposed in future work. Heat demand (p.u.) K Q :

(1) Constants
Te ratio between the reactive power price and energy price P CHu , P DCHu : Charge and discharge efciency of ESS η T , η L , η H : Te efciency of turbine, loss efciency, and heat efciency in the CHP λ: Te price of purchasing electrical energy by the MG from the upstream market ($/MWh) ρ, c: Te prices of electrical and heat energy in the retail market ($/MWh) π: Degree of robustness ξ: Coefcient of the participation rate in the DRP ΔF: Flexibility tolerance (p.u.) Maximum radius of uncertainty u:

(2) Variables
Matrix of uncertainty variables V, δ: Voltage magnitude (p.u.) and voltage angle (rad) x: A binary variable corresponding to the charge/discharge mode of ESS σ: Uncertainty level Δ c : Robustness value of the problem of microgrid (MG) operation

Sets and indices b:
Bus index and set of buses h, OH: Index and set of operation hours k: An auxiliary index representing the bus ref: Slack bus Ξ: Uncertainty set.

Data Availability
Data will be available upon request from the corresponding author.

Conflicts of Interest
Te authors declare that there are no conficts of interest.