To ensure the EPS’s ability to meet load demands, a detailed power flow analysis is essential. This analysis employs Equations 1, 2 to establish equality constraints, which are further refined in Equations 3, 4. Equation 5 extends the power flow balance equation to incorporate CHP, P2G, SPV, and wind turbines, with the corresponding operating constraints provided in Equations 6, 7. P G i − P D i = V i ∑ j = 1 n V i Y ij cos δ i − δ j − θ ij , (1) Q G i − Q D i = V i ∑ j = 1 n V i Y ij sin δ i − δ j − θ ij , (2) P Gi min ≤ P G i ≤ P Gi max , (3) V i min ≤ V i ≤ V i max , (4) ∑ k = 1 G P k ⁢ t + P CHP + P Wt + P PV + P ESS − D EB − D P 2 G − D L − D Loss = ∑ j = 1 n V n ⁢ V m ⁢ Y nm ⁢ cos ⁢ δ n − δ n − θ nm , (5) 0 ≤ P WT ≤ P Wt max , (6) 0 ≤ P PV ≤ P PV max . (7) Here, D 1 t is the electric load of the EPS (Dong et al., 2019; Leung, 2015; Yu et al., 2018).

This article investigates the impact of demand-side management (DSM) and hybrid renewable energy storage system (HESS) on the performance of the IES. In this context, a detailed model of these technologies is presented.

DHSs transport thermal energy between sources and consumers using networks of supply and return pipelines. Each node in the DHS has three distinct temperatures: supply, return, and ambient temperatures. The return temperature is always lower than the supply temperature due to energy consumption. Thermal and hydraulic modeling is essential to analyze the DHS. The hydraulic model ensures the continuity of mass flow, and the thermal model governs heat balance, temperature drop, and temperature mixing at the nodes. The modeling of the district heating system is based on Equations 24–30. Equation 26 quantifies the temperature reduction within a pipe, which is influenced by the heat transfer coefficient and the pipe diameter. It is important to acknowledge that heat loss along the pipe inevitably occurs because of the disparity between the water temperature and the ambient temperature. In this study, the total heat loss of the DHS is defined as the summation of individual heat losses across all pipes (Ghaffarpour, 2020; Klinkel, 2020; Liu, 2013). The operational limits are determined using Equations 29, 30. ∑ m ∈ Λ n m m n , t = 0 , ∀ m , n ∈ Λ D H S , ∀ t ∈ T , (24) H CHP + H H S i n − H H S O U T + H GB + H EB − H L − H L o s s = C P . m m n . τ m n i n − τ m n o u t , ∀ m , n ∈ Λ D H S , (25) τ m n i n − τ m n a = e λ m n . L m n C P . m m n . τ m n o u t − τ m n a , ∀ m , n ∈ Λ D H S , (26) τ m n o u t ∑ m ∈ Λ n m m n = ∑ m ∈ Λ n m m n . τ m n i n , ∀ m , n ∈ Λ D H S , (27) τ m n i n / o u t min ≤ τ m n i n o u t ≤ τ m n i n o u t max , (28) H CHP min ≤ H C H P ≤ H CHP max , (29) m mn min ≤ m m n ≤ m mn max . (30)

