Multi-objective operation of distributed generations and thermal blocks in microgrids based on energy management system

The optimal multi-objective operation of microgrids with distributed generations and thermal block based on the energy management system is presented. In the thermal block, the combined heat and power system and boiler and thermal storage system supply the load of the block thermal. Therefore, the proposed strategy minimises three objectives of an microgrid’s operation, that is, cost, energy loss and voltage deviation functions. Also, it is subject to AC power ﬂow equations, system operation limits and distributed generations and thermal block constraints. The problem is modelled as a multi-objective optimisation method; then is converted to a single-objective problem model using the ε -constraint-based Pareto optimisation. Moreover, the proposed strategy includes uncertainties of the load, energy price and active power of renewable distributed generations; hence, a stochastic programming based on the hybrid Mont Carlo simulation and point estimate method is used to model these uncertain parameters. Then, the teaching–learning-based optimisation and ﬁreﬂy algorithm are utilised to solve the problem to achieve a reliable optimal solution. Finally, the proposed strategy is simulated on the IEEE 69-bus microgrid, and thus, the capabilities of the proposed strategy are extracted. , Maximum capacity of the distribution station (p.u.); S Lmax , Maximum capacity of the distribution line (p.u.); t , ST, Indices and set of simulation time; V, θ , Voltage magnitude (p.u.), voltage angle (rad); V min , V max , Minimum and maximum allowed voltage magnitudes (p.u.); V ref , Voltage magnitude of the slack bus (p.u.); w, S , Indices and set of scenario sample


