This is a self-archived – parallel published version of this article in the publication archive of the University of Vaasa. It might differ from the original. Optimal scheduling of CCHP-based resilient energy distribution system considering active microgrids' multi-carrier energy transactions Author(s): Armioun, Majid; Nazar, Mehrdad Setayesh; Shafie-khah, Miadreza; Siano, Pierluigi Title: Optimal scheduling of CCHP-based resilient energy distribution system considering active microgrids' multi-carrier energy transactions Year: 2023 Version: Accepted manuscript Copyright ©2023 Elsevier. This manuscript version is made available under the Creative Commons Attribution–NonCommercial–NoDerivatives 4.0 International (CC BY–NC–ND 4.0) license, https://creativecommons.org/licenses/by-nc-nd/4.0/ Please cite the original version: Armioun, M., Nazar, M. S., Shafie-khah, M. & Siano, P. (2023). Optimal scheduling of CCHP-based resilient energy distribution system considering active microgrids' multi-carrier energy transactions. Applied Energy 350, 121719. https://doi.org/10.1016/j.apenergy.2023.121719 Optimal scheduling of CCHP-based resilient energy distribution system considering active microgrids’ multi-carrier energy transactions Majid Armioun a, Mehrdad Setayesh Nazar a, Miadreza Shafie-khah b, Pierluigi Siano c,d a Faculty of Electrical Engineering, Shahid Beheshti University, Tehran, Iran b School of Technology and Innovations, University of Vaasa, 65200 Vaasa, Finland c Department of Management & Innovation Systems, University of Salerno, Fisciano, SA 84084, Italy d Department of Electrical and Electronic Engineering Science, University of Johannesburg, Johannesburg 2006, South Africa Abstract This paper introduces a two-stage two-level optimization method for optimal day-ahead and real-time scheduling of multicarrier en- ergy distribution systems and microgrids. The model considers the incentive-based and pricebased demand response programs to encourage microgrids to transact electrical, heating, and cooling energy carriers with the energy distribution system, which is named hereafter as the energy system. Further, the model formulates the resilient operation of the energy system considering the en- ergy transactions with the electrical, heating, and cooling markets. The main contribution of this paper is the integration of de- mand response procedures of microgrids in energy transactions with the energy system considering the switching of electrical switches and heating and cooling control valves. The optimization process is another contribution of this paper that is decomposed into two stages that consist of day-ahead and real-time horizons. The first stage is also decomposed into two levels that deter- mine the optimal scheduling of the energy system and microgrids in day-ahead markets. The second stage is comprised of two levels that commit the energy system and microgrids resources. A resiliency index is proposed to assess the resiliency of the energy sys- tem in shock conditions. The proposed method was simulated for the 123-bus test system. Different types of microgrids, incentive- based and price-based demand response processes were considered. Simulation results confirmed that the proposed method can re- duce the costs of residential, industrial, and commercial microgrids by about 4.47%, 3.88%, and 5.47% concerning only the real- time pricing process. Further, the model can increase the aggregated benefits of the energy system in the day-ahead and real-time markets by about 0.608 Million Monetary Units (MMUs) and 1.10 MMUs, respectively. 1. Introduction 1.1. Literature review The optimal operation of multi-carrier energy distribution systems is an important issue considering the energy transactions of these systems with the multi-carrier energy markets and microgrids. The energy sys- tem can be equipped with multi-carrier Distributed Energy Resources (DERs) and can utilize different demand response programs to encour- age the Microgrids (MGs) to change their load and energy transaction profiles. The multi-carrier energy transactions of microgrids mailto:psiano@unisa.it www.sciencedirect.com/science/journal/03062619 https://www.elsevier.com/locate/apenergy https://doi.org/10.1016/j.apenergy.2023.121719 https://doi.org/10.1016/j.apenergy.2023.121719 https://doi.org/10.1016/j.apenergy.2023.121719 Nomenclature Abbreviation ARIMA Autoregressive Integrated Moving Average CCHP Combined Cool, Heat and Power. CCP Chance-Constrained Programming CVaR Conditional Value-at-Risk DA Day-Ahead DER Distributed Energy Resources DRP Demand Response Program DSO Distribution System Operator ESS Electrical Energy Storage GA Genetic Algorithm IBDR Incentive-Based Demand Response IBEDRP Incentive-Based Electrical Demand Response Program IBHDRP Incentive-Based Heating Demand Response Program IBCDRP Incentive-Based Cooling Demand Response Program MG Micro Grid MILP Mixed Integer Linear Programming MINLP Mixed Integer Non-Linear Programming MIQP Mixed Integer Quadratic Programming MPEC Mathematical Programming with Equilibrium Constraints MCESRI Multi-Carrier Energy System Resiliency Index PBEDRP Price-Based Electrical Demand Response Program PBHDRP Price-Based Heating Demand Response Program PBCDRP Price-Based Cooling Demand Response Program PHEV Plug-in Hybrid Electrical Vehicle PSO Particle Swarm Optimization Sets J Set of MGs that participated in IBEDRP J’ Set of MGs that participated in IBHDRP J" Set of MGs that participated in IBCDRP G Set of MGs that participated in PBEDRP G’ Set of MGs that participated in PBHDRP G" Set of MGs that participated in PBCDRP Parameters APV Area of photovoltaic array η Photovoltaic array conversion efficiency I Solar irradiation t0 Outside air temperature γDS(t) Electrical generation costs of energy system SUi DS(t) Startup costs of the energy system γ’DS(t) Heating energy generation costs of energy system γ’’DS(t) Cooling energy generation costs of the energy system EENSCDS(t) DSO’s electrical energy not supplied costs HENSCDS(t) DSO’s heating energy not supplied costs CENSCDS(t) DSO’s cooling energy not supplied costs ENPDG CO2(t), ENPDG SO2(t), ENPDG NOx(t) Volume of CO2, SO2, and NOx pollution of the distributed generation unit ENPMarket CO2 Market(t), ENPMarket SO2 Market(t), ENPMarket NOx Market(t) Volume of CO2, SO2, and NOx pollution of electrical wholesale market generation unit PL DS(t) Active power of DSO’s loads HL DS(t) Heating power of DSO’s loads RL DS(t) Cooling power of DSO’s loads γMG(t) Cost of energy generation of MG γ’’MG(t) Cost of cooling energy generation of MG γ’MG(t) Cost of heating energy generation of MG γMG(t) Cost of electrical energy generation of MG EENSCS MG(t) Electrical energy not supplied costs of MG HENSCS MG(t) Heating energy not supplied costs of MG CENSCS MG(t) Cooling energy not supplied costs of MG ENPDG CO2(t), ENPDG SO2(t), ENPDG NOx(t) Volume of CO2, SO2, and NOx pollution of MG’s distributed generation unit HL MG(t) MG’s heating power of loads HACH MG (t) MG’s heating energy consumption of absorption chiller RL m,s MG MG’s cooling power of loads ECICDS(t) Electrical energy not supplied costs of energy system HCICDS(t) Heating energy not supplied costs of energy system CCICDS(t) Cooling energy not supplied costs of energy system PCritical,HCritical,RCritical Critical electrical loads, critical heating loads, and critical cooling loads ECICMG(t) MG’s electrical consumer interruption costs HCICMG(t) MG’s heating consumer interruption costs CCICMG(t) MG’s heating consumer interruption costs vci Cut-in wind velocity vr Nominal wind velocity vco Cut-out wind velocity vwind wind velocity T Number of operational scheduling hours of the energy system DSDAS Number of energy system day-ahead operating scenarios DSRTS Number of energy system real-time operating scenarios NDG Number of energy system DGs NHDER Number of energy system heating DERs NCDER Number of energy system’s cooling DERs NIDG Number of energy system intermittent DGs NEL Number of electrical load buses of the energy system NHL Number of heating load buses of the energy system NCL Number of cooling load buses of the energy system NRMG Number of residential MGs NCMG Number of commercial MGs NIMG Number of industrial MGs NEIBDR Number of electrical loads participating in IBEDRP NHIBDR Number of heating loads participating in IBHDRP NCIBDR Number of cooling loads participating in IBCDRP MGDAS Number of MG’s day-ahead operating scenarios MGRTS Number of MG’s real-time operating scenarios NDG’ Number of MGS’ DGs NHRMG Number of heating energy resources of MGs NCRMG Number of cooling energy resources of MGs NIDG’ Number of MGs’ intermittent DGs NEL’ Number of electrical load buses of MGs NHL’ Number of electrical load buses of MGs NCL’ Number of electrical load buses of MGs NRTSS Number of energy system real-time simulation steps NNCSW Number of energy system electric switches NNCHV Number of energy system heating energy control valve NNCCV Number of energy system cooling energy control valve NCEL Number of critical electrical load buses of the energy system NCHL Number of critical heating load buses of the energy system NCCL Number of critical cooling load buses of the energy system Variables AELRC(r, t) Actual electrical load reduction of residential MG SELRC(r, t) Submitted value of electrical load reduction of residential MG AELCC(r, t) Actual electrical load reduction of commercial MG SELCC(r, t) Submitted value of electrical load reduction of commercial MG AELIC(r, t) Actual electrical load reduction of industrial MG SELIC(r, t) Submitted value of electrical load reduction of industrial MG AHLRC(r, t) Actual heating load reduction of residential MG SHLRC(r, t) Submitted value of heating load reduction of residential MG AHLCC(r, t) Actual heating load reduction of commercial MG SHLCC(r, t) Submitted value of heating load reduction of commercial MG AHLIC(r, t) Actual heating load reduction of industrial MG SHLIC(r, t) Submitted value of heating load reduction of industrial MG ACLRC(r, t) Actual cooling load reduction of residential MG SCLRC(r, t) Submitted value of cooling load reduction of residential MG ACLCC(r, t) Actual cooling load reduction of commercial MG SCLCC(r, t) Submitted value of cooling load reduction of commercial MG ACLIC(r, t) Actual cooling load reduction of industrial MG SCLIC(r, t) Submitted value of cooling load reduction of industrial MG ℤDA DS (t) Objective function of the energy distribution system PEDER i DS (t) Electrical power generation of DSO’s distributed generation units Ki DS(t) Binary decision variable commitment of DSO’s electrical DERs HHDER DS (t) Heating power generation of DSO’s DERs K’DS(t) Binary decision variable commitment of DSO’s heating DERs RCDER DS (t) Cooling power generation of DSO’s DERs K’’DS(t) Binary decision variable commitment of DSO’s cooling DERs CIBEDR DA MG (t) DSO’s costs of IBEDRP paid to MGs CIBHDR DA MG (t) DSO’s costs of IBHDRP paid to MGs CIBCDR DA MG (t) DSO’s costs of IBCDRP paid to MGs CPBEDR DA MG (t) DSO’s costs of PBEDRP paid to MGs CIBHDR DA MG (t) DSO’s costs of PBHDRP paid to MGs CIBCDR DA MG (t) DSO’s costs of PBCDRP paid to MGs CIMPORT DA DS (t) DSO’s costs of electrical energy purchased from the electricity market CPHEVPL DA DS (t) DSO’s costs of electrical energy purchased from PHEV parking lots BPHEVPL DA DS (t) Benefits of electrical energy sold to PHEV parking lots CE IMPORT DA DS (t) DSO’s costs of imported electricity to the electricity wholesale market BE EXPORT DA DS (t) DSO’s benefits of exported electricity from the electricity wholesale market CH IMPORT DA DS (t) DSO’s costs of imported heating energy to the heating energy market BH EXPORT DA DS (t) DSO’s benefits of exported electricity from the heating energy market CC IMPORT DA DS (t) DSO’s costs of imported electricity to the cooling energy market BC EXPORT DA DS (t) DSO’s benefits of exported electricity from the cooling energy market CIDG DS (t) Operating costs of DSO’s intermittent electricity facilities VCPBEDR DA MG (t) Variable costs of PBEDRP paid to MGs VCPBHDR DA MG (t) Variable costs of PBHDRP paid to MGs VCPBCDR DA MG (t) Variable costs of PBCDRP paid to MGs PE Market DS (t) Exported/imported power to/from the electricity market PCCHP DS (t) Generated electrical power of CCHP PESS DS (t) Transacted electrical power with electrical energy storages PCCH DS (t) Electricity consumption of the compressed chiller HCCHP DS (t) Generated heating power of DSO’s CCHP units HH Market DS (t) Exported/imported heating power to/from the heating energy market HBoiler DS (t) Generated heating power of the boiler HTES DS (t) Heating power transactions of thermal energy storages HACH DS (t) Heating energy consumption of absorption chiller HPBDR DS (t) PBHDRP heating loads RCCHP DS (t) Generated cooling power of DSO’s CCHP units RC Market DS (t) Exported/imported power to/from the cooling market RACH DS (t) Generated cooling power of absorption chillers RCCH DS (t) Generated cooling power of compressed chillers RPBDR DS (t) PBCDRP cooling loads PMG(t) Generated electrical power of MG’s DERs KMG(t) Binary decision variable of unit commitment of MG’s electrical DERs HHDER MG (t) Generated heating power of MG’s DERs K’MG(t) Binary decision variable of unit commitment of MG’s heating DERs RHCDER MG (t) Generated heating power of MG’s DERs K’’MG(t) Binary decision variable of unit commitment of MG’s cooling DERs BIBEDR DA MG (t) MGs benefits of IBEDRP BIBHDR DA MG (t) MGs benefits of IBHDRP BIBCDR DA MG (t) MGs benefits of IBCDRP BPBEDR DA MG (t) MGs benefits of PBEDRP BPBHDR DA MG (t) MGs benefits of PBHDRP BPBCDR DA MG (t) MGs benefits of PBCDRP CIMPORT DA MG (t) MG’s costs of electrical energy purchased from the distribution system BEXPORT DA MG (t) MG’s benefit of electrical, heating, and cooling energy sold to the distribution system CIDG MG (t) operating costs of MG’s intermittent electricity facilities PenaltyIBEDR DA MG (t) Penalty costs of MG’s IBEDRP PenaltyIBHDR DA MG (t) Penalty costs of MG’s IBHDRP PenaltyIBCDR DA MG (t) Penalty costs of MG’s IBCDRP PDG MG(t) Generated electrical power of distributed generation units of MG PDSO MG (t) MG’s exported (imported) power to (from) the electrical distribution system PCCHP MG (t) Generated electrical power of MG’s CCHP PPHEV MG (t) Transacted electrical power with MG’s PHEV parking lots PPBDR MG (t) Price-based demand response programs MG’s electrical load PCCH MG (t) Electricity consumption of MG’s compressed chiller HCCHP MG (t) generated heating power of MG’s CCHP units HBoiler MG (t) generated heating power of MG’s boiler HPBDR MG (t) MG’s PBHDRPs heating loads RCCHP s MG (t) Generated cooling power of MG’s CCHP units RACH s MG(t) Generated cooling power of absorption chillers RCCH s MG(t) Generated cooling power of compressed chillers RPBDR b,s MG(t) MG’s PBCDRPs cooling loads P0,H0,R0 Unserved loads of electrical, heating, and cooling systems XEL,XHCV ,XCCV Switching status of electrical switches, heating and cooling system control valves with the energy system can change the volume and availability of DERs’ of the energy system in normal and shock conditions. Further, the switching of electrical switches and heating and cooling valves can change the energy flow of the energy system’s resources, which should be considered a resilient operational scheduling tool. The optimal operational scheduling of energy systems is highly considered in recent papers and the literature as described in the following. Ref. [1] proposed a bi-level active and reactive power scheduling process. The model considered the locational marginal prices. The upper-level problem determined the optimal dispatch of Combined Cooling, Heating, and Power (CCHP) units; meanwhile, the lower-level problem performed the market-clearing process. The model utilized the Big-M and duality theories to convert the bi-level Mathematical Pro graming with Equilibrium Constraints (MPEC) model to a single-level mixed-integer second-order cone-programming model. The case studies presented that the energy efficiency and renewable energy contributions were increased. Ref. [2] introduced a joint chance constraint method to minimize the operating costs of power-to-gas fa cilities and capture carbon dioxide. The uncertainties of intermittent electricity generation units, electricity market price, and loads were modeled. The proposed model reduced carbon emissions by about 50%. However, the power-to-gas process increased the operating costs of the system by about 10%. Refs. [1,2] did not consider the resiliency of the energy system and the switching process of control valves. Ref. [3] explored an optimization process for power, cooling, heat ing, and hydrogen energy systems. The model assessed a Mixed-Integer Nonlinear Programming (MINLP) model to minimize the capital costs, operating costs, maintenance costs, and environmental emissions. The Pareto optimal solution was found using the proposed model, in which the degree of objective function satisfaction was about 0.75. Further, the method revealed that the utilization of energy storage reduced envi ronmental emissions by about 45%. Ref. [4] assessed an adiabatic compressed air energy storage technology for CCHP dispatch. The model considered the temperature dynamic behavior of the energy storage facilities and Mixed-Integer Linear Programming (MILP) was used to solve the problem. The proposed model reduced the operating costs of the system by about 13.7% if compared with the conventional method optimization process. Further, the system costs were reduced by about 48.35%, when the compressed air energy storage was utilized. Refs. [3,4] did not model microgrid energy transactions with the energy system and heating and cooling markets. Ref. [5] proposed a multi-objective problem model to minimize operating costs and emission pollution of a CCHP-based microgrid. The model utilized a scenario-based approach and robust optimization to consider the uncertainties. The output results indicated that the model reduced the emission cost by about 55%. However, the optimization process increased the operating costs by about 11% to achieve zero-level risk. Ref. [6] showed that the incentive-based participation of multi- energy microgrids in energy markets reduced the operating costs by about 4.3$. The model considered energy storage and intermittent power generation facilities and Chance-Constrained Programming (CCP) was used to solve the problem. The optimization method increased the operating costs by about 4.3% to guarantee the reliable performance of the system with a probability of 98%. Refs. [5,6] did not model the resilient operation of the system and the switching process of control valves. Ref. [7] explored different Incentive-Based Demand Response (IBDR) alternatives for microgrids. The resilient response of the system was simulated and different economy and reliability indices were studied for normal and resilient conditions. The model considered the commitment of microgrids, which minimized the microgrids’ energy procurement costs; meanwhile, maximized their profits. The proposed model increased the expected profit of the system by about 4% and 2.7% for normal and resilient conditions, respectively. Ref. [8] utilized an incentive-based trading mechanism for a CCHP-based MG cluster utilizing a decentralized solution model. An asymmetric Nash bargai ning theory was used to determine the benefit of microgrids. The model utilized a contribution rate that was defined as the economic value of exchanged energy. Then, based on the contribution rates of microgrids, their benefits were determined. Refs. [7,8] did not assess the heating and cooling energy transactions of MGs with the energy system. Ref. [9] introduced a coordinated operation of a CCHP-based in dustrial energy park to minimize operation costs and CO2 emissions. A hybrid stochastic/information gap decision theory was used to solve the problem. The uncertainties of intermittent energy generations, loads, and prices were considered. The results revealed that total operation cost and emission pollution were reduced by about 3.1%, and 2.2%, respectively. Ref. [10] explored an integrated model for scheduling intermittent power generation facilities, demand response alternatives, energy storage, and electric vehicles. The model considered gas, elec trical, and heat markets and a robust optimization method was used to solve the problem. The process reduced operational costs by about 14.2% considering the electric vehicle contributions. Further, the model reduced the operating costs by about 13% using demand response pro cedures. Refs. [9,10] did not consider the resiliency operation of the energy system. Ref. [11] proposed an optimization procedure for a CCHP-based system that utilized compressed air energy storage and an internal combustion engine. The model used a bi-level optimization process, in which the upper level optimized the active equipment ca pacity of the system and the lower level minimized the primary energy saving ratio and energy consumption cost. The model reduced the pri mary energy saving ratio by 12.07%. Ref. [12] evaluated a real-time demand response process with optimal scheduling of CCHP-based microgrids. The objective functions minimized operating costs and carbon emissions. The model considered real-time demand response processes, energy hubs, and intermittent energy generation facilities. The method found the optimal Pareto solutions using the max-min fuzzy and weighted sum approach. The output revealed that the carbon emission and operation costs were reduced by 2.26% and 3.97%, respectively. Refs. [11,12] did not model heating and cooling energy markets and the resiliency of energy systems. Ref. [13] proposed an optimization process that preserved informa tion privacy and considered the independent scheduling of microgrids. The uncertainties of multi-carrier energy loads and intermittent energy generations were modeled and a CCP method was used to solve the problem. Further, the cost allocation was modeled to increase the sta bility of the microgrids coalition. Ref. [14] presented a meta-heuristic optimization method for optimal design and energy management of DERs in microgrids. The model considered energy storage, boiler, and intermittent power generation facilities. A Particle Swarm Optimization (PSO) process was carried out. The process reduced fuel consumption and emission pollution by about 6.2% and 46.8%, respectively. Refs. [13,14] did not model the heating and cooling markets, resilience operating conditions of the energy system, and switching of control valves. Ref. [15] introduced a two-stage CCHP microgrid energy manage ment process that comprised by economic dispatch and real-time adjustment processes. The first-stage and second-stage problems were solved using MILP and Mixed Integer Quadratic Programming (MIQP) solvers, respectively. The proposed model reduced the system’s costs by about 16.51% and 11.05% concerning the following thermal load mode and the following electric load mode, respectively. Ref. [16] proposed an optimization process for active power, reactive power, and topology of the system to minimize operating costs, energy loss, and emission pollutants. The optimization process utilized a hybrid big crunch ma chine to optimize the problem. The proposed method reduced the operating costs of the system by about 49.03% concerning the custom optimal power flow method. Refs. [15,16] did not explore the optimal resilient operational scheduling of energy systems considering the heating and cooling markets. Ref. [17] assessed a dynamic robust optimal scheduling process that 5 Table 1 Comparison of the proposed method with other papers. References 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 Proposed Approach Heating market and/or cooling market ✘ ✔ ✘ ✘ ✘ ✘ ✘ ✘ ✔ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ Switching of control valves ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ Switching of electrical switches ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ Resilient operation of system ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ Method MILP ✔ ✘ ✘ ✔ ✔ ✘ ✔ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ ✘ ✘ ✘ ✘ ✘ ✘ MINLP ✘ ✘ ✔ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ ✘ ✘ ✘ ✘ ✘ ✘ ✘ Heuristic ✘ ✔ ✘ ✘ ✔ ✘ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✘ ✘ ✔ ✔ ✔ ✔ ✔ ✔ Model Determ. ✔ ✘ ✔ ✔ ✘ ✘ ✔ ✔ ✘ ✘ ✔ ✔ ✔ ✔ ✔ ✘ ✘ ✔ ✘ ✘ ✔ ✔ ✘ Stoch. ✘ ✔ ✘ ✘ ✔ ✔ ✘ ✘ ✔ ✔ ✘ ✘ ✘ ✘ ✘ ✔ ✔ ✘ ✔ ✔ ✘ ✘ ✔ Objective Function Revenue ✘ ✘ ✘ ✘ ✘ ✘ ✔ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ Gen. Cost ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ Storage Cost ✔ ✔ ✔ ✔ ✔ ✔ ✘ ✔ ✔ ✔ ✔ ✘ ✘ ✘ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ Secu. Costs ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ PHEV cost ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ DRP costs ✘ ✔ ✘ ✘ ✘ ✘ ✔ ✘ ✔ ✘ ✘ ✔ ✔ ✔ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ WT ✘ ✔ ✔ ✘ ✔ ✔ ✘ ✔ ✘ ✔ ✘ ✔ ✔ ✘ ✘ ✔ ✔ ✔ ✔ ✘ ✘ ✘ ✔ PV ✔ ✔ ✔ ✘ ✔ ✔ ✘ ✔ ✘ ✘ ✘ ✔ ✔ ✔ ✘ ✔ ✔ ✔ ✔ ✘ ✘ ✘ ✔ DA-Market ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ scheduling RT- Market ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ ✘ ✘ ✘ ✔ ✘ ✘ ✘ ✔ scheduling Uncertainty Model PHEV ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ DRP ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ DA Market price ✘ ✔ ✘ ✘ ✔ ✔ ✘ ✘ ✔ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ ✘ ✔ ✘ ✘ ✘ ✔ RT Market price ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ Contingency ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✔ Loads ✘ ✔ ✘ ✘ ✔ ✔ ✘ ✘ ✔ ✘ ✘ ✘ ✘ ✘ ✘ ✔ ✔ ✘ ✔ ✘ ✘ ✘ ✔ Inter. Electricity ✘ ✔ ✘ ✘ ✔ ✔ ✘ ✘ ✔ ✘ ✘ ✘ ✘ ✘ ✘ ✔ ✔ ✘ ✔ ✘ ✘ ✘ ✔ determines the optimal operating strategy of the system for day-ahead and intra-day horizons. The objective functions include the operating cost, the cost of energy loss, the energy purchasing cost, and the emis sion pollutant cost. The method reduced the system costs by about 4.77% concerning the worst-case dispatch case. Ref. [18] studied the effectiveness of a collaborative energy management model to control heating, cooling, and electric loads. A Q-learning method was utilized to solve the problem. The model results were compared with the Genetic Algorithm (GA) and PSO results. The proposed method presented the lowest regulation cost, balanced power, and load fluctuations concern ing the PSO and GA results. The model reduced the operating costs of the system by about 14.77% and 15.01% concerning the PSO and GA out puts, respectively. Refs. [17,18] did not assess the optimal scheduling of energy systems considering the heating and cooling markets. Ref. [19] modeled the optimal scheduling of the system for day- ahead, intra-day, and real-time horizons using a two-stage robust opti mization model. A CCG method was utilized to change the problem into a min–max-min problem. The results revealed that the intra-day prob lem reduced the costs of the worst-case scenario and the system reached the secure operating point in this case. Further, based on the proposed model, when the uncertain budget parameter of the system was increased to its maximum value, the increase in the system cost was the minimum. Ref. [20] proposed a robust energy management process for microgrids utilizing the Conditional Value-at-Risk (CVaR) model. An iterative sequential optimization process was used to solve the problem. The model utilized the linear decision rule to assess the uncertainty of energy flows. The method reduced the operating costs of the system by about 9.76% concerning the robust joint chance-constrained method. Refs. [19,20] did not consider the resilient operational scheduling of multi-carrier energy systems and switching of electrical switches and control valves. Ref. [21] introduced a model to peak-load reduction and economic optimization. The linearization process was used to change the non- linear problem into a MILP model. Then, an ε-constraint method was utilized to solve the problem. Finally, the Pareto optimal set was generated. The computation time was reduced by about 71.2%; mean while, the computational accuracy was improved by about 12.73%. The outputs of the MILP optimization process were compared with the out puts of GA, PSO, mutation PSO, and immune GA. The MILP had ad vantages in both computation time and accuracy in solving the problem. Ref. [22] evaluated an optimization model for a solar-driven CCHP system that determined the parameters and configuration of facilities. The sustainability indicator, annual total costs, and energy-saving ratio were utilized to assess the system’s alternatives. Then, the multi- objective problem was modeled as a Pareto front optimization process. The case study was performed for a hotel and office building. The out puts presented that the annual costs of hotel and office buildings were reduced by about 7.69% and 7.33%, respectively. Refs. [21,22] did not Fig. 1. Energy distribution system with its distributed energy resources. model the resilient operational scheduling of the energy system, switching of control valves, and real-time scheduling of the system. 1.2. Novel contribution Based on the Literature Review subsection, Refs. [1–22] did not assess the energy transactions of multi-carrier MGs with the energy system considering the heating and cooling markets. Further, Refs. [1–22] did not optimize the electrical, heating, and cooling systems’ topology considering the contingent condition of the system and the optimal commitment of MGs in this condition. Thus, as shown in Table 1, a framework that considers the simulta neous optimal operation of multi-carrier MGs and their energy trans actions with the energy system considering the heating and cooling markets, optimal switching of electrical switches, and heating and cooling control valves is less frequent in the recent literature. An integrated model for incentive-based and price-based demand response processes that explores the impacts on the resilient operation of the energy system considering the switching of multi-carrier systems has been thus developed in this paper. The main novel contributions of this paper can be presented as: • The proposed method considers the electrical, heating, and cooling energy transactions of the energy system with multi-carrier energy markets in the day-ahead and real-time horizons. Further, the multi- carrier energy transactions of MGs with the energy systems are modeled. • The optimal switching of electrical switches and heating and cooling control valves in shock conditions is modeled. • A resiliency index is proposed to assess the resiliency level of the energy system. • The proposed solution process is decomposed into two-stage two- level optimization processes that optimize the energy system and MGs operational scheduling for day-ahead and real-time horizons. 2. Problem modeling and formulation Fig. 1 shows the multi-carrier energy transactions of the energy distribution system with energy markets, microgrids, and system loads. The energy system has multiple electrical, heating, and cooling energy distributed energy resources, Plug-in Hybrid Electric Vehicles (PHEVs) parking lots, and energy storage facilities. It is assumed that the optimal scheduling of distributed energy re sources and load commitment is carried out by the Distribution System Operator (DSO) through smart grid infrastructure [23]. Further, the DSO utilizes two categories of demand response programs [24,25]: 1) incentive-based demand response programs, and 2) price-based demand response programs. In the incentive-based demand response programs, the MGs gain profit based on their energy consumption reduction based on the multiplier of real-time energy market price [25]. Each MG submits a maximum value for its energy consumption reduction. The price-based demand response programs consist of 1) fixed tariff, 2) Time-Of-Use (TOU) tariff, and 3) real-time tariff. The energy consumption reduction of each consumer in TOU and real-time pricing processes is rewarded by the value of TOU and real- time market prices multiplied by the reduced energy consumption, respectively. However, the TOU and real-time tariff model do not utilize the submission process as defined for the incentive-based demand response process. 2.1. The wind turbine and photovoltaic arrays models The Riley model is utilized to model wind speed. Eq. (1) presents the electrical power of a wind turbine: Pw(vwind) = ⎧ ⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎩ 0 vwind < vci PR (vwind − vci) (vr − vci) vci ≤ vwind < vr PR vr ≤ vwind ≤ vco 0 vwind ≥ vco (1) where, vci,vr, vco, and vwind are cut-in wind velocity, nominal wind ve locity, cut-out wind velocity, and wind velocity, respectively [23]. The maximum power output of photovoltaic arrays is written as (2) [23]: PPV = APV .η.I.(1 − 0.005×(t0 − 25) ) (2) where, APV , I, t0, ηare the area of photovoltaic array, solar irradiation, outside air temperature, and photovoltaic array conversion efficiency, respectively. The CCHPs, boilers, and energy storage facilities models are avail able in [23] and are not presented for the sack of space. 2.2. The electrical demand response process model Based on the above description of the incentive-based process, the Incentive-Based Electrical Demand Response Program (IBEDRP) model considers that the aggregated electrical energy consumption reduction of each MG should be less than its submitted bid [25]. Based on this model, each microgrid can submit a load reduction bid for the next day- ahead operational scheduling process. Then, the DSO optimizes the energy system’s day-ahead operational scheduling problem and de termines the optimal values of load reduction bids of microgrids. Thus, Eq. (3), Eq. (4), and Eq. (5) can be presented for IBEDRP of residential, commercial, and industrial microgrids, respectively. AELRC(r, t) = SELRC(r, t).αr,tSELRC(r, t) ≤ SELRCmax t (3) where AELRC(r, t),SELRC(r, t),αr,t ,SELRCmax t are the actual electrical load reduction of residential MG, the submitted value of electrical load reduction of residential MG, a pre-defined multiplier, and the maximum value of electrical load reduction of residential MG, respectively. AELCC(c, t) = SELCC(c, t).αc,tSELCC(c, t) ≤ SELCCmax t (4) where AELCC(c, t), SERCC(c, t), αc,t , SERCCmax t are the actual electrical load reduction of commercial MG, the submitted value of electrical load reduction of commercial MG, a pre-defined multiplier, and the maximum value of electrical load reduction of commercial MG, respectively. AELIC(i, t) = AELIC(i, t).αi,tSELIC(i, t) ≤ SELICmax t (5) where AELIC(i, t), SELIC(i, t), αi,t , SELICmax t are the actual electrical load reduction of industrial MG, the submitted value of electrical load reduction of industrial MG, a pre-defined multiplier, and the maximum value of electrical load reduction of industrial MG, respectively. The energy management system of MGs can submit different values of electrical load reduction steps for the operational scheduling process. 2.3. The heating demand response process model The Incentive-Based Heating Demand Response Program (IBHDRP) model formulation is the same as IBEDRP. Thus, Eq. (6), Eq. (7), and Eq. (8) can be presented for IBHDRP of residential, commercial, and in dustrial microgrids, respectively. AHLRC(r, t) = SHLRC(r, t).α’r,tSHLRC(r, t) ≤ SHLRCmax t (6) where AHLRC(r, t),SHLRC(r, t),α’r,t ,SHLRCmax t are the actual heating load reduction of residential MG, the submitted value of heating load reduction of residential MG, a pre-defined multiplier, and the maximum value of heating load reduction of residential MG, respectively. AHLCC(c, t) = SHLCC(c, t).α’c,tSHLCC(c, t) ≤ SHLCCmax t (7) where AHLCC(c, t), SHRCC(c, t), α’c,t, SHRCCmax t are the actual heating load reduction of commercial MG, the submitted value of heating load reduction of commercial MG, a pre-defined multiplier, and the maximum value of heating load reduction of commercial MG, respec tively. AHLIC(i, t) = AHLIC(i, t).α’i,tSHLIC(i, t) ≤ SHLICmax t (8) where AHLIC(i, t), SHLIC(i, t), α’i,t , SHLICmax t are the actual heating load reduction of industrial MG, the submitted value of heating load reduc tion of industrial MG, a pre-defined multiplier, and the maximum value of heating load reduction of industrial MG, respectively. 2.4. The cooling demand response process model The Incentive-Based Cooling Demand Response Program (IBCDRP) model formulation is the same as IBEDRP. Thus, Eq. (9), Eq. (10), and Eq. (11) can be presented for IBCDRP of residential, commercial, and industrial microgrids, respectively. ACLRC(r, t) = SCLRC(r, t).α’’r,tSCLRC(r, t) ≤ SCLRCmax t (9) where ACLRC(r, t), SCLRC(r, t),α’’r,t , SCLRCmax t are the actual cooling load reduction of residential consumers, the submitted value of cooling load reduction of residential consumers, a pre-defined multiplier, and the maximum value of cooling load reduction of residential consumers, respectively. ACLCC(c, t) = SCLCC(c, t).α’’c,tSCLCC(c, t) ≤ SCLCCmax t (10) where ACLCC(c, t), SCRCC(c, t), α’’c,t , SCRCCmax t are the actual cooling load reduction of commercial consumers, the submitted value of cooling load reduction of commercial consumers, a pre-defined multiplier, and the maximum value of cooling load reduction of commercial consumers, respectively. ACLIC(i, t) = ACLIC(i, t).α’’i,tSCLIC(i, t) ≤ SCLICmax t (11) where ACLIC(i, t), SCLIC(i, t), α’’i,t , SCLICmax t are the actual cooling load reduction of industrial consumers, the submitted value of cooling load reduction of industrial consumers, a pre-defined multiplier, and the maximum value of cooling load reduction of industrial consumers, respectively. 2.5. The objective function of DSO for day-ahead operational scheduling problem (First level of day-ahead optimization process) As shown in Fig. 1, the DSO transacts electrical, heating, and cooling energy carriers with the wholesale electricity market, heating energy market, and cooling energy market, respectively in the day ahead and real-time horizons. The day-ahead optimization process determines the optimal scheduling of multi-carrier energy systems and MGs’ DERs considering the forecasted values of multi-carrier energy consumptions, energy market prices, and intermittent energy generations. However, based on the stochastic behavior of intermittent energy generation and volatility of energy carriers’ prices of electricity, heating, and cooling markets there are mismatches between the day-ahead optimal sched uling of DERs and the real-time operating conditions. Thus, a revised scheduling of energy systems and MGs’ DERs should be carried out to Fig. 2. The proposed framework for the optimization process. assess the impacts of energy generation and consumption mismatches. Hence, the optimization process should have two-stage optimization procedures, and these processes determine the optimal scheduling of the energy system and MGs in the day-ahead and real-time phases, respec tively. The optimization framework is presented in Fig. 2. At the first stage optimization process, the DSO and MGs are optimizing their day- ahead objective functions in the first and second levels, respectively. Then, the real-time data are updated and in the second stage of opti mization, the DSO and MGs optimize their real-time objective functions in the first and second levels, respectively. The optimization processes of the optimization framework are presented in the following subsections. Based on Fig. 1 and Fig. 2, the DSO simultaneously transacts multi- carrier energy with microgrids and energy markets. Thus, the DSO can gain profit by selling energy carriers to MGs and energy markets. Thus, the distribution system operator should maximize the system reliability; meanwhile, minimizing the system operating costs. Thus, the objective function of the distribution system operator can be formulated as (12): Min ℤDA DS (t) = ∑T t=1 GDA DS (t)+ ∑T t=1 ∑DSDAS s=1 prs⋅ADA s DS(t)+ ∑T t=1 ∑DSRTS k=1 prk⋅ℤRT DS (12) where, GDA DS (t) consists of the following terms. Eq. (13) consists of the operating costs of DSO electrical distributed energy resources (PEDER i DS (t)⋅γi DS(t)⋅Ki DS(t)), DSO’s start-up costs (SUi DS(t)⋅[Ki DS(t) − Ki DS(t − 1) ]), operating costs of DSO heating en ergy distributed energy resources (HHDER i’DS (t)⋅γ’i’DS(t)⋅K’i’DS(t)), operating costs of DSO cooling energy resources (RCDER i’’DS (t)⋅γ’’i’’DS(t)⋅K’’i’’DS(t)), the costs of IBEDRP paid to MGs (CIBEDR DA j MG (t)), the costs of IBHDRP paid to MGs (CIBHDR DA j’MG (t)), the costs of IBCDRP paid to MGs (CIBCDR DA j˝MG (t)), the costs of Price-Based Electrical Demand Response Program (PBEDRP) paid to MGs (CPBEDR DA g MG (t)), the costs of Price-Based Heating Demand Response Program (PBHDRP) paid to MGs (CPBHDR DA g’MG (t)), the costs of Price-Based Cooling Demand Response Program (PBCDRP) paid to MGs (CPBCDR DA g’’MG (t)), costs of electrical energy purchased from wholesale electricity market (CIMPORT DA DS (t)), costs of electrical energy purchased from PHEV parking lots (CPHEVPL DA DS (t)), and benefits of electrical energy sold to PHEV parking lots (BPHEVPL DA DS (t)). CE IMPORT DA DS (t), BE EXPORT DA DS (t)are the DSO’s imported electricity costs and exported electricity benefits, respectively. CH IMPORT DA DS (t), BH EXPORT DA DS (t)are the DSO’s imported heating energy costs and expor ted heating energy benefits, respectively. CC IMPORT DA DS (t), BC EXPORT DA DS (t)are the DSO’s imported cooling energy costs and expor ted cooling energy benefits, respectively. Where, Pi DS(t), γi DS(t),Ki DS(t) are DSO’s generated power of distributed energy resources, cost of energy generation, and binary de cision variable of DSO’s DER commitment, respectively. The ADA s DS(t) is decomposed into the following terms. Eq. (14) consists of the operating costs of DSO’s intermittent elec tricity facilities (CIDG (i,s) DS(t)), the variable costs of PBEDRP paid to MGs (VCPBEDR DA MG (t)), the variable costs of PBHDRP paid to MGs (VCPBHDR DA MG (t)), the s variable costs of PBCDRP paid to MGs (VCPBCDR DA MG (t)), the electrical energy not supplied costs (EENSCS DS(t)), the heating energy not supplied costs (HENSCS DS(t)), the cooling energy not supplied costs (CENSCS DS(t)). The proposed second objective function can be presented as (15): BDA DS (t) = TENPDG DS(t)+ TENPMarket DS(t) (15) Eq. (15) consists of the volume of environmental pollution of DSO’s distributed generation units and the fossil-fuelled electricity generation of wholesale market units, respectively. Where, TENPDG DS(t) and TENPMarket(t) are formulated as Eq. (16) and Eq. (17), respectively. GDA DS (t) = ∑NDG i=1 [ PEDER i DS (t)⋅γi DS(t)⋅Ki DS(t) + SUi DS(t)⋅[Ki DS(t) − Ki DS(t − 1) ] ] + ∑NHDER i′=1 [ HHDER i DS (t)⋅γ′i′DS(t)⋅K′i′DS(t) ] + ∑NCDER i′′=1 [ RCDER i DS (t)⋅γ’′i′′DS(t)⋅K’′i′′DS(t) ] + ∑J j=1 CIBEDR DA j MG (t) + ∑J′ j′=1 CIBHDR DA j′MG (t) + ∑J˝ j′′=1 CIBCDR DA j˝MG (t) + ∑G g=1 CPBEDR DA g MG (t)+ ∑G′ g′=1 CPBHDR DA g′MG (t) + ∑G′′ g′′=1 CPBCDR DA g′′MG (t) + CPHEVPL DA DS (t) − BPHEVPL DA DS (t)+ CE IMPORT DA DS (t) − BE EXPORT DA DS (t) + CH IMPORT DA DS (t) − BH EXPORT DA DS (t)+ CC IMPORT DA DS (t) − BC EXPORT DA DS (t) (13) ADA s DS(t) = ∑NIDG i=1 CIDG (i,s) DS(t) + ∑G g=1 VCPBEDR DA (g,s) MG (t) + ∑G′ g′=1 VCPBHDR DA (g′,s) MG (t) + ∑G′′ g′′=1 VCPBCDR DA (g′′,s) MG (t)+ EENSCS DS(t) + HENSCS DS(t) + CENSCS DS(t) (14) TENPDG DS(t) = ∑NDG i=1 [ ENPDG CO2 DS(t) +ENPDG SO2 DS(t) +ENPDG NOx DS(t) ) ⋅PDG DS (t) ] (16) where, ENPDG CO2(t),ENPDG SO2(t),ENPDG NOx(t)are the volumes of CO2, SO2, and NOx pollution of distributed generation units, respectively. TENPMarket(t) = ENPMarket CO2 Market(t)+ENPMarket SO2 Market(t)+ENPMarket NOx Market(t) (17) where, ENPMarket CO2 Market(t), ENPMarket SO2 Market(t), ENPMarket NOx Market(t)are the volumes of CO2, SO2, and NOx pollution of wholesale market generation units, respectively. The ℤRT DS is the objective function of DSO’s real-time optimization problem. The proposed objective function has the following constraints. A. Electric power balance constraint The electrical power balance constraint should be satisfied for any simulation interval as (18): ∑NDG i=1 PEDER i,s DS(t) ∓ PE Market s DS (t) + PCCHP s DS (t) ∓ PESS s DS(t) ∓ PPHEVPL s DS (t) = ∑NEL m=1 PL m,s DS(t) − ∑NRMG a=1 AELRCa,s MG(t) − ∑NCMG a′=1 AELCCa′,s MG(t)− ∑NIMG a′′=1 AELICa′′,s MG(t) − ∑NEIBDR b=1 PPBDR b,s DS(t) + PCCH s DS (t) (18) Eq. (18) presents that the aggregated generated electrical power of DSO’s distributed generation units (PEDER i,s DS(t)), the exported (imported) power to (from) electricity market (PE Market s DS (t)), generated power of CCHP (PCCHP s DS (t)), transacted electrical power with electrical energy storages (PESS s DS(t)), and transacted electrical power with PHEV parking lots (PPHEVPL s DS ) should be equal to the aggregated active power of loads PL m,s DS(t), residential MG load reduction AELRCa,s MG(t), commercial MG load reduction (AELCCa’,s MG(t)), industrial MG load reduction (AELICa’’,s MG(t)), electrical loads participating in price-based demand response programs (PPBDR b,s DS(t)), and the electricity consumption of com pressed chiller (PCCH s DS(t)). B. Heating power balance constraint The heating power balance constraint of DSO should be satisfied for any simulation interval as (19): HCCHP s DS (t) ∓ HH Market s DS (t) + HBoiler s DS (t) ∓ HTES s DS(t) = ∑NHL m=1 HL m,s DS(t) − ∑NRMG a=1 AHLRCa,s MG(t) − ∑NCMG a′=1 AHLCCa′,s MG(t)− ∑NIMG a′′=1 AHLICa′′,s MG(t) − ∑NHIBDR b=1 HPBDR b,s DS(t) + HACH s DS (t) (19) Eq. (19) presents that the aggregated generated heating power of DSO’s CCHP units (HCCHP s DS (t)), the exported (imported) power to (from) heating energy market (HH Market s DS (t)), generated heating power of the boiler (HBoiler s DS (t)), and the heating power transactions of the thermal energy storages (HTES s DS(t)) should be equal to the aggregated heating power of loads (HL m,s DS(t)), residential MG heating load reduction (AHLRCa,s MG(t)), commercial MG heating load reduction (AHLCCa’,s MG(t)), industrial MG heating load reduction (AHLICa’’,s MG(t)), PBHDRP heating loads (HPBDR b,s DS(t)), and the heating energy consumption of absorption chiller (HACH s DS(t)). C. Cooling power balance constraint The cooling power balance constraint of DSO should be satisfied for any simulation interval as Eq. (20): RCCHP s DS (t) ∓ RC Market s DS (t) + RACH s DS(t) + RCCH s DS (t) = ∑NCL m=1 RL m,s DS(t) − ∑NRMG a=1 ACLRCa,s MG(t) − ∑NCMG a′=1 ACLCCa′,s MG(t) − ∑NIMG a′′=1 ACLICa′′,s MG(t) − ∑NCIBDR b=1 RPBDR b,s DS(t) (20) Eq. (20) presents that the aggregated generated cooling power of DSO’s CCHP units (RCCHP s DS (t)), the exported (imported) power to (from) cooling energy market (RC Market s DS (t)), and generated cooling power of absorption chillers (RACH s DS(t)) and compressed chillers (RCCH s DS(t)) should be equal to the aggregated cooling power of DSO’s loads (RL m,s DS(t)), resi dential MG cooling load reduction (ACLRCa,s MG(t)), commercial MG cooling load reduction (ACLCCa’,s MG(t)), industrial MG cooling load reduction (ACLICa’’,s MG(t)), and PBCDRP cooling loads (RPBDR b,s DS(t)). The AC load flow constraints for the electrical system, the constraints of facilities’ maximum capacity, charge, and discharge constraints, and mass balance constraints should be considered for each of the simulation intervals and are not presented for the sake of space. 2.6. The objective function of MGs for day-ahead operational scheduling problem (Second level of day-ahead optimization process) As shown in Fig. 1, the MGs transact energy carriers with the DSO. The MG’s operator should minimize the operating costs of MG; mean while, maximizing the energy transactions and demand response pro grams’ contributions profits. Thus, the proposed day-ahead objective functions of MGs consist of the following terms: 1) The fixed and start-up costs of operating costs of MGs’ distributed generation units, demand response program benefits, and costs (benefits) of imported/exported energy from/to energy system, 2) The operating costs of MG’s inter mittent electricity generation considering uncertainties, the variable operating costs of MGs’ distributed generation units, the load commit ment benefits, and energy not supplied costs, 3) the costs of environ mental pollutions, 4) and the expected value of the real-time objective function. The objective function can be presented as (21): Min FDA MG(x) = ∑T t=1 M DA MG(t)+ ∑T t=1 ∑MGDAS s=1 prs⋅ADA s MG(t) + ∑T t=1 SDA MG(t)+ ∑T t=1 × ∑MGRTS k=1 prk⋅FRT MG (21) where, MDA MG(t) is presented as (22): Eq. (21) consists of the MG’s operating costs of electrical distributed generation units (PDG i MG(t)⋅γi MG(t)⋅Ki MG(t)), operating costs of MG’s heat ing energy distributed energy resources (HHDER i’MG (t)⋅γ’i’MG(t)⋅K’i’MG(t)), operating costs of MG’s cooling energy resources (RCDER i’’MG(t)⋅γ’’i’’MG(t)⋅K’’i’’MG(t)), the benefits of IBEDRP (BIBEDR DA j MG (t)), the benefits of IBHDRP (BIBHDR DA j’MG (t)), the benefits of IBCDRP (BIBCDR DA j˝MG (t)), the benefits of PBEDRP (BPBEDR DA g MG (t)), the benefits of PBHDRP (BPBHDR DA g’MG (t)), the benefits of PBCDRP (BPBCDR DA g’’MG (t)), costs of electrical energy purchased from the energy system(CIMPORT DA MG (t)), and benefit of electrical, heating, and cooling energy sold to the energy system (BEXPORT DA MG (t)). Where, Pi MG(t), γi MG(t),Ki MG(t) are generated power of MG distrib uted generation unit, cost of electrical energy generation of MG, and binary decision variable of unit commitment of MG facilities, respec tively. The ADA s MG(t) term consists of Eq. (23) terms. A DA s MG(t) = ∑NIDG′ i=1 CIDG (i,s) MG(t) + ∑J j=1 PenaltyIBEDR DA (j,s) MG (t)+ ∑J′ j′=1 PenaltyIBHDR DA (j′,s) MG (t) + ∑J′′ j′′=1 PenaltyIBCDR DA (j′′,s) MG (t)+ EENSCS MG(t) + HENSCS MG(t) + CENSCS MG(t) (23) Eq. (23) consists of the operating costs of MG’s intermittent elec tricity facilities (CIDG (i,s) MG(t)), the penalty costs of MG’s IBEDRP (PenaltyIBEDR DA (j,s) MG (t)), the penalty costs of MG’s IBHDRP (PenaltyIBHDR DA (j’,s) MG (t)), the penalty costs of MG’s IBCDRP (PenaltyIBCDR DA (j’’,s) MG (t)), the electrical energy not supplied costs of MG (EENSCS MG(t)), the heating energy not supplied costs of MG (HENSCS MG(t)), the cooling energy not supplied costs of MG (CENSCS MG(t)). The third term of the objective function of Eq. (21) can be presented as (24): SDA MG(t) = TENPDG MG(t) (24) Eq. (24) consists of the volume of environmental pollution of distributed generation units of MG, which can be formulated as (25): TENPDG MG(t) = ∑NDG’ i=1 [ ENPDG CO2 MG(t) +ENPDG SO2 MG(t)+ENPDG NOx MG(t) ) ⋅PDG MG(t) ] (25) where, ENPDG CO2(t),ENPDG SO2(t),ENPDG NOx(t)are the volumes of CO2, SO2, and NOx pollution of MG’s distributed generation unit, respectively. The FRT MG is the objective function of MG’s real-time optimization problem. The proposed objective function has the following constraints. A. Electric power balance constraint The electrical power balance constraint of MG should be satisfied for any simulation interval, which can be formulated as (26): ∑NDG′ i=1 PDG i,s MG(t) ∓ PDSO s MG(t) + PCCHP s MG (t) ∓ PESS s MG(t) ∓ PPHEV s MG (t) = ∑NEL′ m=1 PL m,s MG(t) − PPBDR b,s MG(t) + PCCH s MG(t) (26) Eq. (26) presents that the aggregated generated electrical power of distributed generation units of MG (PDG i,s MG), the exported (imported) power to (from) electrical energy system(PDSO s MG(t)), generated electrical power of MG’s CCHP (PCCHP s MG (t)), transacted electrical power with elec trical storages (PESS s MG(t)), and transacted electrical power with PHEVs’ parking lots (PPHEV s MG (t)) should be equal to the aggregated active power of MG’s loads (PL m,s MG), price-based demand response programs MG’s electrical loads (PPBDR b,s MG(t)), and the electricity consumption of com pressed chiller (PCCH s MG(t)). B. Heating power balance constraint The heating power balance constraint should be satisfied for any simulation interval, which can be formulated as (27): HCCHP s MG (t) + HBoiler s MG (t) = ∑NHL′ m=1 HL m,s MG(t) − HPBDR b,s MG(t) + HACH s MG(t) (27) Eq. (27) presents that the aggregated generated heating power of MG’s CCHP units (HCCHP s MG (t)), and generated heating power of MG’s boiler (HBoiler s MG (t)) should be equal to the aggregated MG’s heating power of loads (HL m,s MG(t)), MG’s PBHDRPs heating loads (HPBDR b,s MG(t)), and the MG’s heating energy consumption of absorption chiller (HACH s MG(t)). C. Cooling power balance constraint The cooling power balance constraint of MG should be satisfied for any simulation interval, which can be presented as (28): RCCHP s MG (t) + RACH s MG(t) + RCCH s MG(t) = ∑NCL′ m=1 RL m,s MG(t) − RPBDR b,s MG(t) (28) M DA MG(t) = ∑NDG′ i=1 [ PDG i MG(t)⋅γi MG(t)⋅Ki MG(t) ] + ∑NHRMG i′=1 [ HHDER i MG (t)⋅γ′i′MG(t)⋅K′i′MG(t) ] + ∑NCRMG i′′=1 [ RCDER i MG (t)⋅γ’′i′′MG(t)⋅K’′i′′MG(t) ] + ∑J j=1 BIBEDR DA j MG (t) + ∑J′ j′=1 BIBHDR DA j′MG (t) + ∑J˝ j′′=1 BIBCDR DA j˝MG (t) + ∑G g=1 BPBEDR DA g MG (t)+ ∑G′ g′=1 BPBHDR DA g′MG (t) + ∑G′′ g′′=1 BPBCDR DA g′′MG (t) + CIMPORT DA MG (t) − BEXPORT DA MG (t) (22) Fig. 3. The flowchart of the proposed procedure. Fig. 4. The topology of the modified 123-bus test system. Fig. 5. The estimated energy consumption of electrical, heating, and cooling loads of the 123-bus system. Eq. (28) presents that the aggregated generated cooling power of MG’s CCHP units (RCCHP s MG (t)), and generated cooling power of absorption chillers (RACH s MG(t)) and compressed chillers (RCCH s MG(t)) should be equal to the aggregated MG’s cooling power of loads (RL m,s MG), and MG’s PBCDRPs cooling loads (RPBDR b,s MG(t). The AC load flow constraints for the electrical system, the constraints of facilities’ maximum capacity, charge, and discharge constraints, and mass balance constraints should be considered for each of the simulation Fig. 6. The estimated values of energy generation of wind turbines and photovoltaic arrays of 123-bus system. Fig. 7. The estimated price of day-ahead electrical, heating, and cooling markets. Fig. 8. The estimated price of real-time electrical, heating, and cooling markets. intervals and are not presented for the sake of space. 2.7. The objective function of DSO for the real-time operational scheduling problem (First level of real-time optimization process) As shown in Fig. 2, the energy system and MGs data are updated for real-time operational scheduling. The DSO should minimize the oper ating costs; meanwhile, maximizing the benefit of energy transactions with the electrical, heating, and cooling markets. Further, the DSO should minimize the unserved load considering the shock conditions. Finally, the DSO should optimize the status of electrical switches and heating and cooling valves. Thus, the distribution system operator’s objective function for the real-time operating conditions can be formu lated as (29): Min ℤRT DS(x) Ω = ∑NRTSS k=1 W1⋅ ( GRT k DS(t) + ART k DS(t) + BRT k DS(t)+ W2⋅ [ ∑NEL k′=1 (1 − P0 k′) + ∑NHL k′′=1 (1 − H0 k′′) + ∑NCL k′′′=1 (1 − R0 k′′′) ] + W3⋅ [ ∑NNCSW j=1 XEL j + ∑NNCHV j′=1 XHCV j′ + ∑NNCCV j′′=1 XCCV j′′ ] ∀Ω ∈ {Normal and Contingent Conditions} (29) where the GRT DS(t) can be represented as (30): Eq. (29) decomposed into DSO’s costs ∑NRTST k=1 W1⋅ ( GRT k DS(t)+ ART k DS(t) + BRT k DS(t) ) , the unserved loads of electrical, heating, and cooling systems [ ∑NEL k’=1(1 − P0 k’) + ∑NHL k’’=1(1 − H0 k’’) + ∑NCL k’’’=1 (1 − R0 k’’’) ] , and the switching of electrical switches, heating, and cooling system control valves [ ∑NNCSW j=1 XEL j + ∑NNCHV j’=1 XHCV j’ + ∑NNCCV j’’=1 XCCV j’’ ] . It is assumed that the DSO can switch the boundary lines’ switches and heating and cooling valves in normal and shock conditions to perform preventive/corrective actions. Eq. (30) consists of the real-time operating costs of DSO electrical, heating energy, and cooling energy distributed generation units (PDG i DS(t)⋅γi DS(t),HHDER i DS (t)⋅γ’i’DS(t),RCDER i DS (t)⋅γ’’i’’DS(t)), the costs of IBEDRP paid to MGs (CIBEDR RT j MG (t)), the costs of IBHDRP paid to MGs (CIBHDR RT j’MG (t)), the costs of IBCDRP paid to MGs (CIBCDR RT j˝MG (t)), the costs of PBEDRP paid to MGs (CPBEDR RT g MG (t)), the costs of PBHDRP paid to MGs (CPBHDR RT g’MG (t)), the costs of PBCDRP paid to MGs (CPBCDR RT g’’MG (t)), costs of electricity, heating, and cooling energy purchased from electrical, heating and cooling energy market (CIMPORT RT DS (t)), and benefit of elec trical energy (BE EXPORT RT DS (t)), heating energy (BH EXPORT RT DS (t)), and cooling energy (BC EXPORT RT DS (t)) sold to the electricity market, heating energy market, and cooling energy market, respectively. The second term of Eq. (29) can be presented as (31): GRT DS(t) = ∑NDG i=1 [ PDG i DS(t)⋅γi DS(t) ] + ∑NHDER i′=1 [ HHDER i DS (t)⋅γ′i′DS(t) ] + ∑NCDER i′′=1 [ RCDER i DS (t)⋅γ’′i′′DS(t) ] + ∑J j=1 CIBEDR RT j MG (t) + ∑J′ j′=1 CIBHDR RT j′MG (t) + ∑J˝ j′′=1 CIBCDR RT j˝MG (t)+ ∑G g=1 CPBEDR RT g MG (t) + ∑G′ g′=1 CPBHDR RT g′MG (t) + ∑G′′ g′′=1 CPBCDR RT g′′MG (t) + CIMPORT RT DS (t) − BE EXPORT RT DS (t) − BH EXPORT RT DS (t) − BC EXPORT RT DS (t) (30) ART s DS(t) = ∑NIDG i=1 CIDG (i,s) DS(t) + ∑G g=1 VCPBEDR RT (g,s) MG (t) + ∑G′ g′=1 VCPBHDR RT (g′,s) MG (t) + ∑G′′ g′′=1 VCPBCDR RT (g′′,s) MG (t)+ ECICDS(t) + HCICDS(t) + CCICDS(t) (31) Table 2 The characteristics of distributed energy resources of microgrids. CCHP Electrical efficiency 0.45 Thermal efficiency 0.3 CCHP Maximum electrical power output 20 kW Maximum thermal power output 20 kW Boiler Efficiency 0.9 Maximum thermal power output 40 kW Electrical chiller Efficiency 0.85 Maximum thermal power output 30 kW Absorption chiller Efficiency 0.90 Maximum thermal power output 25 kW Eq. (31) consists of the operating costs of DSO’s intermittent elec tricity facilities (CIDG (i,s) DS(t)), the variable costs of PBEDRP paid to MGs (VCPBEDR RT (j,s) MG (t)), the variable costs of PBHDRP paid to MGs (VCPBHDR RT (j’,s) MG (t)), the variable costs of PBCDRP paid to MGs (VCPBCDR RT (j’’,s) MG (t)), the electrical energy not supplied costs (ECICDS(t)), the heating energy not supplied costs (HCICDS(t)), the cooling energy not supplied costs (CCICDS(t)). Fig. 9. (a) The electrical loads of MGs for one of the reduced scenarios., (b) The heating loads of MGs for one of the reduced scenarios., (c) The cooling loads of MGs for one of the reduced scenarios. Fig. 10. The average values of the electrical power output of residential, commercial, and industrial MGs’ wind turbines for different scenarios. Fig. 11. The average values of the electrical power output of residential, commercial, and industrial MGs’ photovoltaic systems for different scenarios. Fig. 12. The values of fixed and time-of-use prices of energy system for electrical, heating, and cooling loads. The third term of Eq. (29) can be presented as (32): BRT DS(t) = TENPDG DS(t)+ TENPMarket DS(t) (32) Eq. (32) consists of the volume of environmental pollution of DSO’s distributed generation units and the fossil-fueled electricity generation of wholesale market units, respectively. The proposed objective function has the same constraints as the day- ahead problem. A Multi-Carrier Energy System Resiliency Index (MCESRI) is utilized to assess the resiliency level of the energy system. The MCESRI is defined as (33): MCESRI = ∑NCEL k′=1 PCritical /( ∑NEL k′=1 P0 − ∑NCEL k′=1 PCritical ) + ∑NCHL k′′=1 HCritical /( ∑NHL k′′=1 H0 − ∑NCHL k′′=1 HCritical ) + ∑NCCL k′′′=1 RCritical /( ∑NCL k′′′=1 R0 − ∑NCCL k′′′=1 RCritical ) (33) where, PCritical, HCritical, RCriticalare the critical values of electrical loads, heating loads, and cooling loads, respectively. Further, P0,H0,R0are the sum of critical and non-critical electrical loads, heating loads, and cooling loads, respectively. 2.8. The objective function of MGs for real-time operational scheduling problem (Second level of real-time optimization process) As shown in Fig. 2, the MGs can transact energy carriers with the DSO. The MG operator should minimize operating costs and maximize the profits of energy transactions with the DSO. Further, the MG oper ator should minimize the penalty costs of any deviation between day- ahead accepted values and their real-time operating variables. Finally, the interruption costs of MG’s loads should be minimized. The proposed real-time objective function is formulated as (34): Min FRT MG(x) = ∑NRTST k=1 [ M RT MG(t)+A RT MG(t) +SRT MG(t) ] (34) where MRT MG(t) can be presented as (35): Fig. 13. The value of incentive-based electrical, heating, and cooling demand response programs for MGs. Table 3 The feasible switching alternatives of electric switches and control valves of district heating and cooling networks. Topology Number Number of ES&CV 1 1 0 1 0 0 1 0 0 2 1 1 1 0 1 1 0 0 3 1 0 1 0 0 1 0 0 4 0 0 1 0 0 1 0 0 5 0 1 1 0 0 1 0 0 6 1 0 1 0 0 1 1 1 7 1 0 1 0 0 1 1 0 8 1 0 1 0 0 1 0 0 9 1 0 1 0 0 1 0 0 10 1 0 1 1 1 0 1 0 11 0 0 1 0 0 1 0 0 12 1 0 1 1 0 1 0 0 13 1 0 1 1 0 1 0 0 14 1 1 1 0 0 1 0 0 15 1 0 1 0 0 1 0 1 16 1 0 1 0 0 1 1 1 Eq. (35) consists of the operating costs of distributed generation units (PDG i MG(t)⋅γi MG(t), HHDER i MG (t)⋅γ’i’MG(t), RCDER i MG (t)⋅γ’’i’’MG(t)), the benefits of IBEDRP (CIBEDR RT j MG (t)), the benefits of IBHDRP (CIBHDR RT j’MG (t)), the bene fits of IBCDRP (CIBCDR RT j˝MG (t)), the benefits of PBEDRP (CPBEDR RT g MG (t)), the benefits of PBHDRP (CPBHDR RT g’MG (t)), the benefits of PBCDRP (CPBCDR RT g’’MG (t)), costs of electricity, heating and cooling energy purchased from DSO (CIMPORT RT MG ), and benefit of electrical, heating, and cooling energy sold to DSO (BEXPORT RT MG ). The ART s MG(t) term is formulated as (36): A RT s MG(t) = ∑NIDG′ i=1 CIDG (i,s) MG(t) + ∑J j=1 PenaltyIBEDR RT (g,s) MG (t)+ ∑J′ j′=1 PenaltyIBHDR RT (g′,s) MG (t) + ∑J′′ j′′=1 PenaltyIBCDR RT (g′′,s) MG (t)+ ECICMG(t) + HCICMG(t) + CCICMG(t) (36) Eq. (36) consists of the operating costs of MG’s intermittent elec tricity facilities (CIDG (i,s) MG(t)), the penalty costs of MG’s IBEDRP (PenaltyIBEDR RT (j,s) MG (t)), the penalty costs of MG’s IBHDRP Table 4 The nine cases of demand response alternatives that the simulation processes were carried out. M RT MG(t) = ∑NDG′ i=1 [ PDG i MG(t)⋅γi MG(t) ] + ∑NHRMG i′=1 [ HHDER i MG (t)⋅γ′i′MG(t) ] + ∑NCRMG i′′=1 [ RCDER i MG (t)⋅γ’′i′′MG(t) ] + ∑J j=1 BIBEDR RT j MG (t) + ∑J′ j′=1 BIBHDR RT j′MG (t)+ ∑J˝ j′′=1 BIBCDR RT j˝MG (t) + ∑G g=1 BPBEDR RT g MG (t) + ∑G′ g′=1 BPBHDR RT g′MG (t) + ∑G′′ g′′=1 BPBCDR RT g′′MG (t)+ CIMPORT RT MG (t) − BEXPORT RT MG (t) (35) (PenaltyIBHDR RT (j’,s) MG (t)), the penalty costs of MG’s IBCDRP (PenaltyIBCDR RT (j’’,s) MG (t)), the MG’s electrical consumer interruption costs (ECICMG(t)), the MG’s heating consumer interruption costs (HCICMG(t)), the MG’s heating consumer interruption costs (CCICMG(t)). The SRT MG(t) term is the same as SDA MG(t). The constraints of the real- time proposed objective function are the same as the day-ahead objec tive function constraints and are not presented for the sake of space. 3. Solution algorithm As shown in Fig. 3, the proposed model is an iterative two-stage two- level optimization process. The following assumptions are considered in the method: • The Auto-Regressive Integrated Moving Average (ARIMA) models were utilized to generate scenarios for the following stochastic pro cesses: hourly electrical, heating, and cooling demand profiles, day- ahead prices of electricity, heating, and cooling markets, hourly charge and discharge of PHEVs’ parking lots, MGs’ energy trans actions and energy consumptions, wind turbine, and photovoltaic electricity generation, and real-time prices of electricity, heating, and cooling markets. The detailed formulations of ARIMA models for generating scenarios for stochastic processes are available in [26]. • However, solving the optimization problem for the generated sce narios was very time-consuming. Thus, a scenario reduction process was utilized. The forward selection method was used to reduce the generated scenarios [27]. The total number of scenarios was 10,007, Table 5 The average costs of residential, industrial, and commercial microgrids for the twenty-one case studies. RTP (without IBDRP) RTP (with IBDRP1) RTP (with IBDRP2) TOU (without IBDRP) TOU (with IBDRP1) TOU (with IBDRP2) Fixed Price dirgorci MlairtsudnI Minimization of operating costs Operating costs (MUs) 57121 55219 54902 60660 59064 58847 68447 Emission pollution (kg) 146920 143184 142992 144406 141209 141209 161919 Minimizing of emission pollution Operating costs (MUs) 146041 146115 148839 146041 146115 148839 146041 Emission pollution (kg) 110657 108696 108696 110657 108696 108696 110657 Minimization of operating costs and emission pollution Operating costs (MUs) 57121 55219 54902 60660 59064 58847 68447 Emission pollution (kg) 146920 143184 142992 144406 143209 144209 161919 dirgorci Mlaicre m mo C Minimization of operating costs Operating costs (MUs) 146041 146115 148839 146041 146115 148839 146041 Emission pollution (kg) 110657 108696 108696 110657 108696 108696 110657 Minimizing of emission pollution Operating costs (MUs) 69855 66945 66772 71966 69764 69821 100454 Emission pollution (kg) 126267 125661 125192 125501 124661 124212 116834 Minimization of operating costs and emission pollution Operating costs (MUs) 69117 66293 65333 68546 72057 73559 76661 Emission pollution (kg) 174834 161821 158721 167511 168039 163802 202399 dirgorci Mlaitn edise R Minimization of operating costs Operating costs (MUs) 169407 160726 181583 163566 182644 174141 181090 Emission pollution (kg) 125042 122826 128261 128362 119565 127174 138321 Minimizing of emission pollution Operating costs (MUs) 82429 79665 78123 87079 86508 83086 122554 Emission pollution (kg) 147733 157076 138963 148091 145853 150296 136696 Minimization of operating costs and emission pollution Operating costs (MUs) 63976 62949 60941 69153 64970 66497 76671 Emission pollution (kg) 166019 163366 160151 169291 168154 166742 182969 Fig. 14. (a). The aggregated values of optimal day-ahead electrical, distributed energy resources of MGs for IBDRP1 case, (b) the aggregated values of optimal day- ahead heating distributed energy resources of MGs for IBDRP1 case, (c) the aggregated values of optimal day-ahead cooling distributed energy resources of MGs for IBDRP1 case. Fig. 15. (a). The aggregated values of optimal day-ahead electrical, distributed energy resources of MGs for IBDRP2 case, (b) the aggregated values of optimal day- ahead heating distributed energy resources of MGs for IBDRP2 case, (c) the aggregated values of optimal day-ahead cooling distributed energy resources of MGs for IBDRP2 case. Fig. 16. (a). The aggregated values of optimal day-ahead electrical, distributed energy resources of the energy system for IBDRP2 case, (b) the aggregated values of optimal day-ahead heating distributed energy resources of the energy system for IBDRP2 case, (c) the aggregated values of optimal day-ahead cooling distributed energy resources of the energy system for IBDRP2 case. which led to the curse of dimensionality. The scenario reduction method reduced the scenarios to 10 each. • The Monte Carlo simulation process was utilized to estimate the location and intensity of the external shocks [28]. For the Monte Carlo simulation, 300 trials were carried out to simulate the location and intensity of external shocks. For each external shock, the following contingencies were considered as the worst-case external shocks: 1) Triple energy carrier distribution network/pipeline out ages and single DER outage; 2) Triple DER outages; and 3) Combination of described outages. It was assumed that each external shock shut down the described facilities for the entire time of the simulation process and the corrective switching process and dis patching of DERs were performed. • For each external shock condition, the remaining energy system performance was optimized using the first level of the first stage optimization process. It was assumed that the external shock sec tionalized the energy system into on-outage and secured zones [24]. Based on the proposed method, the energy system operator should Fig. 17. (a). The optimal day-ahead dispatch of cooling energy storage of energy system for IBDRP2 case, (b) The optimal day-ahead dispatch of heating energy storage of energy system for IBDRP2 case, (c) The optimal day-ahead dispatch of cooling PHEV parking lots of energy system for IBDRP2 case. change the topology of electrical, heating, and cooling systems to restore the service for the on-outage zones using the DERs of neighbour-secured zones [24]. The detailed process of switching electrical switches, heating, and cooling valves is presented in [24] and is not presented for the sake of space. • The second level of the second-stage optimization process deter mined the optimal contribution of MGs in external shock conditions. The MGs’ DERs were dispatched for restoring the energy system using the switching process of electrical switches, and heating and cooling control valves. • The simulation was carried out on a PC (Intel Core i7–13,700 pro cessor, 128 GB memory, DDR4 3200 MT). The Monte Carlo simula tion time was about 3485 seconds. The detailed process of the Monte- Carlo simulation process is available in [29] and is not presented for the sack of space. The weighting factors are equal to one and the weighted sum method is carried out for the proposed objective functions. The objective functions and their constraints problems are linearized [24]. The method codes are developed in MATLAB and GAMS. Fig. 3 presents the flowchart of the method. As shown in Fig. 3, at the first level of the first stage optimization process, the energy system operator optimizes the day-ahead scheduling of the system’s distributed energy resources and estimates the values of MGs’ incentive-based and price-based demand response contributions. Then, the second level of the first stage optimization process determines the optimal scheduling of MGs’ DERs and incentive-based and price-based demand response var iables. It is assumed that based on the outputs of the second level of the first stage optimization problem, which update the data of MG’s DERs decision variables and incentive-based and price-based demand response variables, the first level of the first stage objective function can be recalculated and the day-ahead optimal scheduling of energy system can be determined. In the real-time scheduling process, if there is not any shock, the process simulates the shock conditions and the first and second levels of second-stage optimization processes determine the optimal dispatch of distributed energy resources of energy systems and MGs, respectively. The simulation of an external shock is performed to assess the behavior of the energy system in contingent conditions and find the optimal preventive scheduling of energy resources and available alternatives for system reconfiguration. This process should be carried out to reduce the impacts of external shock. Further, the available dispatching strategies of the energy system operator against probable shocks can be explored using the simulation process of external shock conditions. Based on the proposed method, in real-time operating conditions, if any shock condition is detected, the real-time optimization processes determine the optimal commitment of the energy system’s DERs, the system’s topology, and MGs’ optimal distributed energy generation commitments, which is one of the contributions of this paper. The optimization process minimizes the impacts of external shocks using the switching process of electrical switches and heating and cooling control valves. 4. Simulation results The simulation was carried out for the modified 123-bus IEEE sys tem. Fig. 4 depicts the system topology. Fig. 5 depicts the electrical, heating, and cooling loads of the energy system. As shown in Fig. 5, the aggregated energy of electrical, heating, and cooling loads are 209.8 MWh, 367.77 MWH, and 146.38 MWH, respectively. Fig. 6 presents the energy generation of wind turbines and photovoltaic systems for one of the reduced scenarios. The estimated values of energy generation of wind turbines and photovoltaic arrays are 30.88 MWh and 47.685 MWh, respectively. Fig. 7 and Fig. 8 show the forecasted price of day-ahead and real-time multi-carrier energy mar kets, respectively. The details of system data are presented in [24] and are not presented for the sake of space. The electrical, heating, and cooling demand of the microgrid was supplied by the microgrid multi-carrier distributed energy resources, Fig. 17. (continued). which transacts electricity, heating, and cooling energy with multi- carrier energy markets. The microgrids are categorized into residen tial, commercial, and industrial microgrids. The characteristics of distributed energy resources of microgrids are presented in Table 2. The MGs data are available in [25] and are not presented for the sake of space. Further, the Eq. (16), Eq. (17), and Eq. (25) data are available in Ref. [25]. The interruption costs of electricity, heating, and cooling loads are considered 0.42 MMUs/kWh [24]. MMUs stands for Million Monetary Units. Fig. 9. (a), (b), and (c) present the electrical, heating, and cooling loads of MGs for one of the reduced scenarios, respectively. Fig. 10 depicts the average values of the electrical power output of residential, commercial, and industrial MGs’ photovoltaic systems for different scenarios. As shown in Fig. 10, the average values of the electrical power output of residential, commercial, and industrial MGs’ photovoltaic systems are 7.195 kWh, 7.169 kWh, 7.00 kWh, 6.901 kWh, 7.190 kWh, 7.277 kWh, and 7.216 kWh, respectively. Further, the aggregated values of the electrical power output of residential, com mercial, and industrial MGs’ photovoltaic systems are 172.700 kWh, 172.121 kWh, 165.636 kWh, 172.573 kWh, 174.670 kWh, and 173.188 kWh, respectively. Fig. 11 shows the average values of the electrical power output of residential, commercial, and industrial MGs’ wind turbines for different scenarios. As shown in Fig. 11, the average values of the electrical power output of residential, commercial, and industrial MGs’ wind turbines are 4.499 kWh, 4.501 kWh, 4.310 kWh, 4.367 kWh, 4.461 kWh 7, 4.325 kWh, and 4.467 kWh, respectively. Further, the aggregated values of the electrical power output of residential, commercial, and industrial MGs’ photovoltaic systems are 107.978 kWh, 108.041 kWh, 103.450 kWh, 104.828 kWh, 107.071 kWh, 103.820 kWh, and 107.213 kWh, respectively. As presented in the problem formulations, different demand response programs are considered. Fig. 12 presents the value of fixed and time-of-use prices of the energy system for electrical, heating, and cooling loads. As shown in Fig. 12, it is assumed that the reward of electrical, heating, and cooling loads in the time-of-use prices of price-based de mand response process are equal to 120% of the average hourly prices of Fig. 18. (a). The day-ahead cost/benefit of energy system transactions with the multicarrier energy markets, and (b) the real-time cost/benefit of energy system transactions with the multicarrier energy markets. electrical, heating, and cooling markets, respectively. Further, the reward of electrical, heating, and cooling loads in the real-time prices of the price-based demand response process are equal to 110% of the real- time price of electrical, heating, and cooling markets, respectively. Fig. 13 depicts the value of incentive-based electrical, heating, and cooling demand response programs for MGs that participate in multi carrier energy markets based on their multi-carrier energy injections into energy system networks. As shown in Fig. 13, based on the volume of the energy injection of MGs into system energy, the incentive-based demand prices are multiplied by the corresponding real-time market. It is assumed that the penalty costs of MGs for any deviation of incentive-based demand response programs in electrical, heating, and cooling markets are equal to 150% of their real-time market prices. Table 3 depicts the feasible switching alternatives of Electric Switches and Control Valves (ES&CV) of district heating and cooling networks. As shown in Table 3, there are sixteen feasible topologies for the 123-bus system that the electrical, heating, and cooling loads can be supplied. The number of the external shock was forty cases and their impacts on the energy system were simulated. The worst-case scenario was considered as the following process. The first shock took place in energy corridor 160 at 15:00. The second shock came about in energy corridor 135 at 19:00. Twenty-one simulation cases were carried out to assess the impacts of different demand response alternatives on the optimal scheduling of the energy system and microgrids, which Table. 4 depicts these cases. Table 5 presents the average costs of residential, industrial, and commercial microgrids for the twenty-one case studies. As shown in Table 5, the real-time demand response integrated with the second package of incentive-based electrical, heating, and cooling demand response programs presented the minimum costs for microgrids considering simultaneous minimization of operating costs and emission pollution. The optimal incentive-based demand response program reduced the costs of residential microgrids by about 4.47% and 3.19% concerning the RTP and RTP (IBDRP1) cases, respectively. Further, the costs of industrial microgrids were reduced by about 3.88% and 0.57% concerning the RTP and RTP (IBDRP1) cases, respectively. Finally, the costs of commercial microgrids were reduced by about 5.47% and 1.44% concerning the RTP and RTP (IBDRP1) cases, respectively. Based on the Table. 5 results, the IBDRP2 and IBDRP1 alternatives were considered the first-best and second-best demand response program al ternatives, respectively. Based on the Table. 5 results, the IBDRP2 and IBDRP1 were Fig. 19. The resiliency index of the energy system for the worst-case shock conditions considering different demand response programs and without the pro posed method. Fig. 20. The optimal topology number of the energy system for different demand response alternatives and the corresponding worst-case shock. (c) Fig. 21. (a). The optimal real-time dispatch of the electrical energy system for IBDRP2 case and worst-case shock, (b) The optimal real-time dispatch of the heating energy system for IBDRP2 case and worst-case shock, (c) The optimal real-time dispatch of the cooling energy system for IBDRP2 case and worst-case shock. considered the first-best and second-best demand response alternatives, respectively. Thus, the MGs’ commitments to these cases are presented. Fig. 14 (a), (b), and (c) present the aggregated values of optimal day- ahead electrical, heating, and cooling distributed energy resources of MGs for the IBDRP1 case, respectively. Based on Fig. 14 (a), the aggregated electricity transactions of MGs with the energy system was 1494 kWh. Further, the aggregated electricity consumptions of MGs without and with IBDRP1 were 48,384 kWh and 48,384 kWh, respec tively. Thus, the IBDRP1 process only changed the electrical load profile and did not reduce the energy consumption of electrical loads. The aggregated electricity generations of non-intermittent and intermittent electricity generation facilities were 46,439 kWh and 3135 kWh, respectively. As shown in Fig. 14 (b), the aggregated heating energy transactions of MGs with the energy system were 537 kWh. Further, the aggregated heating energy consumptions of MGs without and with IBDRP1 were 19,353 kWh and 19,140 kWh, respectively. Thus, the IBDRP1 reduced the heating energy consumption by about 1.1%. The boilers produced heating energy without and with the IBDRP1 process by about 17,954 kWh and 17,805 kWh, respectively. Further, the ab sorption chillers consumed heating energy by about 10,120 kWh and 9648 kWh without and with the IBDRP1 process, respectively. Based on Fig. 14 (c), the aggregated cooling energy transactions of MGs with the energy system was 482 kWh. Furthermore, the aggregated cooling energy consumptions of MGs without and with IBDRP1 were 12,096 kWh and 11,912 kWh, respectively. It can be concluded that the IBDRP1 reduced the cooling energy consumption by about 1.51%. The compression chillers produced 1975 kWh and 2746 kWh of cooling energy without and with the IBDRP1 process, respectively. However, the cooling energy generation of absorption chillers was reduced by about 4.67% concerning without the IBDRP1 case. Fig. 15 (a), (b), and (c) depict the aggregated values of optimal day- ahead electrical, heating, and cooling distributed energy resources of MGs for the IBDRP2 case, respectively. As shown in Fig. 15 (a), the aggregated electricity transaction of MGs with the energy system was 1500 kWh. Thus, the IBDRP2 process increased the electricity trans actions of MGs with the energy system by about 0.43% concerning the IBDRP1 electricity transactions. The aggregated electricity consump tions of MGs without and with IBDRP2 were 48,384 kWh and 48,384 kWh, respectively. The IBDRP2 process did not reduce the energy con sumption of electrical loads. The aggregated electricity generations of non-intermittent and intermittent electricity generation facilities were 46,326 kWh and 3174 kWh, respectively. Thus, the IBDRP2 process reduced the energy generation of non-intermittent electricity generation facilities by about 0.24% and increased the commitment of energy storage-based intermittent electricity generation facilities by about 1.23% concerning the IBDRP1 case. As shown in Fig. 15 (b), the aggregated heating energy transaction of MGs with the energy system was 822 kWh. Thus, the IBDRP2 process increased the heating energy transactions of MGs with the energy system by about 53.02% concerning the IBDRP1 case. Furthermore, the aggregated heating energy consumptions of MGs without and with IBDRP2 were 19,353 kWh and 19,572 kWh, respectively. Thus, the IBDRP2 increased the heating energy consumption by about 1.12%, based on the fact that the proposed demand response tariff encouraged MGs to consume more heating energy at off-peak periods. The boilers produced heating energy without and with the IBDRP2 process by about 17,954 kWh and 18,522 kWh, respectively. Based on Fig. 15 (c) results, the aggregated cooling energy trans actions of MGs with the energy system was 684 kWh. Thus, the IBDRP2 process increased the cooling energy transactions of MGs with the en ergy system by about 42.08% concerning the IBDRP1 case. The aggre gated cooling energy consumptions of MGs without and with IBDRP2 were 12,096 kWh and 12,048 kWh, respectively. Thus, the IBDRP2 reduced the cooling energy consumption by about 0.39%. The compression chillers produced 1975 kWh and 3084 kWh of cooling energy without and with the IBDRP2 process, respectively. Thus, the cooling energy generation of compression chillers was increased by about 56.15% concerning without the IBDRP1 case. Fig. 16 (a), (b), and (c) present the aggregated values of optimal day- ahead electrical, heating, and cooling distributed energy resources of the energy system for the IBDRP2 case, respectively. The aggregated elec tricity transactions of MGs with the energy system before and after the IBDRP2 process were 177 MWh and 146 MWh, respectively, as shown in Fig. 16 (a). Thus, the IBDRP2 process decreased the electricity trans actions of the energy system with the wholesale electricity market by about 17.7% concerning the no-demand response case. The aggregated electricity consumptions of the energy system without and with IBDRP2 were 294 MWh and 286 MWh, respectively. The IBDRP2 process decreased the energy consumption of energy system electrical loads by about 2.75%. The aggregated electricity generations of non-intermittent electricity generation facilities were 170 MWh and 146 MWh, respec tively. Thus, the IBDRP2 mechanism reduced the energy generation of non-intermittent electricity generation facilities by about 14%. Based on Fig. 16 (b) results, the aggregated heating energy trans actions of the energy system with the heating energy market before and after the IBDRP2 process were 11.28 MWh and 5.58 MWh, respectively. Thus, the IBDRP2 process decreased the heating energy transactions of the energy system with the heating energy market by about 50.53%. Further, the aggregated heating energy consumptions of the energy system without and with the IBDRP2 process were 367.77 MWh and 348.03 MWh, respectively. Thus, the IBDRP2 decreased the heating energy consumption by about 5.53%. The boilers produced heating energy without and with the IBDRP2 process by about 183.35 MWh and 169.31 MWh, respectively. The aggregated cooling energy transactions of the energy system with the cooling energy market before and after the IBDRP2 process were 12.18 MWh and 9.33 MWh, respectively, as shown in Fig. 16 (c). It can be concluded that the IBDRP2 process decreased the cooling energy transactions of the energy system with the cooling energy market by about 23.39%. Furthermore, the aggregated cooling energy consump tions of the energy system without and with the IBDRP2 process were 146.38 MWh and 139.46 MWh, respectively. Thus, the IBDRP2 decreased the cooling energy consumption by about 4.72%. Fig. 17 (a), (b), and (c) present the optimal day-ahead dispatch of cooling energy storage, heating energy storage, and PHEV parking lots of the energy system for the IBDRP2 case, respectively. The maximum and minimum values of heating energy storage heating energy transactions with the energy system were 96 kWh and 32 kWh, respectively. Further, the average value of day-ahead heating energy transactions of heating energy storage facilities was 80 kWh. The maximum and minimum values of cooling energy storage cooling en ergy transactions with the energy system were 96 kWh and 30 kWh, respectively. Further, the average value of day-ahead cooling energy transactions of cooling energy storage facilities was 76 kWh. Finally, the maximum and minimum values of PHEV parking lots’ electricity transactions with the energy system were 444 kWh and 163 kWh, respectively. Further, the average value of day-ahead electricity trans actions of PHEV parking lots was 125 kWh. Fig. 18 (a) and Fig. 18 (b) present the day-ahead and real-time cost/ benefit of energy system transactions with the multicarrier energy markets for the IBDRP2 process, respectively. As shown in Fig. 18 (a), the energy system’s benefits for transactions with the day-ahead elec tricity, heating, and cooling markets were 0.635 MMUs, − 8910.1 MUs, and − 18512.3 MUs, respectively. Further, as shown in Fig. 18 (b) the energy system benefits for transactions with the real-time electricity, heating, and cooling markets were 1.27 MMUs, − 0.139 MMUs, and − 36715.2 MUs, respectively. Thus, the aggregated benefits of the energy system in the day-ahead and real-time markets were 0.61 MMUs and 1.10 MMUs, respectively; meanwhile, the energy system operational costs were about 0.49 MMUs and 0.63 MMUs in the day-ahead and real- time markets, respectively concerning no energy transactions with en ergy markets. Hence, the proposed method increased the benefits of the energy system using energy transactions with the day-ahead and real- time markets. Fig. 19 depicts the resiliency index of the energy system for the worst-case shock conditions considering different demand response programs and without the proposed method. As shown in Fig. 19, the IBDRP2 process provided a higher resiliency index. The resiliency index of the energy system tends to infinity using the proposed method for the worst-case shock scenario. Fig. 20 shows the optimal topology number of the energy system for the real-time optimization process for different demand response alter natives and the corresponding worst-case shock. Fig. 21 (a), (b), and (c) present the optimal real-time dispatch of electrical, heating, and cooling distributed energy systems of the energy system for IBDRP2, respectively considering the worst-case shock impact on the system. As mentioned earlier, the worst-case scenario occurred in energy corridor 160 at 15:00 and the second shock happened in energy corridor 135 at 19:00. The CHPs of the energy system continuously produced electricity by about 7.214 MW before and after external shocks, as shown in Fig. 21 (a). However, the electricity generation of DGs reduced from 9.17 MW to 7.47 MW for the first shock. Further, when the second shock occurred, the electricity generation of DGs increased from 4.98 MW to 5.51 MW for the second shock. The aggregated and average electricity consump tions of the system were 970.86 MWh and 10.12 MW, respectively for normal conditions. However, the aggregated and average electricity consumptions of the system were 928.30 MWh and 9.57 MW, respec tively. Thus, the IBDRP2 decreased the aggregated and average elec tricity consumption of the system by about 4.38% and 5.24%, respectively. Further, the aggregated and average electricity trans actions of the energy system with the wholesale electricity market were 724.65 MWh and 7.47 MW, respectively for normal conditions. Finally, the aggregated and average electricity transactions of the energy system with the wholesale electricity market were 758.72 MWh and 7.82 MW, respectively for shock conditions. Thus, the shock conditions increased the electricity transactions of the energy system with the wholesale electricity market by about 4.7%. Based on Fig. 21 (b) results, the aggregated and average heating energy consumptions of the system were 1243.73 MWh and 12.82 MW, respectively for normal conditions. However, the aggregated and average heating energy consumptions of the energy system were 1163.30 MWh and 11.99 MW, respectively. Thus, the proposed process decreased the aggregated and average heating energy consumption of the system by about 6.47% and 6.42%, respectively. Further, the aggregated and average heating energy transactions of the energy sys tem with the heating energy market were 31.67 MWh and 0.33 MW, respectively for normal conditions. Finally, the aggregated and average heating energy transactions of the energy system with the heating en ergy market were 66.71 MWh and 0.69 MW, respectively for shock conditions. Thus, the shock conditions increased the heating energy transactions of the energy system with the heating energy market by about 52.53%. The aggregated and average cooling energy consumptions of the system were 491.32 MWh and 5.06 MW, respectively for normal con ditions, based on Fig. 21 (c) results. The aggregated and average cooling energy consumptions of the system were 461.61 MWh and 4.75 MW, respectively. Thus, the proposed model decreased the aggregated and average cooling energy consumption of the system by about 6.04% and 6.12%, respectively. Further, the aggregated and average cooling energy transactions of the energy system with the cooling energy market were 344.64 MWh and 3.55 MW, respectively for normal conditions. Finally, the aggregated and average cooling energy transactions of the energy system with the cooling energy market were 287.49 MWh and 2.96 MW, respectively for shock conditions. Thus, the shock conditions decreased the cooling energy transactions of the energy system with the cooling energy market by about 16.58%. In conclusion, the proposed model successfully determines the optimum commitment of multi-carrier MGs, energy system DERs scheduling, and the topology of the system for normal and external shock conditions. 5. Conclusion This paper proposed a framework for day-ahead and real-time scheduling of microgrids and the energy distribution system that transacted multi-carrier energy with each other in normal and shock conditions. The energy system also transacted electrical, heating, and cooling energy carriers with wholesale electricity market, heating en- ergy market, and cooling markets, respectively. The proposed method calculated the resiliency index of the energy system and determined the optimal status of electrical switches, heating, and cooling control valves to enhance the resiliency of the system. The two-stage two-level optimization model determines the scheduling of the energy system and microgrids considering the incentive-based and price-based demand response procedures. The proposed method was simulated for the 123-bus test system. The proposed incentive-based demand response process successfully reduced the aggregated and average electricity consumption of the system by about 4.38% and 5.24%, respectively. Further, the proposed process decreased the aggregated heating and cooling energy consumption of the system by about 6.47% and 6.04%, respectively. The authors are working on the modeling of a hydrogen energy carrier system to integrate with the proposed model. References [1] Jiang Tao, Dong Xinru, Zhang Rufeng, Li Xue. Strategic active and reactive power scheduling of integrated community energy systems in day-ahead distribution electricity market. Appl Energy 2023;336:120558. [2] Xiao Min, Smaisim Ghassan Fadhil. Joint chance-constrained multi-objective optimal function of multi-energy microgrid containing energy storages and carbon recycling system. J Energy Stor 2022;55:105842. [3] Pang Kang Ying, Liew Peng Yen, Woon Kok Sin, Ho Wai Shin, Alwi Sharifah Rafidah Wan, Klemes Jirí Jaromír. Multi-period multi-objective optimisation model for multi-energy urban-industrial symbiosis with heat, cooling, power and hydrogen demands. Energy 2023;262:125201. [4] Li Yaowang, Yao Fuxing, Zhang Shixu, Liu Yuliang, Miao Shihong. An optimal dispatch model of adiabatic compressed air energy storage system considering its temperature dynamic behavior for combined cooling, heating and power microgrid dispatch. J Energy Stor 2022;51:104366. [5] Zhang Zhiyang, Altalbawy Farag MA, Al-Bahrani Mohammed, Riadi Yassine. Regret-based multi-objective optimization of carbon capture facility in CHP-based microgrid with carbon dioxide cycling. J Clean Prod 2023;384:135632. [6] Mianaei Peyman Khorshidian, Aliahmadi Mohammad, Faghri Safura, Ensaf Mohammad, Ghasemi Amir, Abdoos Ali Akbar. Chance-constrained http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0005 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0005 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0005 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0010 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0010 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0010 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0015 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0015 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0015 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0015 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0020 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0020 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0020 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0020 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0025 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0025 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0025 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0030 http://refhub.elsevier.com/S0306-2619(23)01083-8/rf0030 programming for optimal scheduling of combined cooling, heating, and power- based microgrid coupled with flexible technologies. Sustain Cities Soc 2022;77: 103502. [7] Vahedipour-Dahraie Mostafa, Rashidizadeh-Kermani Homa, Anvari- Moghaddam Amjad, Siano Pierluigi, Catalao Joao PS. Short-term reliability and economic evaluation of resilient microgrids under incentive-based demand response programs. Electr Power Energy Syst 2022;138:107918. [8] Guo Jiqun, Tan Jinjing, Li Yang, Haifei Gu, Liu Xin, Cao Yang, et al. Decentralized incentive-based multi-energy trading mechanism for CCHP-based MG cluster. Electr Power Energy Syst 2021;133:107138. [9] Guo Qun, Nojavan Sayyad, Lei Shi, Liang Xiaodan. Economic-environmental evaluation of industrial energy parks integrated with CCHP units under a hybrid IGDT-stochastic optimization approach. J Clean Prod 2021;317:128364. [10] Lekvan Amir Aris, Habibifar Reza, Moradi Mehran, Khoshjahan Mohammad, Nojavan Sayyad, Jermsittiparsert Kittisak. Robust optimization of renewable-based multi-energy micro-grid integrated with flexible energy conversion and storage devices. Sustain Cities Soc 2021;64:102532. [11] Li Ke, Wei Xingguo, Yan Yi, Zhang Chenghui. Bi-level optimization design strategy for compressed air energy storage of a combined cooling, heating, and power system. J Energy Stor 2020;31:101642. [12] Saberi Kasra, Pashaei-Didani Hamed, Nourollahi Ramin, Zare Kazem. Optimal performance of CCHP based microgrid considering environmental issue in the presence of real time demand response. Sustain Cities Soc 2019;45:596–606. [13] Zhou Xiaoqian, Ai Qian. Distributed economic and environmental dispatch in two kinds of CCHP microgrid clusters. Electr Power Energy Syst 2019;112:109–26. [14] Tooryan Fatemeh, HassanzadehFard Hamid, Dargahi Vahid, Jin Shuangshuang. A cost-effective approach for optimal energy management of a hybrid CCHP microgrid with different hydro