The NGS network comprises loads, pipelines, gas compressors (GCs), storage units, and P2G units. The following equations were employed for the NGS modeling process. Equation 31 represents the pipeline flow equation under a constant gas pressure condition. Additionally, the gas flow balance at each node is governed by Equation 32. Equation 33 quantifies the gas consumption of the GC, and the constraints governing the NGS variables are defined in Equations 34, 35. P n , t 2 − P m , t 2 = Z m n S m n 2 , ∀ m , n ∈ Λ N G S , ∀ t ∈ T , (31) ∑ g ∈ Ω G S Q g , t G S + Q G S i n − Q G S o u t Q G B ⁢ Q C H P − Q G C + Q P 2 G − Q L = ∑ m ∈ Λ N G S S m n , t , ∀ m , n ∈ Λ N G S , ∀ t ∈ T , (32) H k G C = B k S k . P o u t P i n Z b − 1 , (33) P n , t 2 min ≤ P n , t 2 ≤ P n , t 2 max , (34) S mn min ≤ S m n ≤ S mn max . (35) Here, S_mn is the gas flow in a pipe of length m-n; P_m and P_n are the pressures at nodes m and n, respectively; and Z_mn is a fixed number; in addition, Z_k and B_k are constants related to the compression factor and the working conditions of the compressor, respectively. S_k is the gas flow through the compressor k, and P_in and P_out are the compressor inlet and outlet pressures, respectively. For GCs, τ represents the gas flow consumed by compressor k and can be expressed as follows (Equation 36): τ k = α + β H k + γ H k 2 . (36) Here, α , β , and γ are compressor consumption coefficients (Wei and Wang, 2020; Li et al., 2003; Woldeyohannes and Abd Majid, 2021).

The MOO emerges as a method for solving an optimization problem with more than one objective. It is based on Equation 51: Minimize F x = min f 1 x , f 2 x , … , f n x   h i x = 0 , i = 1 , 2 , … , m , g i x ≤ 0 , i = 1 , 2 , … , p . (51) Here, x = [x₁, x₂, . . ., x_k]^T is the vector of decision variables, g_i is the ith unequal condition, h_i is the ith equality condition, m is the number of unequal conditions, and p is the number of equality conditions (Bhesdadiya et al., 2016; Dhiman et al., 2020; Zhang et al., 2021). The objective functions are formulated as follows:

This study formulates the optimization problem for minimizing losses and total cost in the optimal operation of the IES as a multi-objective problem. Different scenarios for loads and generating units are considered. The formulations for minimizing total losses and total costs are presented in Shaabani et al. (2017), Bhesdadiya et al. (2016), Yuan et al. (2020), Chen and Wang (2020), and Pansota et al. (2021).

The proposed EMS introduces a robust decision-making framework for optimizing the operation of an IES 24 h in advance. It effectively addresses the uncertainties inherent in RESs, electrical loads, and demand response. The proposed model encompasses two objective functions: minimizing the overall energy cost and emissions and minimizing the overall energy loss. Due to the model’s highly nonlinear nature, especially concerning thermal and gas system variables, an iterative approach utilizing the NSGA-III algorithm is employed to address the complex optimization problem. The optimal solution is obtained via solution reproduction during the iterative process.

In this paper, a clustering method called the Sim&Corrloss algorithm is used in scenario analysis to avoid unreliable decision-making. Unreliability results from the interaction between random variables and their correlations. The clustering algorithm used is based on the similarity and dependence of the saved scenarios after their reduction. The clustering method must be designed to preserve the correlation between the random variables. In the Sim&Corrloss clustering method, the objective function is developed based on the correlation loss and similarity functions as Equations 68–70: max ∑ m = 1 N l ∼ ∑ i ∈ I m S i m ⁡ ξ i − ξ m − β C o r r l o s s ⁡ ξ i − ξ m s . t ⋃ m = 1 N l = Ω l , ∑ i ∈ Ω l p i = 1 I m ∩ I m ′ = 0 , ∀ m ≠ m ′ , m , m ′ ∈ 1 , 2 , … , N l , (68) s i m ξ i , ξ j = 1 d ∑ k = 1 d 1 − p i p j p i + p j ξ i k − , ξ j k max ξ l k − min ξ l k + ε , (69) c o r l o s s ξ i , ξ j = ∑ i = 1 d − 1 ∑ j = i + 1 d Δ ρ i j 2 . (70)