INTRODUCTION
Due to the increased environmental pollution arising from mismanagement of fossil fuel consumption, the number of distributed generations (DGs) is also increasing in the distribution network [1]. Since the actual distribution system cannot control a large number of DGs [2], the network is divided into different microgrids (MGs) based on the smart grid concept [3]. Therefore, the optimal operation method can be implemented on MGs with different DGs if the network utilised a smart telecommunication platform [4]. Moreover, among DGs, the combined heat and power (CHP) system provides high performance and efficiency in power systems thanks to the simultaneous generation of electrical and thermal energies [5]. Hence, it is suitable to employ the model of thermal block including CHP, boiler and thermal storage system (TSS) in the MG operation problem, where boiler and TSS can be used as flexible sources for the CHP.
There is widespread research on the MG's operation in the presence of CHP. The authors in ref. [6] present the mathematical model of the CHP operation in an MG to investigate IET Gener. Transm. Distrib. 2021;15:1451-1462.
wileyonlinelibrary.com/iet-gtd performance of the DG against a set of network frequencies. In ref. [7], the optimal dispatch model of a combined cooling, heating and power (CCHP) system is presented for electrical and gas networks. Also, this model as a stochastic multi-objective method for CCHP in MG is proposed in ref. [8], where the model of the MG's operation cost and pollutant emission in a two-objective optimisation problem is considered and thus it uses the hybrid non-dominated sorting genetic algorithm with co-evolution theory and the beetle antennae search algorithm to solve this formulation. The dispatch model of the CHP in electrical networks is formulated in ref. [9], where the modified group search optimiser algorithm is used to solve the problem. Also, this case is expressed in ref. [10] based on subsidy strategies. The inflexible operation of CHPs in conventional combined electrical and heating networks leads to large wind power curtailments during the winter. The thermal inertia of aggregated buildings can provide additional operational flexibility to promote wind power accommodation. Hence, this statement is considered in ref. [11], where a day-ahead scheduling method is suggested to integrate these networks considering the buildings thermal inertia. In ref. [12], the economic power flow model considering the operation cost of the network and sources is formulated. In the same study, the CHP is used as a flexible source to improve the electrical network flexibility in the presence of wind system uncertainty, and a TSS is utilised in the CHP to improve the flexibility of the source in its heating part. The authors in ref. [13] solve the optimal power and heat flow problem in a CHP MG considering a robust model for selling/purchasing power price uncertainty in the energy market. Moreover, the dispatch model of the CHP with a TSS is formulated in ref. [14], where TSS can be used as a flexible source for the CHP. In ref. [15], a model of CHP/CCHP is used in energy hub that proposes an optimal operation model for electrical, gas and heating networks.
In most of the literature on the CHP management in electrical networks, the electrical model of the CHP is considered. Hence, it is used as a flexible source in the electrical network in the presence of renewable sources [6][7][8][9][10][11][12][13][14][15]. Nonetheless, the heating power of the CHP depends on active power of the CHP, so it is an inflexible source in the thermal network. Thus, it is suitable to use the real model of the CHP in electrical, heating, and gas energy parts in different problems, even in the optimisation problem of electrical power systems. Hence, the flexible energy management strategy can be implemented in this network to achieve a flexible scheme for thermal and gas blocks [12,14]. It is worth noting that the flexibility is defined as "the modification of generation injection and/or consumption patterns in reaction to an external price or activation signal to provide service within an electrical system" [16]. Hence, the flexible energy management refers to an energy management system (EMS) that can manage active power of flexible sources to improve system flexibility and network operation situations in a network with non-dispatchable sources. Moreover, the MG operation model considering economical and technical objectives in the presence of CHPs and TSSs has been less studied. Authors in ref. [8] consider a two-objective function including operation cost and pollutant emission models. Improving a network index may deteriorate another index of the network. For example, it is expected that to improve network economic situation, the local sources and storages in the grid inject high power into this system. Hence, it is expected that the overvoltage appears in the system and power loss increase significantly. Therefore, the multi-objective optimisation model needs to the network operation problem to achieve suitable situations for different indices of the power system.
To cope with these issues and explain a more accurate operation model of the MG, a multi-objective operation of an MG with thermal blocks and DGs based on an EMS, as shown in Figure 1, is proposed in this study. The problem is to minimise the energy loss and voltage deviation functions (ELF and VDF), and operation cost of the MG, including electricity, gas and heating energy costs as multi-objective economical and technical functions. Also, it is subjected to AC optimal power flow (OPF) equations in the presence of DGs and thermal blocks and system operation limits. Moreover, the proposed multi-objective problem is converted into a single-objective optimisation model using the ε-constraintbased Pareto optimisation method to obtain the best compromise solution. Also, the stochastic programing based on the combination of Mont Carlo simulation (MCS) and point estimate method (PEM) is used to model uncertainties of load, energy price, and renewable DGs power. This strategy includes a non-linear model, so the hybrid teaching-learning-based optimisation (TLBO) and firefly algorithm (FA) method is used to solve the proposed model. Briefly, the main contributions of this paper compared to previous works are summarised as follows: -Proposing a flexible scheme to the MG in the presence of renewable sources using CHP and micro-turbine, where the framework is implemented on thermal blocks using a boiler and a TSS; -Extracting a multi-objective problem to the energy management of the MG based on the best compromise solution between economic and technical indices; -Using a stochastic programming based on the combination of MCS and PEM to model the uncertainties of load, energy price and renewable DGs generation power.
The rest of the paper is organised as follows: Section 2 describes the stochastic multi-objective operation model for an MG with DGs and thermal blocks. Section 3 presents the solution method according to the ε-constraint-based Pareto optimisation method. Sections 4 and 5 address numerical simulations and the main conclusions of the paper, respectively.

PROPOSED PROBLEM FORMULATION
In this section, the optimal operation model of the MG in the presence of thermal block including CHP, boiler and TSS FIGURE 1 Proposed optimal operation scheme of the microgrid based on an energy management system is presented. The problem contains three objective functions: The MG operation cost function, ELF, and VDF. Also, it is subject to OPF equations in the presence of DGs, and thermal block constraints that consider the interaction between the CHP, boiler and TSS models, and the power balance equation between sources and the thermal load.

Objective function
The proposed multi-objective functions are modelled as Equations (1)-(3), which respectively refer to minimisation of MG operation cost, ELF and VDF. In Equation (1), the expected operation cost includes four parts, where the electrical network operation cost is expressed in the first part. The fuel cost of the non-renewable energy sources such as MT and CHP and boiler is formulated in the second part; the third part presents the thermal block sources, that is, CHP, boiler and TSS, the revenue from selling heating power to supply the thermal load. Also, the last part refers to DGs revenue that originates from selling and injecting active power into the MG. Note that the operation cost of the TSS is formulated in the third part of Equation (1), so that thermal block revenue will be reduced/increased if the TSS operates in the charging/discharging mode. Equation (2) describes the expected ELF that is equal to the sum of active power loss at the proposed operation horizon, that is, 24 h. Also, active power loss at hour t is related to the difference between active power generation by DGs and upstream network and load active power in that hour. Finally, the minimisation of the expected VDF is modelled in Equation (3), where it is expected that the voltage magnitude of all busses is close to the voltage of the slack bus during all simulation hours [16]. (1) . (2)

