Optimal Operation for Economic and Exergetic Objectives of a Multiple Energy Carrier System Considering Demand Response Program

Abstract

An MECS (multiple energy carrier system) could meet diverse energy needs owing to the integration of different energy carriers, while the distinction of quality of different energy resources should be taken into account during the operation stage, in addition the economic principle. Hence, in this paper, the concept of exergy is adopted to evaluate each energy carrier, and an economic–exergetic optimal scheduling model is formulated into a mixed integer linear programming (MILP) problem with the implementation of a real-time pricing (RTP)-based demand response (DR) program. Moreover, a multi-objective (MO) operation strategy is applied to this scheduling model, which is divided into two parts. First, the ε-constraint method is employed to cope with the MILP problem to obtain the Pareto front by using the state-of-the-art CPLEX solver under the General Algebraic Modeling System (GAMS) environment. Then, a preferred solution selection strategy is introduced to make a trade-off between the economic and exergetic objectives. A test system is investigated on a typical summer day, and the optimal dispatch results are compared to validate the effectiveness of the proposed model and MO operation strategy with and without DR. It is concluded that the MECS operator could more rationally allocate different energy carriers and decrease energy cost and exergy input simultaneously with the consideration of the DR scheme.

1. Introduction

At present, to realize the sustainable development of society, many scholars have been paying attention to energy savings and efficiency improvements due to severe fossil energy crises and environmental problems [1]. In such context, an MECS (multiple energy carrier system), which was formulated in a project “Vision of Future Energy Networks” at ETH Zurich together with other partners, has been regarded as a promising approach for next-generation energy systems [2]. An MECS refers to a system coupled with diverse kinds of energy carriers, such as electricity, gas, heat, and cooling to achieve coordinate energy production and delivery, as well as to meet various energy consumptions in coordination [3].

Modeling and optimal operation are the two key issues in achieving the expected potentials of an MECS. The basic concept of an energy hub (EH) has been introduced from a global perspective to build up the model of the MECS [4], where the core of the EH is a coupling matrix description. A standardized matrix modeling method was proposed, and the coupling matrix was automatically built based on graph theory and the Gaussian elimination technique, wboth hich could be conveniently implemented in a computer to facilitate the modeling effort [5]. Another standardized modeling approach for a complex EH model was introduced in [6]. The coupling matrix of complex EH model was divided into several simple EH models based on nodal arrangement and virtual node insertion methods. Though most researchers have adopted the EH concept to model an MECS [7,8], some scholars have formulated the model from local perspectives. The energy conversion and coupling relationship are separately described for each device and energy carrier, with obvious an distinction compared to the EH concept, which is utilized to establish the relationship of the whole system from a global perspective. The device models for photovoltaic (PV), electrical energy storage (ESS), combined heat and power (CHP) plant, gas boiler (GB), heat pump (HP), and thermal energy storage (TES) have been developed with input–output characteristics, capacity constraints, ramp rate limits and other conditions [9]. In addition, models for each energy carrier (electricity, gas, and heat) have been proposed based on the energy balance principle in the MECS [9,10]. Based on MECS modeling, many scholars have focused on optimal operation strategy to realize rational energy distribution and efficient energy management. For minimizing the day-ahead operation cost, a mixed integer nonlinear programming (MINLP) approach was employed to solve the optimal scheduling problem of building an energy system, which is coupled with electricity and cooling under dynamic electricity pricing [11]. In [12], operation strategies for both grid-connected and islanded modes of a multi-energy microgrid were implemented by the CPLEX solver in General Algebraic Modeling System (GAMS) with the aim of minimizing net operating cost. A multi-objective (MO) optimal electrical and thermal energy management for a residential EH was conducted to acquire the Pareto front when considering the environmental emissions and total energy cost, and the fuzzy method was employed to determine the best possible solution on the Pareto front [13]. However, most researchers have focused on the EH model and the single-objective optimal economic operation in the MECS, in which the local relation of each component might not be adequately considered; additionally the operation goal neglects the characteristics of different energy carriers.

To realize the rational utilization of energy in an MECS, the energy quality should be taken into account in the operation strategy, and the measurement of energy quality for different energy carriers needs to be conducted. Exergy, derived from the second law of thermodynamics, has been employed as a promising tool to analyze energy quality. It is the maximum theoretical work obtainable from the system as it reaches the equilibrium state with its surrounding environment, the so-called “reference environment” [14,15]. Exergy could be regarded as the usable part of energy, and the other part which cannot be utilized is called anergy [16]. Some scholars have conducted considerable research about the exergy issue. An exergetic analysis and a performance evaluation were conducted on some renewable energy resources, including solar energy systems, wind energy systems, geothermal energy systems, and biomass systems [17]. Energy, exergy and environmental analyses of a hybrid combined cooling heating and power (CCHP) system driven by biomass and solar energy have been studied, and the complementarity of the two kinds of energy has been analyzed [18]. To determine the best design parameters of a polygeneration system which was composed of gas turbine cycle, a Rankine cycle, an absorption cooling cycle, and domestic hot water, a multi-objective optimization approach was applied based on a genetic algorithm considering the total cost and exergy efficiency as objective functions [19]. Most studies have concentrated on the system design and performance evaluation of different energy systems using the exergy method, but few researchers have contributed to the optimal operation strategy of the MECS. In [20], the concept of exergy hub, which is analogous to EH, was proposed, and three single-objective optimal dispatch cases of MECSs were investigated for exergetic efficiency maximization, energetic efficiency maximization, and cost minimization. A multi-objective mixed integer linear programming (MOMILP) problem for the operation optimization of a distributed energy system was solved, and the weighted sum approach was employed to convert the MOMILP problem to a single-objective MILP according to [21]. However, these studies for optimal operation lack flexibility, and the interaction between the MECS and load demand side are inadequately considered.

Demand response (DR) could contribute to the rational use of electricity by realizing the peak load shifting and curtailment, benefitting both grid utility and customers. Generally, DR is divided into two categories, price-based DR (PBDR) programs and incentive-based DR (IBDR) programs [22]. In PBDR programs, dynamic electricity pricing is employed to decrease the peak load and increase the valley load, where three major cases are developed—time of use (TOU), critical peak pricing (CPP) and real-time pricing (RTP). On the other hand, direct load control (DLC), emergency demand response program (EDRP), capacity market program (CAP) and other terms are included in IBDR, where participants could receive participation payments [23]. Among these approaches, the RTP program might be the most direct and efficient DR program, a statement supported by by many economists [24]. Many works in the literature about the RTP-based DR program have investigated power systems to measure the economic benefit with a deterministic or stochastic price. A smart home energy system (HES) structure was proposed with a deterministic time-varying price signal in [25]. To realize the optimal energy management for residential appliances, the Monte Carlo simulation and robust optimization technique were adopted to address the uncertainties of RTP in PBDR in [26]. Because of the successful implementation of DR in the power system, some farsighted scholars introduced it into an MECS to study their interaction. Three scenarios were conducted—PBDR programs (TOU, CPP, and TOU and CPP), IBDR programs (EDRP and CAP) and a combination between the PBDR and IBDR programs (TOU and CAP)—to simultaneously minimize EH operation cost and maximize customer benefit, according to [27]. In [28], an RTP model was developed that was calculated by the predicted electricity demand and TOU pricing, and it was implemented in a micro energy grid which was described by the EH concept. However, DR programs are not adequately considered in the operation strategy with the exergetic objective, and the interaction between DR and multi-objectives has been less investigated by most research.