In the Sim&Corrloss clustering algorithm, the parameters are chosen based on a balance between maintaining similarity between the original scenario set and minimizing correlation loss during the scenario reduction process. In this algorithm, the ratio β is crucial as it controls the trade-off between the correlation loss and the similarity during scenario reduction, and the objective is to find an optimal value of β that yields the best performance while ensuring stability in the reduced scenario set. In this paper, β is tested one by one from 0 to 1 with a 0.1 interval through a large number of simulation tests. This method is explained in more detail in Hu and Li (2019). To lessen the scenarios, the average values of solar irradiance, wind speed, and load demand over 24 h are calculated, and then the Sim&Corrloss clustering method is used.

The test IES used to validate the proposed optimal management strategy is shown in Figure 2. All the data used in the study can be found in Yu et al. (2018). The simulation is performed over a period of 1 day, and the results of the energy analysis are obtained in hourly intervals. The parameters of the algorithms are listed in Table 1. The values of the parameters required for the simulation are presented in Tables 2–10. As shown in Figure 2, the test case system contains various energy components, including an electric boiler, a gas boiler transformer, a CHP plant, an ESS, an HSS, a GS, a P2G converter, a wind turbine, an SPV unit, a hot water pipeline, a compressor, and a gas pipeline. Each subsystem in the IES has its own converters with unique energy dispatch coefficients and energy conversion efficiencies. The required parameters of the components are listed in Table 3. In summary, both electric and heat loads are included in the DRP. The maximum generation of the wind turbine is set to 50 MW, and the capacity of the SPV unit is set to 27 MW. Additional data required for these sources can be found in Yu et al. (2018). To demonstrate the validity of the proposed technique, 1,000 scenarios are generated with random variables, which are then reduced to 30 scenarios using a clustering algorithm.

The optimal scheduling of the IES in the most probable scenario, aiming to minimize the total losses and total operating and emission costs, is presented in Figures 3–5. Power generation units such as the SPV unit, wind turbines, and CHP plant are assumed to be capable of injecting both active (P) and reactive (Q) power, and their installation location is considered to be the PQ bus. The optimal Pareto front obtained in the scenario with the highest probability of occurrence is shown in Figure 6, which depicts the total losses and total operating costs of the system over a 24-h period. The statistical analysis of the objective functions defined in the optimization process is presented in Table 11. These results indicate that the changes in the objective functions remain within the specified range when the load demand, wind speed, and solar irradiance change. This suggests that the proposed method can effectively handle changes in input variables and their impact on output variables. Additionally, it emphasizes the importance of considering all aspects of the problem in short-term studies. Finally, the obtained statistical values are compared to the results obtained using MOSOA. The optimal Pareto front presented in Figure 6 shows the compromise point determined using NSGA-III for the optimization problem. The total operating and emission costs are 1.1602 E+05 $, and the total loss value is 0.1188 E+02 MW. The corresponding values for the objective functions obtained using MOSOA are 1.1593 E+05 $ and 0.1189 E+02 MW. When it comes to reducing the loss value for a true operating point, NSGA-III demonstrates superior performance, whereas MOSOA performs better in reducing total costs. As a result, a trade-off must be struck between minimizing loss and minimizing cost in operational optimization.

During the time intervals [01:00, 07:00] and [21:00, 24:00], the SPV plant does not generate any electricity, and the electricity demand is met by other power sources, whereas the ESS is charged. The CHP plant also supplies electricity during these periods. When the SPV plant generates more electricity, the net consumption of the electric boiler and P2G increases to meet the heat and gas demand. During the periods of low wind power generation, the load on the P2G decreases, but the output of thermal power plants increases. The operational results of the heat generation units and storage device are depicted in Figure 4. As shown, the dispatchable unit’s heat generation increases, whereas the gas boiler’s heat generation decreases. This necessitates an expansion of both gas and electricity generation. Heat storage ensures that the heat surplus or deficit throughout the day is compensated.

The power generated by the RES is proportional to wind speed and solar irradiance, and this is taken into account in the optimization process. However, the charging and discharging of storage devices are generally proportional to energy consumption, considering the type of load in each subsystem. This means that storage devices should be charged when energy consumption is low and discharged when consumption is high. Optimal management is implemented with a DRP for electrical and thermal loads. Load shifting is employed to smooth the load profile.