Microgrid-optimal power flow constraints in the presence of distributed generations
This section contains the AC power flow equations, DGs constraints, and system operation limits. The model of AC power flow is presented in Equations (4)-(8) that are nodal active and reactive power balance, Equations (4) and (5), active and reactive power flowed on the distribution line, Equations (6) and (7) and voltage angle of the slack bus, Equation (8) [17]. System operation limits including capacity limitations of the distribution line and the station and voltage limitation of the buses are expressed in Equations (9)-(11), [17]. It is noted that P G and Q G refer to active and reactive power of the distribution station. Furthermore, it is assumed that the station is connected to the slack bus; hence, the values of these variables are zero in other buses. Finally, the operation limits of the MT, WT and PV are formulated in Equations (12)- (14), respectively [18]. It is supposed that DGs provide only active power based on the IEEE 1547 standard, so the capability of reactive power control of these DGs is not modelled in these constraints.

Thermal block constraints
Note that this block includes the CHP, boiler, TSS and thermal load. The CHP injects active power into the MG and boiler provides the thermal load of the block. The TSS is used in this block to increase the flexibility of the block. Finally, the mathematical model of the proposed thermal block is expressed in Equations (15)- (26). Constraint Equation (15) refers to the heating power balance between different thermal sources and the load. The CHP's operation model is formulated in Equations (16)- (19), respectively, representing the heat power generation by the CHP, active and heat power limits of the CHP, and consumption gas power of the CHP [5]. Moreover, constraints Equations (20) and (21) refer to heat power generation limit of the boiler and consumption gas power of the boiler, respectively. Finally, the TSS model is presented in Equations (22)- (26), so that its charging and discharging rate limits are given in Equations (22) and (23), respectively. The stored energy of the TSS at different hours is based on Equations (24) and (25), and its limitation is provided in Equation (26) [6].
H CHPmin H BO min

Stochastic programming
In the proposed problem, Equations (1)- (26), the parameters of active, reactive and thermal loads, that is, P D , Q D and H D , the maximum generation power of the PV and WT, that is, P PVmax and P WTmax , price of electricity, gas and thermal energies, that is, ρ ele , ρ gas and ρ heat , are uncertain. It is noted that the value of these parameters obtain generally from forecasting method. Since this method includes forecasting error, thus, it is expected that these parameters are not included the deterministic situation, hence, they are known as uncertain parameter that is modelled by stochastic programming in this paper. Therefore, the stochastic programing based on the hybrid MCS and PEM is used in this section to model these uncertain parameters. In this approach, the MCS generates many scenario samples to load and energy prices, P PVmax and P WTmax , according to normal, beta and Weibull probability distribution functions, respectively [19]. In the next step, the mean and standard deviation values of different uncertain parameters are calculated. Finally, the PEM based on the (2n + 1) technique extracts the final scenario samples [20]. In this method, n is the total number of uncertain parameters and (2n + 1) refers to the total number of final scenario samples extracted by the PEM. Finally, the model of PEM (2n + 1) is written as follows [20]: Step 1: Determine the total number of uncertain parameters (n) Step 2: E(S h ) = 0, h = 1, 2, where h is the indice of the output torque, and E(f) refers to the expected value of f.
Step 3: Select an uncertain variable Z l .
Step 4: Calculate the skewness, Z l ,3 , kurtosis, Z l ,4 , of Z l using Equations (27) and (28): Terms z l and z l denote the mean and standard deviation values of the variable Z l . Also, E[(z l − z l ) 3 ] and E[(z l − z l ) 4 ]are calculated as follows: where, N is the total number of scenario samples generated by the MCS, and prob shows the probability of each scenario sample.
Step 5: Calculate two standard locations (ξ): Step 6: Obtain the two approximated locations: Step 7: Solve the deterministic problem model considering two approximated locations: Step 8: Calculate the impact weight (ω l ) of the points: Step 9: Update E(S h ): Step 10: Repeat Steps 3-9, while all uncertain parameters are considered in the calculations.
Step 11: Obtain the impact weight of the mean point (ω μ ): Step 12: Update E(S h ): Step 13: Calculate the mean (μ S ) and standard deviation (σ S ) of the output uncertain variable:

SOLUTION METHOD
The Pareto optimisation solves the proposed multi-objective problem in this section. The proposed approach helps to find and compare the set of acceptable solutions and present them to the decision maker who will choose the final solution among them. Moreover, the ε-constraint-based Pareto optimisation is used in this paper to model the proposed problem as a singleobjective function formulation [21]. According to this method, the function Cost in Equation (1) can be used as the objective function in the new problem and functions of ELF and VDF in Equations (2) and (3) are respectively constrained to ε elf and ε vdf as follows: (40) Subject to: Constraints (4) − (26) .
In the above problem, ε elf is between minimum and maximum values of the ELF, that is, ELF min and ELF max , and ε vdf varies between minimum and maximum values of the VDF, that is, VDF min and VDF max . Therefore, different values are found for functions Cost, ELF and VDF based on different values of ε elf and ε vdf , where the set of these points is called the Pareto front [21]. Hence, the fuzzy decision maker method [21] is used in this paper to achieve the best compromise solution. Finally, the details of this approach are as follows, where it follows the flowchart given in Figure 2 [21]: Step 1: Calculate linear fuzzy membership function (f ) values for each objective function based on different values of ε elf and ε vdf as follows: Step 2: Obtain the value ofa r = min(f r Cost ,f r ELF ,f r VDF )., where r = 1, 2, …, and it refers to different steps that are included various value for ε elf and ε vdf .
Step 3: The best compromise solution is related to a point that contains max (a 1 , a 2 , …).

Case study
In this section, the proposed method is implemented on the IEEE 69-bus MG [22]. This network includes a base power of 1 MVA and a base voltage of 12.66 kV, where minimum and maximum allowed voltages are set to 0.9 and 1.05 p.u., respectively [22]. The data of the peak electrical load and characteristics of the distribution lines and the station can be extracted from [22]. Additionally, electrical load values at different hours are equal to the product of the peak electrical load and hourly electrical load factor curve as plotted in Figure 3(a) [18]. In addition, this network includes three thermal blocks, the locations of which are tabulated in Table 1. The block contains a CHP, a boiler and a TSS, so that maximum/minimum active and heat power of the CHP is 200 and 120 kW, respectively. Also, efficiencies

FIGURE 3
Hourly curve of a) load factor, b) maximum active power generation of the photovoltaic and wind turbine and c) prices of various energies [18] t , h , l are set to 0.4, 0.4 and 0.08, respectively. Moreover, the maximum/minimum heat power of the boiler is 100 kW with efficiency of 0.85. The TSS with efficiency of 0.85 includes charging/discharging rate of 30 kW, and minimum and maximum allowed sorted energy of 10 and 100 kWh. The block consists of a thermal load with a peak value of 150 kW and its hourly thermal load factor curve is shown in Figure 3(a) [18].
There are MT, PV and WT in the proposed network, where their locations are expressed in Table 1. Minimum/maximum active power of the MT with efficiency of 0.4 is 100 kW, and the maximum active power generation of the PV and WT is plotted in Figure 3(b) [18]. Finally, the hourly electricity, gas, and heating energy price curves are illustrated in Figure 3(c) [18].