In this paper, the main goal was to realize rational energy distribution with the principles of simultaneously decreasing the operation cost and minimizing the exergy input for an MECS. A grid-connected MECS was built to meet the diverse energy needs of loads. To obtain an optimal energy management strategy, an optimal MO scheduling model considering both the economic and exergetic objectives was introduced from the local perspectives with the implementation of the RTP-based DR scheme, which could flatten the electricity load curve. The Pareto front was obtained by the ε-constraint approach to the optimal MO scheduling model, and the compromise solution on the Pareto front was acquired by the approach of linear programming technique for a multidimensional analysis of preference (LINMAP). Therefore, the novelty and contribution of this paper could be summarized as follows:

• The MO optimal scheduling model of an MECS with economic and exergetic objectives are conducted and formulated into an MILP problem.
• To evaluate the performance of the DR program in an MECS, an RTP-based DR model is introduced with a dynamic pricing rate.
• An MO operation strategy is developed where the ε-constraint method is utilized to obtain the Pareto front, and the LINMAP approach is introduced to make a trade-off between the economic and exergetic objectives.

In the following, the economic–exergetic scheduling model of an MECS is introduced, and the RTP-based DR scheme is formulated in Section 2. In Section 3, the MO operation strategy is presented to obtain the Pareto front, and a trade-off is conducted between these two different objectives. Case studies are conducted in Section 4, and the conclusion for this paper is summarized in Section 5.

2. Problem Formulation

The energy generation, transmission, coupling and consumption of different energy carriers such as coal, solar, natural gas, electricity, cooling, and heat are described in Figure 1 to fulfill the diverse energy demands of this paper. Traditionally, primary energy resources have been utilized through energy generation and energy transmission, and the energy coupling process have not really been taken into account in the energy supply chain. Therefore, the energy coupling segment, which is regarded as an MECS, is focused on in this paper. Some energy conversion devices are considered to realize the transformation between different energy carriers, such as a photovoltaic (PV), a solar thermal collector (STC), a combined heat and power (CHP), a gas boiler (GB), an electrical chiller (EC), an absorption chiller (AC) and a heat exchanger unit (HEU). In addition, energy storage systems (ESSs), such as cooling energy storage (CES) and thermal energy storage (TES), are well considered to boost the flexibility of each energy carrier. The MO optimal operation model of a grid-connected MECS is formulated in this section, which is divided into three parts. In Section 2.1, the economic and exergetic objective models are formulated. The operation model for each energy device, which contains equality constraints and inequality constraints, is developed in Section 2.2. To measure the effect of DR for an MECS, the RTP model for the DR scheme is introduced in Section 2.3.

2.1. Objectives

2.1.1. Economic Objective

In this study, the electricity was not permitted to be sold back to the utility grid, which meant that bi-directional power flow was not considered. The energy cost of an MECS is the aggregation of two terms: the cost of electricity bought from the grid and the cost of natural gas bought from the natural gas pipeline; this relationship is described as follows:

_f1=∑t=124(_ce(t)_Pgride(t)+_cg(t)LH_{Vgas Gg}(t)/3.6)Δt

where f₁ is the economic objective,_ce(t)and_cg(t)are the electricity and gas pricing at period t,_Pgride(t)is the power purchased from the utility grid,_Gg(t)is the natural gas bought from the gas pipeline, andLH_Vgasis the lower heat value (LHV) of natural gas.

2.1.2. Exergetic Objective

Due to the coupling of different energy carriers, energy efficiency might not be the best evaluation index, since it emphasizes energy quantities but neglects energy qualities. To this end, maximizing exergy efficiency of the whole day is a better choice for an evaluation index in the day-ahead operation of an MECS, which is formulated as:

_ψex=∑t=124E_xout(t)Δt/∑t=124E_xin(t)Δt

where_ψexis the exergy efficiency andE_xout(t)andE_xin(t)are the exergy output and input of an MECS, respectively.

According to Figure 1, the exergy used is equal to the exergy needed from the energy consumption segments and is presented as:

E_xout(t)=E_xoute(t)+E_xoutc(t)+E_xouth(t)

whereE_xoute(t),E_xoutc(t), andE_xouth(t)are the electricity, cooling, and heat exergy output of an MECS, respetively.

Electricity energy is pure exergy according to [16], because it can be completely converted to work. Cooling and heat exergy are tightly related to the required temperature of demands and reference temperature, which generally adopt ambient temperature. Energy flowing with different temperatures means different abilities to work; hence, the temperature factor must be taken into account to formulate cooling and heat exergy. The exergy formulation for each energy carrier is introduced as follows:

E_xoute(t)=_Le(t)

E_xoutc(t)=(_Ta(t)/_Treqc−1)_Lc(t)

E_xouth(t)=(1−_Ta(t)/_Treqh)_Lh(t)

where_Le(t),_Lc(t), and_Lh(t)are the load demands of electricity, cooling, and heating at period t, respectively;_Ta(t)is the ambient temperature;_Treqcand_Treqhare the required temperatures of the cooling and heat loads, respectively.

However, it must be noted that diverse energy demands could be obtained by load forecast techniques with high precision, and total exergy needs could be calculated by forecasted load profiles. Therefore, the exergetic objective could be converted from maximizing exergy efficiency to minimizing the exergy input of the MECS according to Equation (2), which is formulated as follows:

_f2=∑t=124E_xinco(t)+E_xins(t)+E_xing(t)

where_f2is the exergetic objective;E_xinco(t),E_xins(t), andE_xing(t)are the exergy input of each primary energy carriers of coal, solar and gas, respectively.

In this paper, the electricity on the grid was assumed to be generated from the power plant, the energy of which was derived from coal, as shown in Figure 1. The exergy in solar energy has been studied in many papers, including various approaches reported in [29,30]. Among these, in this paper, the exergy-to-energy ratio of radiation developed by Petela was adopted to describe solar exergy, which is regarded as thermal radiation at the sun temperature [15,17]. Chemical-specific exergy is often evaluated by the LHV and exergy factor of the energy carrier [20], and the exergy of natural gas could be formulated by specific exergies and mass flow rates. Therefore, the exergy rate for each primary energy carrier could be expressed as:

{E_xinco(t)=E_xine(t)/_ψexpE_xine(t)=_Pgride(t)

{E_xins(t)=I(t)(_SSTC+_SPV)γ(t)/1000γ(t)=1+13^{(_Ta(t)/_Tsun)4}−43(_Ta(t)/_Tsun)

{_πg=LH_{Vgas ξgas}E_xing(t)=_{πg Gg}(t)/3.6

where_ψexpis the exergy efficiency of the power plant;I(t)is the solar radiation;_SSTCand_SPVare the installed area of the STC and the PV;γ(t)is the exergy-to-energy ratio of solar energy;_Tsunis the temperature of the Sun;_πgand_ξgasare the specific exergy and exergy factor of natural gas, respectively.

2.2. Constraints

The establishment of proper constraints is a vital issue to rationally allocate different energy flows on the premise of fulfilling diverse energy needs and guaranteeing the security of energy devices in an MECS. Therefore, relative equality and inequality constraints must be taken into account for each energy carrier and energy device, as shown in Figure 1. Balancing electricity, heat, cooling and natural gas should be considered to meet energy demands and conduct energy distribution between each energy device in the MECS. For every energy equipment, balance and restrictions must be simultaneously involved. Moreover, for the sake of intuitively reflecting the coupling relationship between energy equipment and energy flow, the modeling approach was employed from local perspectives in this paper.

2.2.1. Energy Flows Constraints

Balancing electricity, heat, cooling and natural gas should be considered to meet the energy demands and conduct energy distribution between each energy device in an MECS. The electricity needs of loads and the EC must be jointly satisfied by a PV, a CHP, and a grid at each time interval. Heat and cooling loads are satisfied by an HEU, chilling devices, and ESS, which could adjust the balance of supply and demand and increase the flexibility of the MECS. In addition, in the MECS, the heat flow distribution among the STC, the CHP, the GB, the AC, and the HEU should be considered. Since there is no natural gas need in the energy demands, only the allocation relationship between the CHP and the GB just be taken into account. Therefore, the energy balance model for each energy carrier is addressed as follows:

_Le(t)+_PECe(t)=_Pgride(t)+_PPVe(t)+_PCHPe(t)

_Lh(t)=_PHEUh,out(t)+_PTESh,dch(t)−_PTESh,ch(t)

_Lc(t)=_PECc(t)+_PACc(t)+_PCESc,dch(t)−_PCESc,ch(t)

_PACh(t)+_PHEUh,in(t)=_PSTCh(t)+_PCHPh(t)+_PGBh(t)

_Gg(t)=_GGB(t)+_GCHP(t)

where_PECe(t)is the electricity need of the EC at period t;_PPVe(t)and_PCHPe(t)are the electricity generation of the PV and the CHP;_PHEUh,in(t)and_PHEUh,outare the heat input and output of the HEU;_PTESh,dch(t)and_PTESh,ch(t)are the heat discharg and charge of TES;_PECc(t)and_PACc(t)are the cooling generated by the EC and the AC;_PCESc,dch(t)and_PCESc,ch(t)are the cooling discharged and charged of CES;_PACh(t)is the heat consumption of the AC;_PSTCh(t),_PCHPh(t)and_PGBh(t)are the heat generated by the STC, the CHP and the GB, respectively; and_GGB(t)and_GCHP(t)are the gas demands of the GB and the CHP, respectively.

2.2.2. Energy Device Constraints

In the MECS shown in Figure 1, the PV and the STC are regarded as uncontrollable devices, and their energy output rates depend on solar irradiance and other conditions, such as the conversion efficiency of the device and time-varying ambient temperature. Therefore, according to [13,15], the electricity and heat output rates of the PV and the STC could be formulated as:

_PPVe(t)=_{ηPV SPV}I(t)(1−0.005(_Ta(t)−298.15))/1000

_PSTCh(t)=_{ηSTC SSTC}I(t)/1000

where_ηPVand_ηSTCare the energy conversion efficiency of the PV and the STC, respectively.

The CHP unit is employed to realize the coupling of electricity, heat, and natural gas, which is vital in the optimal operation of an MECS due to its controllability. Hence, the capacity and ramp limit issues must be taken into account in the operating model. The constant efficiency model was adopted to describe the coupling relation of different energy carriers according to [13,14], and the 0–1 variable was introduced to describe the start/stop state because the electricity rate of the CHP could not be decreased to zero due to the technique restrictions, as shown in the following:

{_PCHPe(t)=_{ηCHPe GCHP}(t)LH_Vgas/3.6_PCHPh(t)=_{ηCHPh GCHP}(t)LH_Vgas/3.6_uCHP(t)_PCHP,mine≤_PCHPe(t)≤_uCHP(t)_PCHP,maxeR_PCHP,mine≤_PCHPe(t+Δt)−_PCHPe(t)≤R_PCHP,maxe

where_ηCHPeand_ηCHPhare the electricity and heat efficiency of the CHP, respectively;_PCHP,mineand_PCHP,maxeare the lower and upper limits of capacity, respectively;_uCHP(t)is a binary variable; 1/0 denotes the start/stop state;R_PCHP,mineandR_PCHP,maxerepresent the ramp-down and ramp-up bounds, respectively.

The GB consumes natural gas to generate heat which could be utilized in the AC and the HEU to meet the cooling and heat demands, and the energy balance and capacity restrictions should be considered:

{_PGBh(t)=_{ηGB GGB}(t)LH_Vgas/3.6_uGB(t)_PGB,minh≤_PGBh(t)≤_uGB(t)_PGB,maxh

where_ηGBis the heat efficiency of the GB;_uGB(t)represents the start/stop state;_PGB,minhand_PGB,maxhare, respectively, the lower and upper limits of the GB.

High-temperature heat flow from the STC, the CHP, and the GB need to be exchanged with low-temperature heat in the HEU according to the required temperature of the user. In this paper, it was assumed that the HEU had enough capacity to consume heat to avoid heat waste. Due to heat loss in the heat exchange process, the efficiency of the HEU was developed to establish the conversion model, i.e.:

_PHEUh,out(t)=_{ηHEU PHEUh,in}(t)

where_ηHEUis the efficiency of the HEU.

To fulfill the cooling demands, the EC and AC devices are included in this paper by consuming different energy carriers. The EC needs electricity, which is a kind of high-quality energy to drive, and the AC takes advantage of low-quality heat to generate cooling energy; their energy conversion relationship can be formulated by the coefficient of performance (COP), which is regarded as a constant in this paper, as follows:

{_PECc(t)=CO_{PEC PECe}(t)_uEC(t)_PEC,minc≤_PECc(t)≤_uEC(t)_PEC,maxc

{_PACc(t)=CO_{PAC PACh}(t)_uAC(t)_PAC,minc≤_PACc(t)≤_uAC(t)_PAC,maxc

whereCO_PECandCO_PACare the COP of the EC and the AC, respectively;_PEC,mincand_PEC,maxcare the lower and upper bounds of the EC;_PAC,mincand_PAC,maxcare the lower and upper bounds of the AC; and_uEC(t)and_uAC(t)denote the on/off states of the EC and the AC, respectively, at period t.

To further improve the flexibility of the MECS, TES and CES were introduced to achieve the energy shift during the day with three modes of operation: charging, discharging and still. The general model for ESS was adopted. The energy stored in ESS at time interval t depends on the energy stored at the previous time step as well as the energy variation, both of which directly relate to the mode selection. The energy dissipation rate and charging/discharging efficiency were applied to the ESS model. Furthermore, the capacity of ESS and the charging/discharging restrictions must be taken into account. Therefore, the models for TES and CES were concluded as:

{_ETES(t)=(1−_δTES)_ETES(t−Δt)+(_{ηTESch PTESch}(t)−_PTESdch(t)/_ηTESdch)Δt_ETESmin≤_ETES(t)≤_{ETESmaxuTESch}(t)_PTES,minch≤_PTESch(t)≤_uTESch(t)_{PTES,maxchuTESdch}(t)_PTES,mindch≤_PTESdch(t)≤_uTESdch(t)_{PTES,maxdchuTESch}(t)+_uTESdch(t)≤1

{_ECES(t)=(1−_δCES)_ECES(t−Δt)+(_{ηCESch PCESch}(t)−_PCESdch(t)/_ηCESdch)Δt_ECESmin≤_ECES(t)≤_{ECESmaxuCESch}(t)_PCES,minch≤_PCESch(t)≤_uCESch(t)_{PCES,maxchuCESdch}(t)_PCES,mindch≤_PCESdch(t)≤_uCESdch(t)_{PCES,maxdchuCESch}(t)+_uCESdch(t)≤1

where_ETES(t)and_ECES(t)are the energy stored in TES and CES at time t, respectively;_δTESand_δCESare the energy dissipation rates;_ηTESchand_ηTESdchare the charging and discharging efficiency of TES, respectively;_ηCESchand_ηCESdchare the charging and discharging efficiencies of CES, respectively;_ETESminand_ETESmaxare the lower and upper limits of TES, respectively;_ECESminand_ECESmaxare the lower and upper limits of CES, respectively;_PTES,minch,_PTES,maxch,_PCES,minch, and_PCES,maxchare the charging bounds for TES and CES;_PTES,mindch,_PTES,maxdch,_PCES,mindch, and_PCES,maxdchare the discharging bounds for TES and CES;_uTESch(t),_uTESdch(t),_uCESch(t)and_uCESdch(t)are all binary variables which denote the operation state for TES and CES.

2.3. RTP-Based DR Model

As described in Section 1, the DR program could achieve the interaction between the demand side and the MECS, which could benefit from a decreasing electricity cost and avoidance of capacity investment. The RTP-based DR program as one kind of PBDR is employed in this section. The RTP modeling is introduced by the day-ahead forecasted electricity load and TOU pricing to acquire the dynamic RTP tariff [28]. Furthermore, the electricity load model was conducted based on the RTP tariff.

To acquire the day-ahead RTP tariff, the day-ahead average electricity demand needs to be calculated first, and then the floating factor of RTP could be obtained. At last, the RTP tariff could developed by TOU price rate and floating factor, as shown in the following:

_Leav=∑t=124_Le(t)/24

θ(t)=_Le(t)/_Leav

_cRTP(t)=_cTOU(t)θ(t)

In addition, the upper and lower limit of the RTP tariff should be considered, i.e.:

_cRTP,min≤_cRTP(t)≤_cRTP,max

Based on the RTP tariff, the new load model could be formulated based on self-elasticity coefficient and cross-elasticity coefficient. The self-elasticity factor is always negative, while the cross-elasticity factor is always positive. By evaluating the ability of the RTP-based DR, the new load profile could be acquired according to [23,31], which is shown as follows:

_LeRTP(t)=_Le(t){1+e(t,t)_cRTP(t)−_cTOU(t)_cTOU(t)+∑k=1k≠t24e(t,k)_cRTP(k)−_cTOU(k)_cTOU(k)}

where_Leavis the average electricity load;θ(t)is the floating factor;_cRTP(t)and_cTOU(t)are the electricity price based on the RTP and TOU schemes, respectively;_cRTP,minand_cRTP,maxare the lower and upper bounds of the RTP tariff, respectively;_LeRTP(t)is the new load model based on RTP pricing; e(t,t) is the self-elasticity coefficient; and e(t,k) represents the cross-elasticity coefficient.

3. MO Operation Strategy 3.1. MO Optimization Approach

Based on the optimal scheduling model formulated in Section 2, objectives (Equations. (1) and Equations. (7)) and constraints (Equations (11)–(24)) were all MILP modelled. Therefore, some heuristic evolutionary algorithms were not suitable for this paper’s model, which might suffer from computational demands and inconsistent solutions owing to the high-dimension and dynamic variables in the MECS [12]. In this paper, the ε-constraint method, as a mathematical programming technique, was implemented to cope with the MO optimization problem, as it converts the MO model into a series of single-objective models. In this approach, one of the objective functions was regarded as an inequality constraint, and the other one was optimized, as shown below:

min_f2s.t.{_f1≤εEqs. (11)−(24)

where ε varies from f_1,max to f_1,min by setting the values of i and n according to:

ε=_f1,max−(_f1,max−_f1,min)(i−1)/(_imax−1),i=1,2…_imax

By solving the single-objective optimization problem by minimizing f₁ with Equations (11)–(24), the optimal solution X₁ could be acquired, and then the values of f₁ and f₂ could be obtained, which are indicated by f_1,min and f_2,max. In addition, X₂, f_2,min and f_1,max could be obtained by minimizing f₂ with Equations (11)–(24). Finally, the ε-constraint approach could be applied to the optimal scheduling model, and the Pareto front could then be acquired. 3.2. Preferred Solution Selection Approach It must be noted that the points on the Pareto front are all optimal values, and their main differences are the preferences from the operator. Therefore, in this study, the LINMAP is introduced to assist the operator to make a trade-off among the economic and exergetic objectives.

According to [32], the concepts of the Utopia point and the Nodir point can be introduced to conduct the most preferred solution among Pareto set points. The Utopia point corresponds to the fact that all objectives are at their best values f_k,min, and the Nodir point is a specific point where all objectives are at their worst values f_k,max. Hence, the distance between each Pareto point and Utopia point is employed in the LINMAP, which means that a shorter distance represents a better solution. However, due to the distinct magnitude of the two objective functions, normalization is required before the implementation of this principle, and then the distance between each normalized point and (0,0) can be calculated, as shown below:

_ρkj={1_fkj≥_{fk,max fkj}−_fk,minfk,max−_fk,minfk,min≤_fkj≤_fk,max0_fkj≤_fk,min

_dj=∑k^(_ρkj)2

where_fk,minand_fk,maxare the minimum and maximum value for k-th objective, respectively;_fkjand_ρkjare the optimal value and optimal normalized value of j-th Pareto point for k-th objective, respectively; and_djrepresents the evaluation index of j-th Pareto point.

3.3. MO Operation Rules

Based on the solution for MO optimal scheduling model in Section 3.1 and Section 3.2, the flowchart for the MO optimal operation of the MECS is developed in this part, as shown in the following (Figure 2):

4. Case Study In order to demonstrate the rationality of the economic–exergetic scheduling model with a MILP description and to show the effectiveness of the MO operation strategy, case studies were simulated for this section. Two cases (with and without RTP-based DR) are investigated to highlight the merit of the DR program in the MECS, and the results are discussed in detail to compare the distinction of the scheduling results. In addition, sensitivity experiments were conducted to measure the objective variation caused by diverse implementation strengths of the DR scheme. Furthermore, due to the MILP formulation of the optimal scheduling model of the MECS, the state-of-the-art CPLEX solver was employed to address the MILP issue under the GAMS environment. 4.1. System Description

The energy supply chain of the MECS (described in Figure 1) fulfills the secondary energy carrier needs (electricity load, cooling load, and hot water load) from the energy user side by the coupling of different energy carriers and the interaction of diverse energy conversion units and energy storage devices. However, it must be noted that, in this paper, the MECS is not an energy coupling system in reality. We just demonstrate the optimal economic–exergetic scheduling model and the MO operation strategy by way of simulation. The detailed energy topology and energy device model were discussed in Section 2; therefore, in this subsection, we mainly focus on the energy device configuration and data sets for the optimal energy management of the MECS. Some literature have conducted optimal operation of an MECS for different seasonal days, such as in [12,14]. However, in this paper, data required for a typical day in summer were implemented as a test example to verify the rationality and effectiveness of the optimal scheduling model and solving algorithm, both og which are not only suitable for summers day but also for each seasonal day.

According to [28], the demands of cooling, heat, and electricity without the consideration of the DR program are given in Figure 3. The electricity loads with the DR scheme implementation for each interval which were calculated by the RTP-based DR model with a self-elasticity coefficient -0.2 and a cross-elasticity coefficient 0.01, according to [31], and are presented in Figure 3a. There are some similar and distinct features by the comparison of these two electricity loads curves. Two peaks and two valleys exist for both cases at the same periods in the summer typical day, while the DR implementation could smooth the electricity demand curve to realize peak shaving and valley filling. Moreover, the cooling and heat energy are required to meet the needs of space cooling and domestic hot water and the energy demands for each interval are illustrated in Figure 3b. In addition, to evaluate the exergy efficiency of the MECS, the temperature parameters of cooling and heat energy were considered to be 299.15 and 333.15 K, respectively [14].

To evaluate the distinction between exergy efficiency and energy efficiency of the whole system, a mathematical model of the energy efficiency was required and could be obtained in many references, such as [18,20]. The electricity from the utility grid was assumed to be generated by power plants and supply to the MECS, and the exergy efficiency of the power plant was employed to calculate the exergy input of the coal, which was considered to be 33.5% [16]. Additionally, the energy efficiency of the power plant was assumed to be 32.4% in this paper, a val uewhich was applied to the energy efficiency calculation according to [20]. The dynamic electricity TOU tariff was applied to the scenario without the DR scheme in this example day and was divided into three parts: peak, mid-peak and off-peak, a division which was derived from [33]. Furthermore, the RTP was obtained according to the proposed model with the minimum $0.05 and maximum $0.25. The electricity pricing of the TOU and RTP schemes for each period are described in Figure 4.

Solar energy is utilized in the MECS extensively because it is environmentally friendly and free of charge. To acquire the energy output of the PV and STC systems, which are regarded as uncontrollable units, ambient temperature and solar radiation were retrieved from a building management system installed in Hong Kong and are illustrated in Figure 5 [11] in this typical day, and the installed area of the PV and the STC were set to be 500 and 300 m², respectively. The ambient temperature varied from 300 to 308 K, and the maximum solar radiation was about 900 W/m² at 12:00. Additionally, to evaluate the exergy in solar energy, the temperature of the sun was typically regarded as 6000 K [34].

Natural gas can be utilized in the CHP unit to generate electricity, while waste heat can be recovered to meet the cooling and heat demands of the AC and the HEU, respectively. It was assumed that the states of natural gas, such as temperature, pressure, and heat value were not changed during the operation of the MECS. The exergy factor of the natural gas was set to be 1.04, and its LHV was 50 MJ/kg, according to [14,35]. However, owing to the strong coupling between electricity and heat in the CHP, the operation of which lacks flexibility, the GB was considered to be the adjustable heat unit and TES was integrated to realize the heat charging and discharging during the operation.

To fulfill the cooling demand, the AC and the EC are employed in the MECS with complementary advantages. The AC could utilize low-quality heat to generate cooling energy, while its COP was low compared to EC. An EC with high COP consumes high-quality electricity. In addition, CES saves the received cooling energy and then releases the stored energy to jointly meet the cooling demand with the AC and the EC.

The efficiency, COP, and energy loss rate of relative devices are shown in Table 1. The restriction conditions are summarized in Table 2. The simulation interval Δt was defined as 1 hour, and the loop count i_max was set to be 20 in the ε-constraint method.

4.2. Results Analyses and Discussion

In Case 1, the economic–exergetic scheduling model was established with the TOU tariff, and DR was not considered. In Case 2, the MO scheduling model was integrated with the DR program and the RTP tariff, and both models were solved by the ε-constraint method to obtain the Pareto front, as shown in Figure 6. To determine a trade-off between the economic and exergetic objectives, the LINMAP approach was then applied, and the results are illustrated in Table 3. The details of the ε-constraint method and LINMAP approach were illsustrated in Figure 2. For Case 1, Solution #17 was the compromise scheme, while Solution #16 was the preferred point for Case 2, both of which are marked in red on each Pareto carve in Figure 6. Detailed analyses of these two selected solutions were conducted, the results of which are shown below.

The energy operation costs for Cases 1 and 2 were, respectively, $3515.84 and 3490.15 with corresponding exergy inputs of 74,389.45 and 73,644.09 kWh. The exergy efficiency for each case was 21.91% and 21.86%, and the energy efficiency for each case was 92.98% and 93.71%, respectively. It was found that the energy cost and exergy input both decreased with DR, which means that the more preferred operation strategy was arranged to simultaneously boost these two different goals after the implementation of the RTP-based DR, and the exergy or energy efficiencies had few variations between these two cases. However, it was noted that the energy efficiency for each case was more than 90% which was because the EC unit possessed a high COP. Therefore, the energy efficiency index only denotes the quantity performance, which neglects the quality evaluation of different energy carriers. If the energy efficiency index is adopted, the energy-saving performance might not be reasonably evaluated. Hence, exergy efficiency might be more rational for the assessment of an MECS compared to the energy efficiency index.

The purchased electricity from the grid is described in Figure 7a, and the optimal scheduling results for electricity are presented in Figure 8 for both cases at each time interval. In Case 1, the electricity was not purchased during the periods of 8:00–20:00 when the TOU tariffs were at peak or mid-peak, as shown in Figure 4. However, the electricity was purchased not only during the off-peak periods 1:00–7:00 and 23:00–24:00 but also during the time interval of 21:00–22:00, even if the price rate was at mid-peak, which was caused by the electricity needs of the EC and demand that could not solely be fulfilled by the CHP unit. Therefore, more electricity was required through the grid when the PV output was zero during the night. Compared to Case 1, the electricity was not bought from the grid when the RTP tariff was high in Case 2, especially during the period of 20:00–22:00. This contributes to the fact that the electricity needs and demand of the EC could be met by the CHP unit alone, as shown in Figure 8b. The electricity demand was reduced under the implementation of DR, which is illustrated in Figure 3. Hence, electricity consumption was more rationally and economically arranged. By comparing Figure 8a,b, it can be seen that there were many differences at 3:00–4:00 and 23:00 other than the 21:00–22:00 time periods. The electricity demands varied slightly, while the RTP rates were obviously less than TOU tariffs during 3:00–4:00, as illustrated in Figure 4. Therefore, the power energy purchased from the utility grid in Case 2 was more than in Case 1 during 3:00–4:00. Different results were obtained during 21:00–23:00, when the electricity demands were at peak and curtailed by the DR scheme with higher RTP tariffs in Case 2. Hence, less power energy was purchased from the grid, which was because of the fact that the new price mechanism more rationally modified the electricity consumption, and the operation strategy was optimized in this paper.

The purchased natural gas from the pipeline is presented in Figure 7b, and the distribution ratios for the CHP unit and the GB are illustrated in Figure 9 for each case at each hour. The purchased gas for both cases was distinct in some periods, owing to the DR program. The amount of gas needed in Case 1 was more than in Case 2 during 3:00–4:00, and more natural gas was purchased in Case 2 at 23:00, as described in Figure 7b. Hence, the scheduling results for natural gas and grid electricity were complementary at these periods because more natural gas are purchased and less electricity were needed when the RTP tariffs were high at some periods, and vice versa. From Figure 9, it can be seen the GB served as an auxiliary heat source to adjust the imbalance between heat/cooling demands and electricity demands due to the strong coupling of the CHP. For example, when the heat and cooling demands were at peak and the electricity demands were low during 14:00–15:00, the heat energy produced by the CHP could not meet the thermal requirements of the AC and the HE; therefore, the GB started to make up for the heat deficiency for both cases, as shown in Figure 9. However, it must be noted that the GB operated not only during 14:00–15:00 but also at 13:00 in Case 1 due to the lack of the filling valley effect of the DR scheme. In Case 2, the electricity demands were increased at 13:00 compared to Case 1; hence, the electricity generated by the CHP is higher than that in Case 1, and the heat energy produced by the CHP could fulfill the thermal needs of the AC and the HE. To this end, the GB unit did not have to start up at 13:00 in Case 2.

The optimal scheduling results of heat energy for demands and devices are presented in Figure 10 and Figure 11, and the cooling energy balance graphs are described in Figure 12. Though the DR scheme was applied to the electricity demands, the heat and cooling operation strategies were affected due to the energy coupling. In Case 2, TES released the thermal energy during 3:00–4:00, while it was still in Case 1, as illustrated in Figure 10. This was owed to the fact that the CHP unit generatds less electricity and thermal energy, while the utility grid supplied more electricity in Case 2, as shown in Figure 8 and Figure 11. Hence, the heat energy vacancy was compensated by TES. In addition, the heat energy demands were fulfilled not only by the HE but also by TES discharging in Case 1 at 23:00, while more heat energy was generated by the CHP unit in Case 2, so TES needed charge from the surplus thermal energy as shown in Figure 10 and Figure 11. This distinction of TES operation mode is obviously dependent on the RTP-DR program. Moreover, the cooling energy operation strategies have been modified in case 2 during some periods as well. Due to the low RTP electricity tariffs at 3:00 and 4:00 in case 2, more electricity energy are purchased from the grid and are converted to cooling energy by the EC unit and then are charged in the CES device so that the surplus cooling energy could be utilized at the peak of cooling demands. In addition, in Figure 12 it can be summarized that the EC module is the main cooling device to fulfill the cooling demand compared with the AC unit, which contributes to the fact that the EC device possesses higher COP by consuming high-quality electricity energy which means lower energy cost.

The state variations of TES and CES for the whole day are described in Figure 13. For TES in Figure 13a, it releases the heat energy at the peak during 12:00–14:00 and 21:00–22:00 for both cases to compensate for the heat deficits. In addition, TES discharged the energy during 3:00–4:00 in Case 2, which contributed to the fact that more electricity was purchased from the grid because of the low RTP tariffs, and the CHP unit generated less electricity and heat energy. Similarly, the CES unit released cooling energy during 12:00–14:00 and 20:00–21:00 when the cooling demands were at the peak, as illustrated in Figure 13b. With the coordination of TES and CES systems, the MECS could be economically and efficiently operated to meet diverse energy needs.

4.3. Sensitivity Experiments

In this section, the effect of the DR scheme is investigated with the consideration of different elasticity coefficients. Twenty-one configurations were developed with different combinations of self- and cross-elasticity coefficients, as presented in Table 4. In addition, to more intuitively investigate the effect of DR program, Equation (29) could be reshaped as:

_LeRTP(t)=_Le(t){1+e(t,t)Ps(t)+e(t,k)Pc(t)}

Ps(t)=_cRTP(t)−_cTOU(t)_cTOU(t)

Pc(t)=∑k=1k≠t24_cRTP(k)−_cTOU(k)_cTOU(k)

The values for Ps and Pc at each hour are illustrated in Figure 14, where it can be seen that the values of Pc were always negative while Ps was positive or negative at different time intervals. In Figure 15, the electricity demands with DR implementation for each configuration are illustrated. The electricity loads were shifted from peak to valley, and the demand curves were more smoothly modified with the increase of the absolute value of the self-elasticity coefficient, as shown in Figure 15a. In Figure 15b, the peak loads are distinctly curtailed while the valley loads are changed little, a difference caused by the fact that the value of Pc was away from zero at peak while Pc was close to zero at the valley, as shown in Figure 14.

After the MO operation strategy was applied to each configuration, the energy cost, exergy input, exergy efficiency and energy efficiency of the preferred solution for each case are illustrated in Figure 16. The economic and exergetic objectives were all reduced with the absolute value increase of the elasticity coefficients, as illustrated in Figure 16a,b, which means that a better Pareto front was obtained with the stronger implementation of the DR scheme. In addition, the changing trends of exergy efficiency and energy efficiency were completely opposite, as seen in Figure 16c,d. However, the variation amplitudes for both exergy and energy efficiency were very small, which means that the DR scheme had less influence on exergy and energy efficiency. Therefore, to distinctly promote the effect of RTP-based DR, a better identification of self- and cross-elasticity coefficients must be conducted since the pricing models of RTP play an important role.

5. Conclusions In this paper, an MECS was employed to fulfill the diverse energy demands with the coupling of different energy carriers by energy conversion devices and energy storage instruments. In addition, the economic–exergetic optimal scheduling model with an RTP-based DR scheme for an MECS was introduced to decrease operation costs and achieve energy savings by considering different qualities of energy. The ε-constraint approach was applied to this MO scheduling model to obtain the Pareto front, and a trade-off was conducted between these two objectives by the LINMAP approach. Moreover, sensitivity analyses of self- and cross-elasticity coefficients were explored and showed that the electricity load curve was modified to be smoother with the increase of the absolute value of the elasticity coefficient, and the economic and exergetic objectives were simultaneously decreased, while the exergy and energy efficiencies were not distinctly changed. Therefore, with the consideration of the RTP-based DR program in the optimal scheduling model, the MECS operator could more rationally allocate different energy carriers and simultaneously decrease energy costs and exergy inputs. However, to successfully implement the DR scheme, elasticity coefficients need to be more accurately identified, a need which will be investigated in the future.

Figure 1. An MECS (multiple energy carrier system) coupled with diverse energy carriers.

Figure 2. Flowchart of multi-objective (MO) operation strategy.

Figure 3. Energy demands of a typical summer day: (a) electricity demand and (b) heat and cooling demands.

Figure 4. TOU (time of use) and RTP (real-time pricing) tariffs.

Figure 5. Ambient temperature and solar radiation.

Figure 6. Pareto front with and without demand response (DR).

Figure 7. The purchased electricity and natural gas: (a) purchased electricity; (b) purchased gas.

Figure 8. Optimal scheduling results for electricity: (a) without DR; (b) with DR.

Figure 9. Natural gas distribution ratio: (a) without DR; (b) with DR.

Figure 10. Optimal scheduling results for heat: (a) without DR; (b) with DR.

Figure 11. Heat balance for each device: (a) without DR; (b) with DR.

Figure 12. Optimal scheduling results for cooling: (a) without DR; (b) with DR.

Figure 13. State variation of energy storage devices: (a) state of TES; (b) state of CES.

Figure 14. Values of Ps and Pc at each hour.

Figure 15. Load profile of DR for each configuration: (a) Configurations 1-11 and (b) Configurations 12-21.

Figure 16. Energy cost, exergy input, exergy efficiency and energy efficiency for each configuration. (a) Energy cost and exergy input for Configurations 1-11; (b) energy cost and exergy input for Configurations 12-21; (c) exergy efficiency and energy efficiency for Configurations 1-11; (d) exergy efficiency and energy efficiency for Configurations 12-21.

Figure 16. Energy cost, exergy input, exergy efficiency and energy efficiency for each configuration. (a) Energy cost and exergy input for Configurations 1-11; (b) energy cost and exergy input for Configurations 12-21; (c) exergy efficiency and energy efficiency for Configurations 1-11; (d) exergy efficiency and energy efficiency for Configurations 12-21.

Parameters	Value	Parameters	Value	Parameters	Value
ηPV	0.157	ηHEU	0.9	ηTESdch	0.98
ηSTC	0.8	COPEC	4.2	δCES	0.02
ηCHPe	0.3	COPAC	0.8	ηCESch	0.97
ηCHPh	0.45	δTES	0.02	ηCESdch	0.95
ηGB	0.88	ηTESch	0.98

Parameters	Value/kW	Parameters	Value/kW	Parameters	Value/kW
PCHP,mine	1000	PAC,minc	100	PCES,mindch	0
PCHP,maxe	200	PAC,maxc	1000	PCES,maxdch	300
RPCHP,mine	500	PTES,minch	0	ETES(0)	250
RPCHP,maxe	−500	PTES,maxch	200	ECES(0)	400
PGB,minh	30	PTES,mindch	0	ETESmin	50
PGB,maxh	300	PTES,maxdch	200	ETESmax	500
PEC,minc	150	PCES,minch	0	ECESmin	80
PEC,maxc	1500	PCES,maxch	300	ECESmax	800

	Case 1: Without DR					Case 2: With DR
i	f1($)	f2(kWh)	ρ1	ρ2	d	f1($)	f2(kWh)	ρ1	ρ2	d
1	3704.41	73574.05	1.000	0.000	1.000	3676.15	72847.62	1.000	0.000	1.000
2	3692.63	73618.71	0.947	0.014	0.947	3663.75	72882.55	0.947	0.014	0.947
3	3680.84	73667.91	0.895	0.029	0.895	3651.35	72923.42	0.895	0.029	0.895
4	3667.59	73727.98	0.836	0.048	0.837	3638.95	72959.54	0.842	0.044	0.843
5	3657.27	73766.30	0.789	0.060	0.792	3626.55	73000.49	0.789	0.059	0.792
6	3645.49	73816.30	0.737	0.076	0.741	3614.15	73042.87	0.737	0.076	0.741
7	3633.70	73864.69	0.684	0.091	0.690	3601.75	73086.82	0.684	0.093	0.690
8	3621.91	73913.88	0.632	0.107	0.641	3589.35	73132.19	0.632	0.111	0.641
9	3610.13	73963.08	0.579	0.122	0.592	3576.95	73180.72	0.579	0.129	0.593
10	3598.34	74012.27	0.526	0.137	0.544	3564.55	73229.62	0.526	0.148	0.547
11	3586.56	74061.46	0.474	0.153	0.498	3552.15	73281.93	0.474	0.169	0.503
12	3574.77	74110.66	0.421	0.168	0.453	3539.75	73337.80	0.421	0.191	0.462
13	3562.98	74159.85	0.368	0.184	0.412	3527.35	73401.26	0.368	0.215	0.427
14	3551.20	74209.47	0.316	0.199	0.373	3514.95	73464.81	0.316	0.240	0.397
15	3539.41	74263.74	0.263	0.216	0.341	3502.55	73546.73	0.263	0.272	0.378
16	3527.63	74318.07	0.211	0.233	0.314	3490.15	73644.09	0.211	0.310	0.374
17	3515.84	74389.45	0.158	0.256	0.301	3477.74	73793.13	0.158	0.367	0.400
18	3504.05	74731.77	0.105	0.363	0.378	3465.34	74162.12	0.105	0.511	0.522
19	3492.27	75732.11	0.053	0.677	0.679	3452.94	74642.98	0.053	0.698	0.700
20	3480.48	76761.51	0.000	1.000	1.000	3440.54	75420.63	0.000	1.000	1.000

#	e(t,t)	e(t,k)	#	e(t,t)	e(t,k)
1	0	0.01	12	−0.20	0
2	−0.04	0.01	13	−0.20	0.005
3	−0.08	0.01	14	−0.20	0.015
4	−0.12	0.01	15	−0.20	0.020
5	−0.16	0.01	16	−0.20	0.025
6	−0.20	0.01	17	−0.20	0.030
7	−0.24	0.01	18	−0.20	0.035
8	−0.28	0.01	19	−0.20	0.040
9	−0.32	0.01	20	−0.20	0.045
10	−0.36	0.01	21	−0.20	0.050
11	−0.40	0.01

Author Contributions

Y.H. formulated the research problem; S.L. reviewed the manuscript; Y.H., S.L., P.D. and Y.Z. modified the revised version; K.Y. proposed the mathematical model and algorithms; K.Y. and Y.Z. wrote the code and addressed the manuscript; W.Z. conducted the analysis of the conclusion; all of the authors were involved in writing this paper.

Funding

This work was funded by the Fundamental Research Funds for the Central Universities (Grant No. 2019MS099).

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Steven, C.; Arun, M. Opportunities and challenges for a sustainable energy future. Nature 2012, 488, 294-303.

2. Geidl, M.; Koeppel, G.; Favre-Perrod, P.; Klöckl, B.; Andersson, G.; Fröhlich, K. Energy hubs for the future. IEEE Power Energy Mag. 2007, 5, 24-30.

3. Krause, T.; Andersson, G.; Fröhlich, K.; Vaccaro, A. Multiple-energy carriers: Modeling of production, delivery, and consumption. Proc. IEEE 2011, 99, 15-27.

4. Geidl, M.; Andersson, G. Optimal Power Flow of Multiple Energy Carriers. IEEE Trans. Power Syst. 2007, 22, 145-155.

5. Wang, Y.; Zhang, N.; Kang, C.; Kirschen, D.S.; Yang, J.; Xia, Q. Standardized Matrix Modeling of Multiple Energy Systems. IEEE Trans. Smart Grid 2019, 10, 257-270.

6. Liu, T.; Zhang, D.; Wang, S.; Wu, T. Standardized modelling and economic optimization of multi-carrier energy systems considering energy storage and demand response. Energy Convers. Manag. 2019, 182, 126-142.

7. Moeini-Aghtaie, M.; Abbaspour, A.; Fotuhi-Firuzabad, M.; Hajipour, E. A decomposed solution to multiple-energy carriers optimal power flow. IEEE Trans. Power Syst. 2014, 29, 707-716.

8. Moeini-Aghtaie, M.; Dehghanian, P.; Fotuhi-Firuzabad, M.; Abbaspour, A. Multiagent genetic algorithm: An online probabilistic view on economic dispatch of energy hubs constrained by wind availability. IEEE Trans. Sustain. Energy 2014, 5, 699-708.

9. Yi, J.-H.; Ko, W.; Park, J.-K.; Park, H. Impact of carbon emission constraint on design of small scale multi-energy system. Energy 2018, 161, 792-808.

10. Khorsand, H.; Seifi, A.R. Probabilistic energy flow for multi-carrier energy systems. Renew. Sustain. Energy Rev. 2018, 94, 989-997.

11. Lu, Y.; Wang, S.; Sun, Y.; Yan, C. Optimal scheduling of buildings with energy generation and thermal energy storage under dynamic electricity pricing using mixed-integer nonlinear programming. Appl. Energy 2015, 147, 49-58.

12. Li, Z.; Xu, Y. Optimal coordinated energy dispatch of a multi-energy microgrid in grid-connected and islanded modes. Appl. Energy 2018, 210, 974-986.

13. Brahman, F.; Honarmand, M.; Jadid, S. Optimal electrical and thermal energy management of a residential energy hub, integrating demand response and energy storage system. Energy Build. 2015, 90, 65-75.

14. Di Somma, M.; Yan, B.; Bianco, N.; Graditi, G.; Luh, P.; Mongibello, L.; Naso, V. Operation optimization of a distributed energy system considering energy costs and exergy efficiency. Energy Convers. Manag. 2015, 103, 739-751.

15. Torio, H.; Angelotti, A.; Schmidt, D. Exergy analysis of renewable energy-based climatisation systems for buildings: A critical view. Energy Build. 2009, 41, 248-271.

16. Krause, T.; Kienzle, F.; Art, S.; Andersson, G. Maximizing exergy efficiency in multi-carrier energy systems. In Proceedings of the IEEE Power and Energy Society General Meeting, Providence, RI, USA, 25-29 July 2010; pp. 1-8.

17. Hepabsli, A. A Key Review on Exergetic Analysis and Assessment of Renewable Energy Resources for a Sustainable Future. Renew. Sustain. Energy Rev. 2008, 12, 593-661.

18. Wang, J.; Yang, Y. Energy, exergy and environmental analysis of a hybrid combined cooling heating and power system utilizing biomass and solar energy. Energy Convers. Manag. 2016, 124, 566-577.

19. Ahmadi, P.; Rosen, M.A.; Dincer, I. Multi-objective exergy-based optimization of a polygeneration energy system using an evolutionary algorithm. Energy 2012, 46, 21-31.

20. Elizondo, L.M.R.; Paap, G.; Ammerlaan, R.; Negenborn, R.R.; Toonssen, R. On the energy, exergy and cost optimisation of multi-energy-carrier power systems. Int. J. Exergy 2013, 13, 364.

21. Yan, B.; Di Somma, M.; Bianco, N.; Luh, P.B.; Graditi, G.; Mongibello, L.; Naso, V. Exergy-based operation optimization of a distributed energy system through the energy-supply chain. Appl. Therm. Eng. 2016, 101, 741-751.

22. Luis, A.A.; Edwin, R.; Francisco, S.; Victor, H. A review and analysis of trends related to demand response. Energies 2018, 11, 1617.

23. Aalami, H.; Moghaddam, M.P.; Yousefi, G. Demand response modeling considering Interruptible/Curtailable loads and capacity market programs. Appl. Energy 2010, 87, 243-250.

24. Albadi, M.; El-Saadany, E.; Albadi, M. A summary of demand response in electricity markets. Electr. Power Syst. Res. 2008, 78, 1989-1996.

25. Erdinc, O.; Paterakis, N.G.; Mendes, T.D.P.; Bakirtzis, A.G.; Catalão, J.P.S. Smart Household Operation Considering Bi-Directional EV and ESS Utilization by Real-Time Pricing-Based DR. IEEE Trans. Smart Grid 2015, 6, 1281-1291.

26. Zhi, C.; Lei, W.; Yong, F. Real-Time Price-Based Demand Response Management for Residential Appliances via Stochastic Optimization and Robust Optimization. IEEE Trans. Smart Grid 2012, 3, 1822-1831.

27. Aghamohamadi, M.; Hajiabadi, M.E.; Samadi, M. A novel approach to multi energy system operation in response to DR programs; an application to incentive-based and time-based schemes. Energy 2018, 156, 534-547.

28. Ma, T.; Wu, J.; Hao, L. Energy flow modeling and optimal operation analysis of the micro energy grid based on energy hub. Energy Convers. Manag. 2017, 133, 292-306.

29. Petela, R. Exergy of undiluted thermal radiation. Sol. Energy 2003, 74, 469-488.

30. Candau, Y. On the exergy of radiation. Sol. Energy 2003, 75, 241-247.

31. Aalami, H.; Yousefi, G.R.; Parsa Moghadam, M. A MADM-based support system for DR programs. In Proceedings of the 43rd Internatioal Universities Power Engineering Conference, Padova, Italy, 1-4 September 2008; pp. 1-7.

32. Aghaei, J.; Baharvandi, A.; Rabiee, A.; Akbari, M.A. Probabilistic PMU Placement in Electric Power Networks: An MILP-Based Multi-objective Model. IEEE Trans. Ind. Inform. 2015, 11, 332-341.

33. Ma, T.; Wu, J.; Hao, L.; Li, D. Energy flow matrix modeling and optimal operation analysis of multi energy systems based on graph theory. Appl. Therm. Eng. 2019, 146, 648-663.

34. Chow, T.; Pei, G.; Fong, K.; Lin, Z.; Chan, A.; Ji, J. Energy and exergy analysis of photovoltaic-thermal collector with and without glass cover. Appl. Energy 2009, 86, 310-316.

35. Ahmadi, P.; Dincer, I. Exergoenvironmental analysis and optimization of a cogeneration plant system using Multimodal Genetic Algorithm (MGA). Energy 2010, 35, 5161-5172.