These figures demonstrate that the integration of renewable energy sources and energy storage enhances system sustainability, resilience, and flexibility to address critical conditions.

As shown in Figure 7, by implementing DRPs and shifting loads to off-peak hours, storage devices are expected to discharge during peak demand periods and charge during off-peak hours. This optimizes storage utilization, reduces load fluctuations, and enhances system robustness.

A statistical analysis of the objective functions defined in the optimization process is presented in Table 11. It contains quantitative statistical analysis of the objective functions, total loss and total cost, for the three algorithms: NSGA-III, MOSOA, and GAMS. The headers represent statistical metrics: mean (average value), Std (standard deviation), max (maximum value), and min (minimum value). The associated data entries are the computed values for total loss (in per-unit, Pu, with base 1e5 MW) and total cost (in per-unit, Pu, with base 1e5 $). The minimum and maximum values of the objective functions for total loss and total cost are found to be within the ranges [0.1178–0.1581] and [1.1521–1.5095], respectively. Statistical analysis of the optimization process reveals that the objective functions, total loss and total cost, are constrained within specified ranges. It is important to note that in statistical analysis, the average value represents the central tendency of the data, w the standard deviation quantifies the dispersion of the data around the mean value. Correspondingly, the average values for the loss function and the total cost function are found to be 0.1389 and 1.3319, respectively. The standard deviations for the loss function and the total cost function are 0.0116 and 0.1034, respectively. These results suggest that the changes in the objective functions, namely, total loss and total cost, will be constrained within the specified ranges when the load demand, wind speed, and solar irradiance undergo variations. This indicates the effective capture of input–output relationships by the proposed method. Moreover, the inherent complexity of energy problems underscores the significance of statistical analysis in short-term studies for optimal energy management systems. The obtained statistical values are compared with those obtained via the MOSOA algorithm and GAMS solver.

A comparative analysis was conducted among NSGA-III, MOSOA, and the GAMS solver. The findings clearly indicate that NSGA-III is superior to MOSOA in minimizing the cost function. Conversely, as the cost function decreases, the loss function tends to increase, suggesting a trade-off between these two objectives. However, it is important to note that the performance of both algorithms may vary depending on the specific characteristics of the problem. Table 12 shows this conclusion. As shown in Table 13, in the scenario with the highest probability of occurrence, PAR is decreased by 5.4247% and 0.1691% in the electrical and gas networks, respectively. In the unscheduled operation scenario, where there is no optimization or scheduling, the PAR in both the electric and gas networks is the highest. This indicates greater pressure on the network during peak times, signifying inefficient energy usage and higher operational stress. In comparison, it can be said that with NSGA-III, the PAR decreases in both the electric and gas networks.

In addition, in the energy hub, based on the definition of an index to measure the energy performance of an energy hub, the total energy input is compared to the total energy output. This index is defined as the ratio between the total energy output and the total energy input. It is expressed as an energy efficiency index (EEI), whose higher value indicates higher efficiency and lower total energy loss. This index is evaluated under three operational modes. As shown in Table 14, in the scenario with the highest probability of occurrence, the results indicate that scheduling with both NSGA-III and GAMS algorithms improved the energy efficiency, with NSGA-III showing a more significant improvement.

• Most conventional energy management systems (EMSs) do not fully consider the state variables of the three subsystems in the IES. In contrast, the proposed approach comprehensively captures all the state variables based on load flow analysis, heat flow analysis, and gas flow analysis in the EPS, DHS, and NGS, respectively. This enables the optimization of loss function and cost function.

• The model does not consider the connection to upstream electricity, gas, and heat networks. This simplification was made to focus on the local energy system optimization. However, future research could explore the impact of interactions with larger networks and the potential for energy exchange.

It is important to note that these limitations do not diminish the value of this study but rather highlight potential avenues for future research and the development of more comprehensive models.