Evaluation of the best compromise solution
Results of the Pareto optimisation for the proposed problem, Equations (40)-(43), are presented in this section. In Table 2,  Figure 4 shows the Pareto front for the problem of Equations (40)-(43) according to results of Table 2. Accordingly, the MG's operation cost is reduced when VDF and ELF values are increased, because DGs are expected to inject high amounts of active power into the MG to minimise the operation  cost. Nevertheless, injecting high active power by local sources increases the overvoltage and active power loss. Therefore, voltage deviation and energy loss will be increased if the MG's operation cost is decreased. Finally, the best compromise solution to the proposed method found by FA, TLBO, CSA, KHO and GWO is reported in Table 3. According to Table 3, the optimal compromise point is ε elf = 2.687 MWh and ε vdf = 0.962 p.u. in the FA solver. In this point, the values of the MG's operation cost, ELF and VDF are $1882.4, 2.687 MWh and 0.962 p.u., respectively. Therefore, the MG's operation cost is increased approximately by $55.1 ($1680.3-$1625.2) or 35% ((1680.3-1625.2)/(1781.5-1625.2)) with respect to its minimum value in this condition. This value is 45% and 40% for ELF and VDF, respectively. It is seen that the optimal compromise point in the TLBO, CSA, KHO and GWO solvers is higher than that of the FA solver. This statement is achieved in Figure 5, which plots the convergence curves of the FA, TLBO, CSA, KHO and GWO solvers. Also, FA solver in comparison with these solvers can be obtained the optimal solution in the low convergence iteration (CI) and calculation time (CT) that are reported in Figure 5 for the mentioned solvers. Therefore, it can be said that the FA is more suitable in comparison to the other mentioned solvers, because it obtains a better optimal point than these solvers. To investigate the capability of the proposed solution method and algorithm, Table 3 presents the optimal value of the proposed objectives according to the proposed algorithm in ref. [9], multi-objective solution method given in ref. [8], multi-objective solution based on the normalised objective function solved by the FA [28], and multiobjective solution based on the weighted sum solved by the FA [29]. Accordingly, it is seen that FA can achieve the best optimal solution in comparison with solvers presented in refs. [8,9], and the proposed multi-objective solution in Section 3 of the paper can find the best solution compared with multi-objective methods in refs. [8,28,29].

4.2.2
Analysis of distributed generations' operation in the microgrid Figure 6 shows the daily active power curve of different DGs such as, PV, WT, MT and CHP related to the best compromise  and 6, the PV and the WT inject active power into the MG about their maximum generation capacity at all simulation hours. As these DGs include zero operation cost, the best option to minimise the MG's operation cost is to use renewable DGs in the MG. Also, based on Figure 6, MTs inject their maximum generation active power capacity into the MG at 18:00-24:00. However, in other hours, the DGs transfer low active power toward the MG so that they are switched off during 5:00-9:00. It should be noted that this action of MTs is based on minimising the MG's operation cost and the difference between electricity and gas energy prices. For example, the electricity energy price is less than gas energy price during 5:00-9:00; hence, MTs are switched off. Because this DG is reluctant to obtain high gas operation cost it cannot reduce the electricity operation cost in this condition. Moreover, electricity price is higher than gas price during 18:00-24:00; hence, MTs inject their maximum generation capacity into the network. In addition, the changes trend of the CHP's active power is the same as MTs active power, so that they inject high power at the period 18:00-24:00 than the period 1:00-17:00. The reason is that the difference between electricity and gas energy prices is low at hours 1:00-17:00 in comparison to period 18:00-24:00. Moreover, based on Figure 6, CHPs inject higher power compared to MTs, and this is due to the following facts: In this network, based on Section 4.1, the capacity of CHPs is greater than that of MTs, and CHPs generate active and heat power, while MTs produce only active power. Also, heat energy price is higher than gas energy price at all simulation hours according to Figure 3(c). Therefore, CHPs tend to generate heat power during the operation horizon, where it depends on active power based on equation (16). Hence, by generating active and heat power, the CHP can dedicate itself to different conditions of various energy prices as far as it is switched on. This is one of significant benefits of utilising a CHP rather than an MT.

4.2.3
Analysis of the thermal block operation Daily heat power curves of CHPs, boilers and TSSs in thermal blocks are plotted in Figure 7. Accordingly, CHPs generate high heat power in comparison with boilers and TSSs at more simulation hours. Because heat power of the CHP depends on the CHP's active power, the changes trend of heat and active powers of the CHP are related to the difference between electricity and gas energy prices and the difference between heat and gas energy prices to obtain the MG's minimum operation cost. Therefore, CHPs generate high heat power at periods 1:00-4:00, 8:00-11:00 and 13:00-24:00 compared to boilers, where this statement is due to high heat or electrical energy prices compared to the gas price at these periods. In other hours, boilers include high heat power in comparison with CHPs, because heat/electricity price is greater/less than gas price at these hours. Hence, it is said that a CHP with generation active and heating power is more economical than a boiler. In  addition, the TSS is charged at low thermal load hours, that is, during 1:00-4:00, that is overlapped with low heat energy price hours, where it receives heat power from the CHP and boiler in this condition. However, it injects this power into thermal block to supply thermal load of the block at 8:00-11:00, which is overlapped with peak thermal load and price hours. This action of the TSS is based on minimising gas and heat operation costs. Finally, the daily gas power curves of MTs, CHPs and boilers are plotted in Figure 8. Based on Figures 6 and 8, the changes trend of gas power of CHPs and MTs (P MT /η mt ) is the same as their active power, where this statement is confirmed in Equation (19). Also, the changes trend of the boiler's gas power is the same as its heat power according to Figures 7 and 8, where it is confirmed in Equation (21). In addition, in refs. [12,14], the TSS is used only as a flexible source to improve flexibility of thermal block. Nonetheless, note that this case is suitable if the CHP's capacity is small, because the TSS's capacity is generally low. Hence, boiler and TSS can obtain a suitable scheme to increase the flexibility of the thermal block when the capacity of CHP is large. This statement is based on results shown in Figure 7.  Evaluation of the capabilities of the proposed multi-objective energy management system in the microgrid In this section, the daily curve of the distribution station's apparent power and network's active power loss are shown in Figures 9 and 10, and bus voltages are plotted in Figure 11. Also, two cases (Case I and Case II) are investigated, where Case I/II refers to power flow analysis (MG without DGs)/the proposed problem Equations (40)-(43). According to Figure 9, the daily apparent power curve of the distribution station in Case II is shifted downward with respect to Case I due to active power injection by DGs in Case II. Moreover, the daily active power loss curve of the network, similar to the distribution station's apparent power curve, is shifted downward in Case II compared to Case I, according to Figure 10. This is because of the optimal management of active power injection by DGs. Finally, it is observed that the flat profile voltage in peak load hour, that is, at 20:00, is created in Case II in comparison with Case I. This is resulted from the optimal energy management of DGs based on Figure 11(a). The daily voltage curve of bus 65, which has a lower voltage than other buses according to Figure 11(a), is plotted in Figure 11(b). Accordingly, voltage of bus 65 is close to 1 p.u. in Case II and, based on the proposed multi-objective EMS framework, this statement refers to lower voltage drop in the MG for Case II. In addition, the MG's operation cost in  the problem describe by Equations (40)-(43), as deterministic and stochastic models according to results of Table 3 and FA solver, is expressed in Table 4. Based on this table, considering uncertainty in the proposed problem increases the operation cost compared to that of the deterministic model, so it adds a $96.4 cost to the operation cost in the deterministic method.

CONCLUSION
In this paper, the multi-objective operation of an MG based on an EMS is presented for an MG with DGs and thermal blocks. The thermal block contains a CHP, a boiler and a TSS to supply the thermal load. This problem aims to minimise the electricity, gas and heat operation cost, energy loss and voltage deviation functions as three objection function models. The formulated problem is subject to AC power flow equations, DGs and thermal blocks constraints and system operation limits. In the next step, the problem is converted into a single-objective method using the ε-constraint-based Pareto optimisation. Then, the stochastic programming according to the hybrid MCS and PEM method is employed to model the uncertainties of load, energy price and generation power of renewable DGs. Finally, two FA and TLBO solvers solve the problem; thus, based on obtained simulation results, it is seen that the FA solver includes lower operation cost, ELF and VDF than those of the TLBO algorithm. Additionally, according to the numerical results, the CHP with electricity and thermal energies generation is the more economical than the MT/boiler that generates only active/heating power. The TSS is a suitable option to reduce the MG's operation cost because it can be charged/discharged during low/peak thermal energy prices. Finally, the problem can improve/reduce the operation indices/energy cost using the proposed EMS.