TY - JOUR
AU - Zhang, Fengyang
AB - 1. Introduction In response to escalating global environmental pollution and fossil energy crises, nations worldwide are increasingly prioritizing low-carbon development and establishing emission reduction targets [1]. Achieving the dual objectives of carbon peaking and carbon neutrality necessitates structural transformation and upgrading of energy systems, where decarbonizing the power sector plays a pivotal role in accelerating national carbon reduction goals [2,3]. While this transition has imposed unprecedented pressures on China’s energy structure and economic development, the synergistic evolution of integrated energy systems (IES) and renewable energy sources not only addresses these challenges but also provides innovative pathways to enhance renewable energy penetration within a low-carbon framework [4]. Integrated Energy Systems transcend the operational silos of conventional energy systems, enabling economically sustainable development while ensuring efficient and eco-friendly energy utilization [5].The multi-energy synergy inherent in IES significantly enhances energy efficiency, reduces operational costs, and mitigates carbon emissions [6].Recent advancements demonstrate diverse methodological innovations:Tan et al. developed a dispatch model for park-level IES (PIES) integrating photovoltaic/thermal (PV/T) hydrogen production with an enhanced stepped carbon trading mechanism. Their refined carbon pricing scheme effectively promoted system-wide energy conservation and emission reduction [7]. Cheng et al. emphasized that leveraging complementary interactions among heterogeneous energy carriers in multi-energy systems (MES) enables substantial carbon emission abatement [8].Luo et al. established a carbon-green certificate co-trading market framework under combined carbon emission trading (CET) and green certificate trading (GCT) mechanisms. This approach eliminates traditional market barriers while improving renewable energy accommodation rates and reducing system emissions [9].Liu et al. proposed a day-ahead dispatch model for regional IES (RIES) incorporating information gap decision theory (IGDT).By analyzing renewable output uncertainties through heat demand response (HDR)-integrated geothermal system modeling, their framework demonstrated superior cost-effectiveness and emission reduction capabilities [10].Wang et al. introduced a stepped carbon trading mechanism with incentive-penalty structures and price-signal-driven demand response strategies, effectively balancing stakeholder interests while minimizing emissions in carbon-constrained environments [11].Zhao et al. devised a two-stage distributionally robust optimization model for wind-PV-storage hybrid systems. Incorporating stepped carbon trading, their methodology reduced reliance on conventional generation units, achieving significant reductions in both emissions and economic costs while validating the model’s adaptability and superiority [12]. In recent years, heuristic algorithms have gained widespread attention and application in the optimal dispatch of IES due to their outstanding optimization performance. With their efficient solving capabilities and flexible adaptability, these algorithms have become essential tools for addressing complex optimization problems. Wang et al. developed an economic dispatch model for an electricity-heat-gas coupled regional IES incorporating a tiered carbon trading mechanism. By integrating carbon emission calculations with the carbon trading mechanism, the model incorporates carbon trading costs into the total system cost. The implementation of the fruit fly optimization algorithm demonstrates effective reduction in system operational costs [13].Considering the uncertainties of renewable energy resources, Jamal Raheela et al. proposed a scenario reduction technique based on an enhanced artificial hummingbird algorithm integrated with Monte Carlo simulation. This method generates an appropriate number of scenarios through Monte Carlo simulation reduction and employs the enhanced artificial hummingbird algorithm to solve both stochastic and deterministic optimal reactive power dispatch (ORPD) problems. Simulation results demonstrate that the enhanced artificial hummingbird algorithm effectively reduces active power losses, improves voltage profiles, and enhances voltage stability when addressing ORPD challenges [14]. Kou et al. addressed the susceptibility of wind turbines to failure in harsh offshore wind farm environments by proposing a chaotic simulated annealing genetic algorithm (CSAGA) incorporating asymmetric time considerations. The improved algorithm integrates logistic-tent chaotic mapping with genetic algorithms and simulated annealing techniques. Comparative simulation experiments demonstrated the superior performance of CSAGA over alternative algorithms [15]. Guo et al. developed a novel multi-strategy adaptive artificial bee colony algorithm (MSABC) for solving complex optimization problems. The MSABC incorporates an evolutionary rate index and an innovative elite-guided search strategy. Experimental results comparing MSABC with existing algorithms confirmed its enhanced convergence performance [16]. Ramadan Ashraf et al. accounted for uncertainties in load demand and renewable distributed generator (RDG) output power by proposing an efficient state-of-the-art technique for optimal RDG sizing and placement in radial distribution systems. The implementation of the artificial hummingbird algorithm (AHA) effectively mitigated power losses, reduced voltage deviation, and minimized expected costs [17]. However, single-objective optimization algorithms, constrained by their singular optimal solution output, inherently obscure the intrinsic conflicts and synergies among multidimensional objectives (e.g., economic efficiency and environmental sustainability) while oversimplifying the complexity of multi-objective synergistic optimization. This simplification risks convergence to local optima, hindering globally efficient resource allocation across multidimensional systems. In contrast, multi-objective optimization algorithms overcome these limitations via parallel optimization mechanisms, simultaneously addressing competing objectives without predefined weightings. By generating a Pareto-optimal solution set, they explicitly quantify trade-off relationships between goals such as economic viability and low-carbon performance, thereby eliminating biases from artificial intervention. The diverse non-dominated solutions within the set empower decision-makers with adaptable strategies spanning cost-prioritized to emission-reduction-focused approaches, accommodating dynamic scenario-specific demands. Furthermore, their global search mechanisms circumvent local optima traps, enabling synergistic resource utilization within multidimensional trade-off spaces. Owing to their superior alignment with complex and evolving optimization paradigms, multi-objective optimization algorithms have garnered substantial research attention in recent years [18].Shan et al. developed an energy optimization model that successfully achieves a globally optimal configuration using an enhanced multi-objective gray wolf algorithm [19]. Wu et al. proposed a grid-connected IES that considers the complementarity between geothermal energy, solar energy, and heat storage. A multi-objective non-dominated sorting genetic algorithm was introduced to optimize the system, improving operational efficiency [20]. Wang et al. constructed a dual-level optimization model for regional IES, considering both the quantity and quality of energy. By embedding a tabu search algorithm into a multi-objective genetic algorithm, they were able to effectively solve the model and significantly reduce the system’s economic cost [21]. Li et al. addressed complex models with multi-objective, non-convex, and strongly constrained decision variables by proposing a multi-objective whale optimization algorithm to solve a two-stage robust game model [22]. In summary, multi-objective optimization algorithms offer a critical pathway for addressing the synergistic optimization of multidimensional objectives in integrated energy systems. By generating a Pareto-optimal solution set, these algorithms explicitly quantify trade-off relationships among competing objectives, providing decision-makers with diverse strategy options and significantly enhancing system flexibility in dynamic scenarios. However, existing algorithms still face challenges such as local optima entrapment and slow convergence rates, particularly in complex systems with high-dimensional, nonlinear constraints, where traditional multi-objective methods often exhibit insufficient global search capabilities, leading to degraded solution quality, prolonged optimization time, and compromised real-time scheduling efficacy. To address these limitations, this study proposes a multi-strategy enhanced multi-objective artificial hummingbird algorithm for low-carbon economic dispatch. The model integrates adaptive spiral migration foraging and simplex method strategies to simultaneously avoid local optima and accelerate convergence. The research contributions are threefold: (1). A low-carbon economic dispatch optimization model is developed, minimizing total system costs while balancing economic and carbon emission objectives through the coordinated operation of carbon capture, refined power-to-gas processes, and a stepwise carbon trading mechanism. (2). A multi-strategy enhanced multi-objective artificial hummingbird algorithm is proposed, integrating logistic-sine fused chaotic mapping, elite opposition-based learning, adaptive spiral migration foraging, and simplex method strategies. The logistic-sine fused chaotic mapping and elite opposition-based learning enhance population initialization by improving the uniformity of initial populations and the quality of initial solutions, while adaptive spiral migration foraging and the simplex method refine overall population quality by updating positions of individuals with the lowest fitness, thereby strengthening the algorithm’s local search capability and optimization efficiency. (3). Simulation experiments validate the effectiveness of the proposed method. The improved multi-objective artificial hummingbird algorithm effectively addresses the challenges of local optima traps and slow convergence in high-dimensional objective spaces when applied to complex multi-dimensional problems in integrated energy systems, providing a reliable technical solution for achieving the “dual-carbon” targets in such systems. 2. Optimal scheduling model of integrated energy system in industrial park 2.1 Integrated energy system model The comprehensive energy system can effectively integrate and coordinate the operation of a variety of energy facilities to achieve synergistic and complementary power supply, so as to enhance the economic benefits of the entire energy system, but also better adapt to the diversified and multi-form energy consumption needs of modern society. In an Integrated Energy System, when a Combined Cooling, Heating, and Power unit equipped with carbon capture technology operates, the byproduct CO₂ is supplied as feedstock to Power-to-Gas equipment. The hydrogen produced by electrolyzers reacts with CO₂ through methanation to generate CH₄, which is then delivered to gas loads via the natural gas pipeline. The CCHP unit, serving as a coupling device for electrical, thermal, and cooling subsystems, provides electrical, heating, and cooling loads to users within the system. The thermal energy generated by electric boilers is utilized to meet the thermal load demand of users.The cold energy generated by electric refrigeration is used to meet the cooling load requirements of users. The integrated energy system model of industrial park is shown in Fig 1. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. Industrial park integrated energy system model. https://doi.org/10.1371/journal.pone.0325310.g001 2.2 Model of combined cooling, heating and electricity supply unit Combined cooling, heating, and power (CCHP) is a high-efficiency energy system that generates electricity via gas turbines while recovering waste heat from the power generation process for space heating or to drive absorption chillers for cooling, thereby achieving tri-generation of electricity, heat, and cold, with its internal configuration illustrated in Fig 2. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Internal structure of combined cooling, heat and electricity supply. https://doi.org/10.1371/journal.pone.0325310.g002 The model of the combined cooling, heating, and power supply unit is presented in Equation (1). (1) Where is the heat power output of the waste-heat boiler at time t; is the heat loss coefficient; is the electrical power output of the CCHP unit; is the power generation efficiency of the CCHP unit; is the waste-heat loss efficiency; is the cold power output of the absorption refrigerator is the CCHP refrigeration efficiency; is the heat production efficiency of the CCHP; is the gas power consumed by the CCHP unit. 2.3 Carbon capture model Carbon capture systems are engineered installations that extract CO₂emissions from industrial exhaust streams or energy production processes, preventing direct atmospheric release through integrated compression followed by either geological sequestration or conversion into industrial feedstocks, thereby achieving substantial greenhouse gas emission reductions; in the power generation sector, the integrated application of carbon capture and storage (CCS) technology plays a pivotal role in enhancing decarbonization effectiveness, making the retrofitting of conventional power plants with carbon capture systems a critical development trend in coming years [23], while the captured CO₂can additionally serve as a carbon source for power-to-gas (P2G) conversion processes [24]. A carbon capture device is installed in the CCHP unit. The electrical power output of the CCHP unit is composed of two parts: (1) the power consumption of the carbon capture device, , and (2) the remaining power consumption, referred to as the net output power of the CCHP unit, . The total electrical power output of the CCHP unit at time t is given by Equation (2). (2) The total energy consumption power of the carbon capture device consists of two parts: (1) the basic energy consumption power , which can be considered a constant, and (2) the operational energy consumption power , which is related to the amount of CO2 captured during operation. This relationship is represented by Equation (3). (3) Where is the carbon emission of the CCHP unit; is the energy consumption for CO₂ capture; is the carbon capture efficiency; is the unit carbon emission intensity; and is the amount of CO2 captured by the CCHP unit. 2.4 Power to gas two-stage operation model Power-to-gas (P2G) systems, serving as critical bidirectional coupling interfaces between electrical grids and natural gas networks, have garnered increasing research attention in IES due to their system flexibility enhancement potential [25]. The two-stage P2G conversion process enables efficient transformation of surplus renewable electricity into chemically storable and pipeline-transportable synthetic natural gas through sequential technological pathways: (1) hydrogen production via proton exchange membrane electrolysis and (2) subsequent catalytic methanation by reacting the produced hydrogen with captured carbon dioxide [26], with the complete technical schematic of this process illustrated in Fig 3. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. P2G two-stage operation process. https://doi.org/10.1371/journal.pone.0325310.g003 The electrolyzer converts surplus electrical energy into hydrogen through water electrolysis, with its mathematical model represented by Equation (4). (4) Where is the hydrogen generation power of the EL device; is the EL conversion efficiency; is the electrical power consumption of the EL device. The hydrogen fuel cell (HFC) is an electrochemical device that directly converts the chemical energy of hydrogen and oxygen into both electrical energy and thermal energy, with its mathematical formulation presented in Equation (5). (5) Where are the output thermal power of the HFC and the electrical power output of the HFC, respectively; is the input hydrogen power for the HFC; , are the efficiencies of heat and electricity conversion, respectively. The methanation reactor (MR) catalytically synthesizes methane through the Sabatier reaction between hydrogen produced from water electrolysis and captured carbon dioxide, with its governing equations mathematically formulated in Equation (6). (6) Where is the gas power output of MR; is MR conversion efficiency; is MR input hydrogen power; is the amount of CO2 consumed; is the density of CO2; and are the volumes of CO2 consumed and methane produced, respectively; is the conversion coefficient of electricity generation into heat supply; and is the calorific value of natural gas combustion. 2.5 Electric refrigerator model The model for the electric refrigerator (ER) is given in Equation (7). (7) Where is the output cold power of the ER; is the refrigeration coefficient of the ER; is the electrical power consumed by the ER. 2.6 Electric boiler model The electric boiler model is given by Equation (8). (8) Where is the output thermal power of the electric boiler; is the conversion efficiency of the electric boiler; and is the electric power consumed by the electric boiler. 2.7 Energy storage model The models of different types of energy storage devices in an integrated energy system are similar. In this paper, the two types of energy storage devices are unified, as shown in Equation (9). (9) Where represent the energy storage capacity of the i type energy storage device at times t and t-1; is the energy loss rate of the i type energy storage device; are the charging and discharging efficiencies of the i type energy storage device; are the charging and discharging power of the i type energy storage device at time t. 2.1 Integrated energy system model The comprehensive energy system can effectively integrate and coordinate the operation of a variety of energy facilities to achieve synergistic and complementary power supply, so as to enhance the economic benefits of the entire energy system, but also better adapt to the diversified and multi-form energy consumption needs of modern society. In an Integrated Energy System, when a Combined Cooling, Heating, and Power unit equipped with carbon capture technology operates, the byproduct CO₂ is supplied as feedstock to Power-to-Gas equipment. The hydrogen produced by electrolyzers reacts with CO₂ through methanation to generate CH₄, which is then delivered to gas loads via the natural gas pipeline. The CCHP unit, serving as a coupling device for electrical, thermal, and cooling subsystems, provides electrical, heating, and cooling loads to users within the system. The thermal energy generated by electric boilers is utilized to meet the thermal load demand of users.The cold energy generated by electric refrigeration is used to meet the cooling load requirements of users. The integrated energy system model of industrial park is shown in Fig 1. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. Industrial park integrated energy system model. https://doi.org/10.1371/journal.pone.0325310.g001 2.2 Model of combined cooling, heating and electricity supply unit Combined cooling, heating, and power (CCHP) is a high-efficiency energy system that generates electricity via gas turbines while recovering waste heat from the power generation process for space heating or to drive absorption chillers for cooling, thereby achieving tri-generation of electricity, heat, and cold, with its internal configuration illustrated in Fig 2. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Internal structure of combined cooling, heat and electricity supply. https://doi.org/10.1371/journal.pone.0325310.g002 The model of the combined cooling, heating, and power supply unit is presented in Equation (1). (1) Where is the heat power output of the waste-heat boiler at time t; is the heat loss coefficient; is the electrical power output of the CCHP unit; is the power generation efficiency of the CCHP unit; is the waste-heat loss efficiency; is the cold power output of the absorption refrigerator is the CCHP refrigeration efficiency; is the heat production efficiency of the CCHP; is the gas power consumed by the CCHP unit. 2.3 Carbon capture model Carbon capture systems are engineered installations that extract CO₂emissions from industrial exhaust streams or energy production processes, preventing direct atmospheric release through integrated compression followed by either geological sequestration or conversion into industrial feedstocks, thereby achieving substantial greenhouse gas emission reductions; in the power generation sector, the integrated application of carbon capture and storage (CCS) technology plays a pivotal role in enhancing decarbonization effectiveness, making the retrofitting of conventional power plants with carbon capture systems a critical development trend in coming years [23], while the captured CO₂can additionally serve as a carbon source for power-to-gas (P2G) conversion processes [24]. A carbon capture device is installed in the CCHP unit. The electrical power output of the CCHP unit is composed of two parts: (1) the power consumption of the carbon capture device, , and (2) the remaining power consumption, referred to as the net output power of the CCHP unit, . The total electrical power output of the CCHP unit at time t is given by Equation (2). (2) The total energy consumption power of the carbon capture device consists of two parts: (1) the basic energy consumption power , which can be considered a constant, and (2) the operational energy consumption power , which is related to the amount of CO2 captured during operation. This relationship is represented by Equation (3). (3) Where is the carbon emission of the CCHP unit; is the energy consumption for CO₂ capture; is the carbon capture efficiency; is the unit carbon emission intensity; and is the amount of CO2 captured by the CCHP unit. 2.4 Power to gas two-stage operation model Power-to-gas (P2G) systems, serving as critical bidirectional coupling interfaces between electrical grids and natural gas networks, have garnered increasing research attention in IES due to their system flexibility enhancement potential [25]. The two-stage P2G conversion process enables efficient transformation of surplus renewable electricity into chemically storable and pipeline-transportable synthetic natural gas through sequential technological pathways: (1) hydrogen production via proton exchange membrane electrolysis and (2) subsequent catalytic methanation by reacting the produced hydrogen with captured carbon dioxide [26], with the complete technical schematic of this process illustrated in Fig 3. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. P2G two-stage operation process. https://doi.org/10.1371/journal.pone.0325310.g003 The electrolyzer converts surplus electrical energy into hydrogen through water electrolysis, with its mathematical model represented by Equation (4). (4) Where is the hydrogen generation power of the EL device; is the EL conversion efficiency; is the electrical power consumption of the EL device. The hydrogen fuel cell (HFC) is an electrochemical device that directly converts the chemical energy of hydrogen and oxygen into both electrical energy and thermal energy, with its mathematical formulation presented in Equation (5). (5) Where are the output thermal power of the HFC and the electrical power output of the HFC, respectively; is the input hydrogen power for the HFC; , are the efficiencies of heat and electricity conversion, respectively. The methanation reactor (MR) catalytically synthesizes methane through the Sabatier reaction between hydrogen produced from water electrolysis and captured carbon dioxide, with its governing equations mathematically formulated in Equation (6). (6) Where is the gas power output of MR; is MR conversion efficiency; is MR input hydrogen power; is the amount of CO2 consumed; is the density of CO2; and are the volumes of CO2 consumed and methane produced, respectively; is the conversion coefficient of electricity generation into heat supply; and is the calorific value of natural gas combustion. 2.5 Electric refrigerator model The model for the electric refrigerator (ER) is given in Equation (7). (7) Where is the output cold power of the ER; is the refrigeration coefficient of the ER; is the electrical power consumed by the ER. 2.6 Electric boiler model The electric boiler model is given by Equation (8). (8) Where is the output thermal power of the electric boiler; is the conversion efficiency of the electric boiler; and is the electric power consumed by the electric boiler. 2.7 Energy storage model The models of different types of energy storage devices in an integrated energy system are similar. In this paper, the two types of energy storage devices are unified, as shown in Equation (9). (9) Where represent the energy storage capacity of the i type energy storage device at times t and t-1; is the energy loss rate of the i type energy storage device; are the charging and discharging efficiencies of the i type energy storage device; are the charging and discharging power of the i type energy storage device at time t. 3. Modeling of cascade carbon trading mechanism To achieve the dual objectives of carbon neutrality and peak carbon emissions, carbon trading markets are undergoing rapid development, with two predominant trading mechanisms emerging: conventional carbon trading and stepped carbon trading. Comparative analyses demonstrate that stepped carbon trading mechanisms exhibit superior efficacy in both emissions control and mitigation promotion relative to traditional approaches [27]. The stepped mechanism establishes a nonlinear relationship between carbon emission costs and emission volumes through phased pricing tiers and dynamic quota adjustments, operating via a structured framework: (1) initial quota allocation based on sector-specific emission caps, followed by (2) mandatory purchase of additional allowances at progressively increasing price tiers when emissions exceed allocated quotas, thereby creating market-driven resource reallocation that reduces aggregate societal abatement costs [28]. The actual carbon emission model is given by Equation (10). (10) Where is the actual carbon emission of IES; is the actual carbon emission from the CCHP units; is the actual carbon emissions associated with electricity purchased from the grid; is the total amount of CO2 captured by the carbon capture devices; is the carbon emissions generated by the CCHP units; is the carbon emission per unit of electricity generation; is the electric power purchased from the grid; T represents a time period, typically 24 h. The carbon emission allowance quota model is shown in Equation (11). (11) Where is the carbon emission quota for the IES; is the initial free carbon emission quota for the CCHP units; is the initial free carbon emission quota for electricity purchased from the grid; is the carbon emission quota per unit of electricity generated by the CCHP unit; and is the carbon emission quota per unit of electricity purchased from the grid. Step carbon trading mechanism model The stepped carbon trading mechanism divides the carbon emission trading volume into several intervals. The carbon trading price increases progressively with higher carbon emissions. Additionally, a compensation coefficient is introduced so that when the carbon emission allowance exceeds actual emissions, the excess carbon emission allowance can be sold. This is represented by Equations (12) and (13). (12)(13) Where is the total carbon emissions in the stepped carbon trading mechanism; is the cost of step-type carbon trading; is the base price of the cascade carbon trading; is the compensation coefficient; is the length of each carbon emission interval in the cascade carbon trading; is the growth rate of the step carbon trading price. 4. Objective functions and constraints 4.1 Objective function In this paper, we define the comprehensive cost of the system, denoted as , and carbon emissions, represented by , as our objective functions. The comprehensive cost includes various components such as electricity purchase costs (), gas purchase costs (), and operational expenses (), as illustrated in Equation (14). (14) The cost of purchasing electricity is given by Equation (15). (15) Where is the electricity purchase price; is the amount of electricity purchased from the grid. The gas purchase cost is expressed in Equation (16). (16) Where is the gas purchase price; is the amount of gas purchased. The operating costs are shown in Equation (17). (17) Where , are the operation and maintenance cost coefficients for the respective equipment. The detailed equipment parameters are provided in Table 1. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. Equipment parameters. https://doi.org/10.1371/journal.pone.0325310.t001 4.2 Constraint condition The electrical power balance constraint is given by Equation (18). (18) Where are the electrical power output from wind power and the electrical power output from photovoltaic systems, respectively; is the electrical power demand; are the energy storage charging power and the energy storage discharging power, respectively. The thermal power balance constraint is given by Equation (19). (19) Where is the thermal power demand. The cold power balance constraint is given by Equation (20). (20) Where is the cold power demand. The gas power balance constraint is given by Equation (21). (21) Where is the gas power demand. The hydrogen power balance constraints are given by Equation (22). (22) Where and are the energy storage charging power of the hydrogen storage device and the energy storage discharging power of the hydrogen storage device, respectively. The constraints for the CCHP units are shown in Equation (23). (23) Where , , and are the upper limit of the electrical, thermal, and cold power output of the CCHP unit, respectively. The constraints for the ER are given by Equation (24). (24) Where is the upper limit of the cold power output from the ER. The constraints for the CCS are shown in Equation (25). (25) Where is the upper limit of the electrical power consumption of the CCS. The constraint for the EL is shown in Equation (26). (26) Where is the upper limit of the electrical power consumed by the EL. The MR constraint is given by Equation (27). (27) Where is the upper limit of the hydrogen power input for the MR. The constraints for the HFC are shown in Equation (28). (28) Where is the upper limit of the hydrogen power input for the HFC. The EB constraint is shown in Equation (29). (29) Where is the upper limit of the thermal power output of the EB. The operation constraints for controllable units are given by Equation (30). (30) Where is the output of the i controlled unit at time t; and are the minimum and maximum output limits of the i controlled unit; and are the lower and upper limits of the climbing speed for the i controlled unit. The energy storage constraints are given by Equation (31). (31) Where and are the minimum and maximum charge states of the i type energy storage device; is the total capacity of the i type energy storage device; and are the maximum charging and discharging efficiencies of the i type energy storage device; and are the lower and upper power limits of the i type energy storage device; and are the state variables of energy storage and discharge for the i type energy storage device at time t; and are the lower and upper limits of the capacity of the i type energy storage device. 4.1 Objective function In this paper, we define the comprehensive cost of the system, denoted as , and carbon emissions, represented by , as our objective functions. The comprehensive cost includes various components such as electricity purchase costs (), gas purchase costs (), and operational expenses (), as illustrated in Equation (14). (14) The cost of purchasing electricity is given by Equation (15). (15) Where is the electricity purchase price; is the amount of electricity purchased from the grid. The gas purchase cost is expressed in Equation (16). (16) Where is the gas purchase price; is the amount of gas purchased. The operating costs are shown in Equation (17). (17) Where , are the operation and maintenance cost coefficients for the respective equipment. The detailed equipment parameters are provided in Table 1. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. Equipment parameters. https://doi.org/10.1371/journal.pone.0325310.t001 4.2 Constraint condition The electrical power balance constraint is given by Equation (18). (18) Where are the electrical power output from wind power and the electrical power output from photovoltaic systems, respectively; is the electrical power demand; are the energy storage charging power and the energy storage discharging power, respectively. The thermal power balance constraint is given by Equation (19). (19) Where is the thermal power demand. The cold power balance constraint is given by Equation (20). (20) Where is the cold power demand. The gas power balance constraint is given by Equation (21). (21) Where is the gas power demand. The hydrogen power balance constraints are given by Equation (22). (22) Where and are the energy storage charging power of the hydrogen storage device and the energy storage discharging power of the hydrogen storage device, respectively. The constraints for the CCHP units are shown in Equation (23). (23) Where , , and are the upper limit of the electrical, thermal, and cold power output of the CCHP unit, respectively. The constraints for the ER are given by Equation (24). (24) Where is the upper limit of the cold power output from the ER. The constraints for the CCS are shown in Equation (25). (25) Where is the upper limit of the electrical power consumption of the CCS. The constraint for the EL is shown in Equation (26). (26) Where is the upper limit of the electrical power consumed by the EL. The MR constraint is given by Equation (27). (27) Where is the upper limit of the hydrogen power input for the MR. The constraints for the HFC are shown in Equation (28). (28) Where is the upper limit of the hydrogen power input for the HFC. The EB constraint is shown in Equation (29). (29) Where is the upper limit of the thermal power output of the EB. The operation constraints for controllable units are given by Equation (30). (30) Where is the output of the i controlled unit at time t; and are the minimum and maximum output limits of the i controlled unit; and are the lower and upper limits of the climbing speed for the i controlled unit. The energy storage constraints are given by Equation (31). (31) Where and are the minimum and maximum charge states of the i type energy storage device; is the total capacity of the i type energy storage device; and are the maximum charging and discharging efficiencies of the i type energy storage device; and are the lower and upper power limits of the i type energy storage device; and are the state variables of energy storage and discharge for the i type energy storage device at time t; and are the lower and upper limits of the capacity of the i type energy storage device. 5. Improved multi-objective artificial hummingbird algorithm 5.1 Artificial hummingbird algorithm The Artificial Hummingbird Algorithm (AHA), proposed by Zhao et al. (2022) as a novel metaheuristic optimization technique, mimics the unique flight capabilities and intelligent foraging strategies observed in natural hummingbirds to achieve superior search performance [29]. Distinguished from conventional metaheuristics, AHA demonstrates three salient characteristics: (1) minimal parameter requirements, (2) stable convergence rates, and (3) enhanced optimization capacity, with its bio-inspired framework – particularly the emulation of hummingbirds’ specialized search and nectar exploration mechanisms – conferring distinct advantages in solving specific optimization challenges [30]. The algorithm computationally replicates hummingbird foraging behavior through an innovative visitation table system that implements memory functions by: 1) recording food source quality grades, 2) determining visitation priority levels based on elapsed time since last visit (longer intervals yielding higher priority), and 3) selecting targets with either the highest priority level or, when levels are equal, the maximal nectar replenishment rate (quantified by current fitness values). Initial population. Place N hummingbirds on each food source, and initialize their positions as defined in Equation (32): (32) Where represents the position of the i food source;and are the upper and lower limits of the search space respectively; is a random number between [0,1]. The food source access table is initialized via Equation (33). (33) Where i = j indicates that hummingbirds forage for specific food sources; i ≠ j indicates that the current j food source was visited by the i hummingbird. During flight, hummingbirds exhibit three flight modes: diagonal, axial, and omnidirectional. Each individual selects among these modes with equal probability. Since the flight process is not the main research content of this paper, its specific model will not be introduced in detail. Guided foraging. When rand < 0.5, hummingbirds conduct guided foraging according to their current flight skills, and obtain candidate food sources by visiting target food sources. The candidate food source location is updated as defined in Equation (34): (34) Where represents the position of the m-th candidate food source in t + 1 iteration; represents the position of the m-th candidate food source in the t iteration; indicates the location where the m-th hummingbird will visit the target food source; is a guiding foraging factor that follows the standard normal distribution. The food source location is updated via Equation (35): (35) Where represents the fitness value of the objective function. Regional feeding. When rand≥0.5, hummingbirds do regional foraging in areas adjacent to the current food source, with location updates defined by Equation (36). (36) Where is the regional foraging factor that follows the standard normal distribution. Migrate for food. AHA defines the migration coefficient M. When the number of iterations t reaches the migration coefficient M, hummingbirds migrate for food and replace the existing food source with other new food sources with the worst nectar replenishment rate. The updated food source locations are defined by Equation (37). (37) Where represents the location of the food source with the worst nectar replenishment rate in the population. 5.2 Multi-target artificial hummingbird algorithm In the multi-objective optimization problem, one objective is optimized at the expense of other objectives. In multi-objective problems, it is impossible to obtain a set of solutions with the minimum of each objective, but only a set of Pareto solutions can be obtained. The screening mechanism of MOAHA is as follows: MOAHA introduces an external archive to store a fixed number of non-dominated solutions; The congestion distance method based on dynamic elimination is adopted to maintain the diversity of solutions in the whole optimization process; Adopt the solution update mechanism based on non-dominant sorting to update the scheme. 5.3 Improved multi-objective artificial hummingbird algorithm 5.3.1 Integrating logistic-sine chaotic mapping with elite reverse learning to initialize the population. Most heuristic optimization algorithms are highly dependent on the initial population position. A randomly distributed initial population often leads to uneven distribution, which can adversely affect the algorithm’s accuracy. The AHA also uses a random population initialization, resulting in an uneven distribution of individuals and a lack of population diversity. To address this, chaotic mapping is frequently used to improve population initialization in meta-heuristic algorithms [31]. Among the various chaotic models, Tent and Logistic are the most commonly used. Xia et al. applied logistic-sine chaotic mapping to generate the initial population, enhancing both the uniformity of the initial population and the quality of the initial solutions [32]. In this paper, we propose a chaotic reverse learning strategy based on a combination of logistic-sine chaotic mapping and elite reverse learning to initialize the population. The mathematical expression of logistic-sine chaotic mapping is given by Equation (38). (38) Where is the chaos coefficient, which is set in this paper =0.5. First, N D-dimensional initial solutions (i = 1,2,…,N;j = 1,2,…,D) are generated by logistic-sine chaotic mapping. Secondly, the initial solutions are sorted, and the corresponding extreme point for each individual is selected as the elite individual =(,,...,). The chaos elite reverse solution =(,,...,) is then generated according to Equation (39). Where the dynamic boundary is set by Equation (40) to regulate the position points that cross the boundary: (39)(40) In the formula, the elite reverse coefficients ∈(0,1), =min() and =max() Finally, all the initial solutions and elite reverse solutions generated by the chaotic mapping are combined and sorted, and the best N solutions are selected as the initial population.The combined strategy integrating Logistic-Sine chaotic mapping with elite opposition-based learning for population initialization effectively addresses the limitations of conventional stochastic initialization methods. This synergistic approach prevents uneven distribution of population individuals while enhancing spatial coverage characteristics, thus mitigating premature diversity loss that critically impacts evolutionary optimization algorithms. The dual mechanisms ensure sufficient initialization heterogeneity through chaotic ergodicity and directional guidance from elite exemplars, establishing a robust foundation for subsequent global exploration phases. Figs 4 and 5 illustrate the process of integrating logistic-sine chaotic mapping with elite reverse learning for population initialization. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Integrating logistic-sine chaotic mapping with elite reverse learning to initialize the population. https://doi.org/10.1371/journal.pone.0325310.g004 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Integrating logistic-sine chaotic mapping with elite reverse learning to initialize the population. https://doi.org/10.1371/journal.pone.0325310.g005 5.3.2 Adaptive spiral migration foraging. In the AHA, migratory foraging behavior drives the least-fit hummingbird individuals to move randomly in search of better solutions. Inspired by the spiral search concept in the Whale Optimization Algorithm (WOA) [33], we propose an improvement to the migration-foraging behavior. The improved mathematical expression for migration and foraging is shown in Equation (41). (41) Where is the adaptive inertia weight; is an adaptive spiral parameter; and is a random number ranging from 0 to 1. is the maximum number of iterations. At the beginning of the iteration, is smaller, meaning that the influence of the optimal hummingbird on the worst-positioned hummingbird is limited. At the same time, is larger, which results in a larger spiral shape, allowing for a more thorough exploration of the solution space. In the later stages of the iterations, becomes larger, increasing the influence of the best hummingbirds on the worst hummingbirds. Concurrently, decreases, reducing the size of the spiral and causing the worst-positioned hummingbird to move closer to the optimal position, thereby accelerating the convergence rate.In the adaptive spiral migration-foraging strategy, when relocating the hummingbird individual with the poorest fitness to a new position for stochastic search during the migration-foraging process, the algorithm initiates with a large spiral pattern to facilitate global exploration. Subsequently, it gradually reduces the spiral radius to enable local refinement. This dual-phase mechanism effectively guides the inferior individuals toward the optimal solution region while maintaining population diversity, thereby accelerating convergence rates through balanced exploration-exploitation dynamics. Fig 6 shows adaptive spiral migration foraging. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 6. Adaptive spiral migration foraging. https://doi.org/10.1371/journal.pone.0325310.g006 5.3.3 Simplex method. The Simplex method is a polyhedral search algorithm known for its fast search speed and simple principles. It primarily determines whether the direction vector of the worst vertex g moves towards the best solution through iteration. The movement is controlled by reflecting, expanding, contracting, and reducing the search space. The Simplex method is depicted in Fig 7. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 7. Simplex method. https://doi.org/10.1371/journal.pone.0325310.g007 The Simplex method is implemented to examine the best advantage , the secondary advantage , and the maximum disadvantage , and subsequently calculate the position of the center point . A reflection operation is performed on the worst vertex , as shown in Equation (42). (42) Where is the reflection coefficient, and = 1. If >, the reflection direction is correct. In this case, an expansion operation is performed, as shown in Equation (43). (43) Where is the expansion coefficient, and = 2. If >, the expansion point replaces the maximum disadvantage point . Otherwise, the reflection point replaces the maximum disadvantage point . If > and the reflection direction is incorrect, an external contraction operation is performed, as shown in Equation (44). (44) Where is the external contraction coefficient, and = 0.5. If >, the external contraction point replaces the maximum disadvantage point . If <<, an internal contraction operation is performed, as shown in Equation (45). (45) Where is the internal shrinkage coefficient, and = 0.5. If >, the internal contraction point replaces the maximum disadvantage point . Otherwise, the reflection point replaces the maximum disadvantage point . By incorporating the Nelder-Mead simplex method into the iterative framework, this methodology systematically enhances population quality through three core operations. The reflection operation allows the worst-performing individual to explore all feasible solutions while avoiding premature convergence. Inner/outer contraction operations facilitate escaping from suboptimal positions while maintaining solution space coverage. Most critically, the expansion operation empowers elite solutions to escape local minima by conducting counter-directional searches away from the worst individual’s position. This systematic approach not only strengthens local search capability but also balances exploitation-exploration dynamics, ultimately improving global optimization performance. Fig 8 shows the procedure of the simplex method. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 8. Simplex method. https://doi.org/10.1371/journal.pone.0325310.g008 The flowchart for the improved MOAHA is shown in Fig 9. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 9. Flow chart of improved multi-objective artificial hummingbird algorithm. https://doi.org/10.1371/journal.pone.0325310.g009 5.4 Algorithm performance verification In order to verify the performance of the proposed algorithm, MOP test function is solved by using the proposed algorithm, MOAHA algorithm, MOGWO algorithm and MODA algorithm respectively. The parameters of the four algorithms are the same, the population size is 100, the number of iterations is 500, and the external archive is 100. Figs 10–13 show the simulation results. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 10. Simulation results of MOP standard test function ZDT1. https://doi.org/10.1371/journal.pone.0325310.g010 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 11. Simulation results of MOP standard test function ZDT2. https://doi.org/10.1371/journal.pone.0325310.g011 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 12. Simulation results of MOP standard test function ZDT3. https://doi.org/10.1371/journal.pone.0325310.g012 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 13. Simulation results of MOP standard test function ZDT6. https://doi.org/10.1371/journal.pone.0325310.g013 The performance of the algorithm is evaluated using two metrics: the IGD index and the Spacing index. These metrics are used to analyze the Pareto frontiers obtained from the four multi-objective problem-solving methods discussed above. The IGD index measures the distance between the non-dominated solution set produced by the algorithm and the actual Pareto frontier. It serves as an indicator of both the convergence and diversity of the algorithm’s performance. Specifically, the IGD index is calculated as the average Euclidean distance from each point on the true Pareto front to the closest point in the non-dominated set. A lower IGD value suggests that the generated solution set is more aligned with the true Pareto frontier. The mathematical expression for the IGD index is provided in Equation (46). (46) Where P is the Pareto solution set obtained by the algorithm; is the optimal solution set of the multi-objective problem; is the symbol for calculating the average value. The Spacing index evaluates the average distance between solutions within the non-dominated set generated by the algorithm, providing insight into the solution distribution. This index is determined by calculating the standard deviation of the minimum distances from each solution to the others in the set. A smaller Spacing value indicates that the solutions are more closely packed and more uniformly distributed. The formula for calculating the Spacing index is given by Equation (47). (47) Where is the average of all ; denotes the shortest distance from the ith solution to all other solutions in the solution space. The performance evaluation metrics for the four algorithms were calculated, and the results are shown in Table 2. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. Algorithm performance evaluation index results. https://doi.org/10.1371/journal.pone.0325310.t002 As summarized in Table 2, in the multi-objective optimization problem represented by the four test functions, the IGD and Spacing indices for the proposed algorithm (MOAHA) are significantly better than those for the other three algorithms. This indicates that the solution set obtained by the MOAHA is not only closer to the real Pareto frontier but also more evenly distributed and diverse compared to the solutions obtained by the other algorithms. In conclusion, the results demonstrate that the improved MOAHA proposed in this paper is highly effective in solving multi-objective optimization problems. 5.1 Artificial hummingbird algorithm The Artificial Hummingbird Algorithm (AHA), proposed by Zhao et al. (2022) as a novel metaheuristic optimization technique, mimics the unique flight capabilities and intelligent foraging strategies observed in natural hummingbirds to achieve superior search performance [29]. Distinguished from conventional metaheuristics, AHA demonstrates three salient characteristics: (1) minimal parameter requirements, (2) stable convergence rates, and (3) enhanced optimization capacity, with its bio-inspired framework – particularly the emulation of hummingbirds’ specialized search and nectar exploration mechanisms – conferring distinct advantages in solving specific optimization challenges [30]. The algorithm computationally replicates hummingbird foraging behavior through an innovative visitation table system that implements memory functions by: 1) recording food source quality grades, 2) determining visitation priority levels based on elapsed time since last visit (longer intervals yielding higher priority), and 3) selecting targets with either the highest priority level or, when levels are equal, the maximal nectar replenishment rate (quantified by current fitness values). Initial population. Place N hummingbirds on each food source, and initialize their positions as defined in Equation (32): (32) Where represents the position of the i food source;and are the upper and lower limits of the search space respectively; is a random number between [0,1]. The food source access table is initialized via Equation (33). (33) Where i = j indicates that hummingbirds forage for specific food sources; i ≠ j indicates that the current j food source was visited by the i hummingbird. During flight, hummingbirds exhibit three flight modes: diagonal, axial, and omnidirectional. Each individual selects among these modes with equal probability. Since the flight process is not the main research content of this paper, its specific model will not be introduced in detail. Guided foraging. When rand < 0.5, hummingbirds conduct guided foraging according to their current flight skills, and obtain candidate food sources by visiting target food sources. The candidate food source location is updated as defined in Equation (34): (34) Where represents the position of the m-th candidate food source in t + 1 iteration; represents the position of the m-th candidate food source in the t iteration; indicates the location where the m-th hummingbird will visit the target food source; is a guiding foraging factor that follows the standard normal distribution. The food source location is updated via Equation (35): (35) Where represents the fitness value of the objective function. Regional feeding. When rand≥0.5, hummingbirds do regional foraging in areas adjacent to the current food source, with location updates defined by Equation (36). (36) Where is the regional foraging factor that follows the standard normal distribution. Migrate for food. AHA defines the migration coefficient M. When the number of iterations t reaches the migration coefficient M, hummingbirds migrate for food and replace the existing food source with other new food sources with the worst nectar replenishment rate. The updated food source locations are defined by Equation (37). (37) Where represents the location of the food source with the worst nectar replenishment rate in the population. Initial population. Place N hummingbirds on each food source, and initialize their positions as defined in Equation (32): (32) Where represents the position of the i food source;and are the upper and lower limits of the search space respectively; is a random number between [0,1]. The food source access table is initialized via Equation (33). (33) Where i = j indicates that hummingbirds forage for specific food sources; i ≠ j indicates that the current j food source was visited by the i hummingbird. During flight, hummingbirds exhibit three flight modes: diagonal, axial, and omnidirectional. Each individual selects among these modes with equal probability. Since the flight process is not the main research content of this paper, its specific model will not be introduced in detail. Guided foraging. When rand < 0.5, hummingbirds conduct guided foraging according to their current flight skills, and obtain candidate food sources by visiting target food sources. The candidate food source location is updated as defined in Equation (34): (34) Where represents the position of the m-th candidate food source in t + 1 iteration; represents the position of the m-th candidate food source in the t iteration; indicates the location where the m-th hummingbird will visit the target food source; is a guiding foraging factor that follows the standard normal distribution. The food source location is updated via Equation (35): (35) Where represents the fitness value of the objective function. Regional feeding. When rand≥0.5, hummingbirds do regional foraging in areas adjacent to the current food source, with location updates defined by Equation (36). (36) Where is the regional foraging factor that follows the standard normal distribution. Migrate for food. AHA defines the migration coefficient M. When the number of iterations t reaches the migration coefficient M, hummingbirds migrate for food and replace the existing food source with other new food sources with the worst nectar replenishment rate. The updated food source locations are defined by Equation (37). (37) Where represents the location of the food source with the worst nectar replenishment rate in the population. 5.2 Multi-target artificial hummingbird algorithm In the multi-objective optimization problem, one objective is optimized at the expense of other objectives. In multi-objective problems, it is impossible to obtain a set of solutions with the minimum of each objective, but only a set of Pareto solutions can be obtained. The screening mechanism of MOAHA is as follows: MOAHA introduces an external archive to store a fixed number of non-dominated solutions; The congestion distance method based on dynamic elimination is adopted to maintain the diversity of solutions in the whole optimization process; Adopt the solution update mechanism based on non-dominant sorting to update the scheme. 5.3 Improved multi-objective artificial hummingbird algorithm 5.3.1 Integrating logistic-sine chaotic mapping with elite reverse learning to initialize the population. Most heuristic optimization algorithms are highly dependent on the initial population position. A randomly distributed initial population often leads to uneven distribution, which can adversely affect the algorithm’s accuracy. The AHA also uses a random population initialization, resulting in an uneven distribution of individuals and a lack of population diversity. To address this, chaotic mapping is frequently used to improve population initialization in meta-heuristic algorithms [31]. Among the various chaotic models, Tent and Logistic are the most commonly used. Xia et al. applied logistic-sine chaotic mapping to generate the initial population, enhancing both the uniformity of the initial population and the quality of the initial solutions [32]. In this paper, we propose a chaotic reverse learning strategy based on a combination of logistic-sine chaotic mapping and elite reverse learning to initialize the population. The mathematical expression of logistic-sine chaotic mapping is given by Equation (38). (38) Where is the chaos coefficient, which is set in this paper =0.5. First, N D-dimensional initial solutions (i = 1,2,…,N;j = 1,2,…,D) are generated by logistic-sine chaotic mapping. Secondly, the initial solutions are sorted, and the corresponding extreme point for each individual is selected as the elite individual =(,,...,). The chaos elite reverse solution =(,,...,) is then generated according to Equation (39). Where the dynamic boundary is set by Equation (40) to regulate the position points that cross the boundary: (39)(40) In the formula, the elite reverse coefficients ∈(0,1), =min() and =max() Finally, all the initial solutions and elite reverse solutions generated by the chaotic mapping are combined and sorted, and the best N solutions are selected as the initial population.The combined strategy integrating Logistic-Sine chaotic mapping with elite opposition-based learning for population initialization effectively addresses the limitations of conventional stochastic initialization methods. This synergistic approach prevents uneven distribution of population individuals while enhancing spatial coverage characteristics, thus mitigating premature diversity loss that critically impacts evolutionary optimization algorithms. The dual mechanisms ensure sufficient initialization heterogeneity through chaotic ergodicity and directional guidance from elite exemplars, establishing a robust foundation for subsequent global exploration phases. Figs 4 and 5 illustrate the process of integrating logistic-sine chaotic mapping with elite reverse learning for population initialization. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Integrating logistic-sine chaotic mapping with elite reverse learning to initialize the population. https://doi.org/10.1371/journal.pone.0325310.g004 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Integrating logistic-sine chaotic mapping with elite reverse learning to initialize the population. https://doi.org/10.1371/journal.pone.0325310.g005 5.3.2 Adaptive spiral migration foraging. In the AHA, migratory foraging behavior drives the least-fit hummingbird individuals to move randomly in search of better solutions. Inspired by the spiral search concept in the Whale Optimization Algorithm (WOA) [33], we propose an improvement to the migration-foraging behavior. The improved mathematical expression for migration and foraging is shown in Equation (41). (41) Where is the adaptive inertia weight; is an adaptive spiral parameter; and is a random number ranging from 0 to 1. is the maximum number of iterations. At the beginning of the iteration, is smaller, meaning that the influence of the optimal hummingbird on the worst-positioned hummingbird is limited. At the same time, is larger, which results in a larger spiral shape, allowing for a more thorough exploration of the solution space. In the later stages of the iterations, becomes larger, increasing the influence of the best hummingbirds on the worst hummingbirds. Concurrently, decreases, reducing the size of the spiral and causing the worst-positioned hummingbird to move closer to the optimal position, thereby accelerating the convergence rate.In the adaptive spiral migration-foraging strategy, when relocating the hummingbird individual with the poorest fitness to a new position for stochastic search during the migration-foraging process, the algorithm initiates with a large spiral pattern to facilitate global exploration. Subsequently, it gradually reduces the spiral radius to enable local refinement. This dual-phase mechanism effectively guides the inferior individuals toward the optimal solution region while maintaining population diversity, thereby accelerating convergence rates through balanced exploration-exploitation dynamics. Fig 6 shows adaptive spiral migration foraging. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 6. Adaptive spiral migration foraging. https://doi.org/10.1371/journal.pone.0325310.g006 5.3.3 Simplex method. The Simplex method is a polyhedral search algorithm known for its fast search speed and simple principles. It primarily determines whether the direction vector of the worst vertex g moves towards the best solution through iteration. The movement is controlled by reflecting, expanding, contracting, and reducing the search space. The Simplex method is depicted in Fig 7. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 7. Simplex method. https://doi.org/10.1371/journal.pone.0325310.g007 The Simplex method is implemented to examine the best advantage , the secondary advantage , and the maximum disadvantage , and subsequently calculate the position of the center point . A reflection operation is performed on the worst vertex , as shown in Equation (42). (42) Where is the reflection coefficient, and = 1. If >, the reflection direction is correct. In this case, an expansion operation is performed, as shown in Equation (43). (43) Where is the expansion coefficient, and = 2. If >, the expansion point replaces the maximum disadvantage point . Otherwise, the reflection point replaces the maximum disadvantage point . If > and the reflection direction is incorrect, an external contraction operation is performed, as shown in Equation (44). (44) Where is the external contraction coefficient, and = 0.5. If >, the external contraction point replaces the maximum disadvantage point . If <<, an internal contraction operation is performed, as shown in Equation (45). (45) Where is the internal shrinkage coefficient, and = 0.5. If >, the internal contraction point replaces the maximum disadvantage point . Otherwise, the reflection point replaces the maximum disadvantage point . By incorporating the Nelder-Mead simplex method into the iterative framework, this methodology systematically enhances population quality through three core operations. The reflection operation allows the worst-performing individual to explore all feasible solutions while avoiding premature convergence. Inner/outer contraction operations facilitate escaping from suboptimal positions while maintaining solution space coverage. Most critically, the expansion operation empowers elite solutions to escape local minima by conducting counter-directional searches away from the worst individual’s position. This systematic approach not only strengthens local search capability but also balances exploitation-exploration dynamics, ultimately improving global optimization performance. Fig 8 shows the procedure of the simplex method. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 8. Simplex method. https://doi.org/10.1371/journal.pone.0325310.g008 The flowchart for the improved MOAHA is shown in Fig 9. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 9. Flow chart of improved multi-objective artificial hummingbird algorithm. https://doi.org/10.1371/journal.pone.0325310.g009 5.3.1 Integrating logistic-sine chaotic mapping with elite reverse learning to initialize the population. Most heuristic optimization algorithms are highly dependent on the initial population position. A randomly distributed initial population often leads to uneven distribution, which can adversely affect the algorithm’s accuracy. The AHA also uses a random population initialization, resulting in an uneven distribution of individuals and a lack of population diversity. To address this, chaotic mapping is frequently used to improve population initialization in meta-heuristic algorithms [31]. Among the various chaotic models, Tent and Logistic are the most commonly used. Xia et al. applied logistic-sine chaotic mapping to generate the initial population, enhancing both the uniformity of the initial population and the quality of the initial solutions [32]. In this paper, we propose a chaotic reverse learning strategy based on a combination of logistic-sine chaotic mapping and elite reverse learning to initialize the population. The mathematical expression of logistic-sine chaotic mapping is given by Equation (38). (38) Where is the chaos coefficient, which is set in this paper =0.5. First, N D-dimensional initial solutions (i = 1,2,…,N;j = 1,2,…,D) are generated by logistic-sine chaotic mapping. Secondly, the initial solutions are sorted, and the corresponding extreme point for each individual is selected as the elite individual =(,,...,). The chaos elite reverse solution =(,,...,) is then generated according to Equation (39). Where the dynamic boundary is set by Equation (40) to regulate the position points that cross the boundary: (39)(40) In the formula, the elite reverse coefficients ∈(0,1), =min() and =max() Finally, all the initial solutions and elite reverse solutions generated by the chaotic mapping are combined and sorted, and the best N solutions are selected as the initial population.The combined strategy integrating Logistic-Sine chaotic mapping with elite opposition-based learning for population initialization effectively addresses the limitations of conventional stochastic initialization methods. This synergistic approach prevents uneven distribution of population individuals while enhancing spatial coverage characteristics, thus mitigating premature diversity loss that critically impacts evolutionary optimization algorithms. The dual mechanisms ensure sufficient initialization heterogeneity through chaotic ergodicity and directional guidance from elite exemplars, establishing a robust foundation for subsequent global exploration phases. Figs 4 and 5 illustrate the process of integrating logistic-sine chaotic mapping with elite reverse learning for population initialization. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Integrating logistic-sine chaotic mapping with elite reverse learning to initialize the population. https://doi.org/10.1371/journal.pone.0325310.g004 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Integrating logistic-sine chaotic mapping with elite reverse learning to initialize the population. https://doi.org/10.1371/journal.pone.0325310.g005 5.3.2 Adaptive spiral migration foraging. In the AHA, migratory foraging behavior drives the least-fit hummingbird individuals to move randomly in search of better solutions. Inspired by the spiral search concept in the Whale Optimization Algorithm (WOA) [33], we propose an improvement to the migration-foraging behavior. The improved mathematical expression for migration and foraging is shown in Equation (41). (41) Where is the adaptive inertia weight; is an adaptive spiral parameter; and is a random number ranging from 0 to 1. is the maximum number of iterations. At the beginning of the iteration, is smaller, meaning that the influence of the optimal hummingbird on the worst-positioned hummingbird is limited. At the same time, is larger, which results in a larger spiral shape, allowing for a more thorough exploration of the solution space. In the later stages of the iterations, becomes larger, increasing the influence of the best hummingbirds on the worst hummingbirds. Concurrently, decreases, reducing the size of the spiral and causing the worst-positioned hummingbird to move closer to the optimal position, thereby accelerating the convergence rate.In the adaptive spiral migration-foraging strategy, when relocating the hummingbird individual with the poorest fitness to a new position for stochastic search during the migration-foraging process, the algorithm initiates with a large spiral pattern to facilitate global exploration. Subsequently, it gradually reduces the spiral radius to enable local refinement. This dual-phase mechanism effectively guides the inferior individuals toward the optimal solution region while maintaining population diversity, thereby accelerating convergence rates through balanced exploration-exploitation dynamics. Fig 6 shows adaptive spiral migration foraging. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 6. Adaptive spiral migration foraging. https://doi.org/10.1371/journal.pone.0325310.g006 5.3.3 Simplex method. The Simplex method is a polyhedral search algorithm known for its fast search speed and simple principles. It primarily determines whether the direction vector of the worst vertex g moves towards the best solution through iteration. The movement is controlled by reflecting, expanding, contracting, and reducing the search space. The Simplex method is depicted in Fig 7. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 7. Simplex method. https://doi.org/10.1371/journal.pone.0325310.g007 The Simplex method is implemented to examine the best advantage , the secondary advantage , and the maximum disadvantage , and subsequently calculate the position of the center point . A reflection operation is performed on the worst vertex , as shown in Equation (42). (42) Where is the reflection coefficient, and = 1. If >, the reflection direction is correct. In this case, an expansion operation is performed, as shown in Equation (43). (43) Where is the expansion coefficient, and = 2. If >, the expansion point replaces the maximum disadvantage point . Otherwise, the reflection point replaces the maximum disadvantage point . If > and the reflection direction is incorrect, an external contraction operation is performed, as shown in Equation (44). (44) Where is the external contraction coefficient, and = 0.5. If >, the external contraction point replaces the maximum disadvantage point . If <<, an internal contraction operation is performed, as shown in Equation (45). (45) Where is the internal shrinkage coefficient, and = 0.5. If >, the internal contraction point replaces the maximum disadvantage point . Otherwise, the reflection point replaces the maximum disadvantage point . By incorporating the Nelder-Mead simplex method into the iterative framework, this methodology systematically enhances population quality through three core operations. The reflection operation allows the worst-performing individual to explore all feasible solutions while avoiding premature convergence. Inner/outer contraction operations facilitate escaping from suboptimal positions while maintaining solution space coverage. Most critically, the expansion operation empowers elite solutions to escape local minima by conducting counter-directional searches away from the worst individual’s position. This systematic approach not only strengthens local search capability but also balances exploitation-exploration dynamics, ultimately improving global optimization performance. Fig 8 shows the procedure of the simplex method. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 8. Simplex method. https://doi.org/10.1371/journal.pone.0325310.g008 The flowchart for the improved MOAHA is shown in Fig 9. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 9. Flow chart of improved multi-objective artificial hummingbird algorithm. https://doi.org/10.1371/journal.pone.0325310.g009 5.4 Algorithm performance verification In order to verify the performance of the proposed algorithm, MOP test function is solved by using the proposed algorithm, MOAHA algorithm, MOGWO algorithm and MODA algorithm respectively. The parameters of the four algorithms are the same, the population size is 100, the number of iterations is 500, and the external archive is 100. Figs 10–13 show the simulation results. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 10. Simulation results of MOP standard test function ZDT1. https://doi.org/10.1371/journal.pone.0325310.g010 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 11. Simulation results of MOP standard test function ZDT2. https://doi.org/10.1371/journal.pone.0325310.g011 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 12. Simulation results of MOP standard test function ZDT3. https://doi.org/10.1371/journal.pone.0325310.g012 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 13. Simulation results of MOP standard test function ZDT6. https://doi.org/10.1371/journal.pone.0325310.g013 The performance of the algorithm is evaluated using two metrics: the IGD index and the Spacing index. These metrics are used to analyze the Pareto frontiers obtained from the four multi-objective problem-solving methods discussed above. The IGD index measures the distance between the non-dominated solution set produced by the algorithm and the actual Pareto frontier. It serves as an indicator of both the convergence and diversity of the algorithm’s performance. Specifically, the IGD index is calculated as the average Euclidean distance from each point on the true Pareto front to the closest point in the non-dominated set. A lower IGD value suggests that the generated solution set is more aligned with the true Pareto frontier. The mathematical expression for the IGD index is provided in Equation (46). (46) Where P is the Pareto solution set obtained by the algorithm; is the optimal solution set of the multi-objective problem; is the symbol for calculating the average value. The Spacing index evaluates the average distance between solutions within the non-dominated set generated by the algorithm, providing insight into the solution distribution. This index is determined by calculating the standard deviation of the minimum distances from each solution to the others in the set. A smaller Spacing value indicates that the solutions are more closely packed and more uniformly distributed. The formula for calculating the Spacing index is given by Equation (47). (47) Where is the average of all ; denotes the shortest distance from the ith solution to all other solutions in the solution space. The performance evaluation metrics for the four algorithms were calculated, and the results are shown in Table 2. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. Algorithm performance evaluation index results. https://doi.org/10.1371/journal.pone.0325310.t002 As summarized in Table 2, in the multi-objective optimization problem represented by the four test functions, the IGD and Spacing indices for the proposed algorithm (MOAHA) are significantly better than those for the other three algorithms. This indicates that the solution set obtained by the MOAHA is not only closer to the real Pareto frontier but also more evenly distributed and diverse compared to the solutions obtained by the other algorithms. In conclusion, the results demonstrate that the improved MOAHA proposed in this paper is highly effective in solving multi-objective optimization problems. 6. Simulation example analysis In order to verify the effectiveness of the proposed algorithm in solving the multi-objective optimization problem of the integrated energy system, the integrated energy system of an industrial park in North China is taken as an example for simulation analysis. Electricity price peak period (11:00–14:00,18:00–21:00) is 0.83 yuan/KWH; The normal period (7:00–10:00,15:00–17:00) is 0.51 yuan/KWH; The valley time (1:00–6:00,22:00–24:00) is 0.3 yuan/KWH; The price of natural gas is 3.5 yuan/m3. The load of electricity, heat, cold and gas and the output of new energy are shown in Fig 14. The parameters of each unit in the system are shown in Table 3. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 3. System parameters of each unit. https://doi.org/10.1371/journal.pone.0325310.t003 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 14. Electricity, heat, cold, gas load and new energy output. https://doi.org/10.1371/journal.pone.0325310.g014 6.1 Model optimization and solution To evaluate the performance of the proposed model in terms of carbon emission reduction and cost efficiency, three distinct scenarios are analyzed: Scenario 1: The model does not consider the carbon capture device but includes P2G equipment and the stepped carbon trading mechanism; Scenario 2: The model does not consider the carbon capture device but adopts the refined two-stage P2G operation and the stepped carbon trading mechanism; Scenario 3: This model incorporates the carbon capture device, refined two-stage P2G operation, and the stepped carbon trading mechanism, representing the optimization model proposed in this study; Scenario 3 represents the proposed optimization model, with its optimized operational outcomes illustrated in Figs 15–19 and solved using the improved multi-objective artificial hummingbird algorithm for case studies under a 24-hour operational cycle, where comparative analysis of four algorithms yields the Pareto front shown in Fig 20, configured with a population size of 100, 500 iterations, and an external archive capacity of 50. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 15. Electric power optimization results. https://doi.org/10.1371/journal.pone.0325310.g015 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 16. Thermal power optimization results. https://doi.org/10.1371/journal.pone.0325310.g016 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 17. Cold power optimization results. https://doi.org/10.1371/journal.pone.0325310.g017 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 18. Gas power optimization results. https://doi.org/10.1371/journal.pone.0325310.g018 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 19. Hydrogen power optimization results. https://doi.org/10.1371/journal.pone.0325310.g019 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 20. Four Pareto frontiers of algorithms. https://doi.org/10.1371/journal.pone.0325310.g020 6.2 Comparative analysis of different scenarios Fig 20 depicts the Pareto frontiers of the proposed integrated energy system model generated by four comparative algorithms. Simulation results indicate that the developed methodology achieves superior Pareto-optimal solutions in multi-objective scheduling optimization for integrated energy systems, outperforming benchmark approaches in simultaneously balancing economic efficiency and carbon emission mitigation. The simulation results for the three scenarios are presented in Table 4. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 4. Results of different scenarios. https://doi.org/10.1371/journal.pone.0325310.t004 Compared to Scenario 1, where the two-stage P2G operation process is absent, the addition of this process in Scenario 2 (replacing the basic P2G equipment) results in a 19.22% increase in the operation and maintenance costs. However, by generating H2 through the electrolyzer, part of which is supplied to the HFC for electricity and heat generation, and part of which is used by the MR to produce natural gas, the system reduces overall energy loss. This leads to a 4.03% reduction in the comprehensive system cost, demonstrating that incorporating the detailed two-stage P2G operation process improves the system’s economic efficiency. In the context of an Integrated Energy System, the introduction of carbon capture equipment in Scenario 3, compared to Scenario 2, led to a 5.11% increase in operation and maintenance costs. However, this addition resulted in a significant reduction of 10,680.27 kg in carbon emissions, corresponding to a 79.13% decrease. Moreover, by utilizing the captured CO₂ as a feedstock for the methanation process, the overall system cost was reduced by 13.83%. This demonstrates that incorporating carbon capture equipment into an IES can yield substantial economic and environmental benefits. 6.1 Model optimization and solution To evaluate the performance of the proposed model in terms of carbon emission reduction and cost efficiency, three distinct scenarios are analyzed: Scenario 1: The model does not consider the carbon capture device but includes P2G equipment and the stepped carbon trading mechanism; Scenario 2: The model does not consider the carbon capture device but adopts the refined two-stage P2G operation and the stepped carbon trading mechanism; Scenario 3: This model incorporates the carbon capture device, refined two-stage P2G operation, and the stepped carbon trading mechanism, representing the optimization model proposed in this study; Scenario 3 represents the proposed optimization model, with its optimized operational outcomes illustrated in Figs 15–19 and solved using the improved multi-objective artificial hummingbird algorithm for case studies under a 24-hour operational cycle, where comparative analysis of four algorithms yields the Pareto front shown in Fig 20, configured with a population size of 100, 500 iterations, and an external archive capacity of 50. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 15. Electric power optimization results. https://doi.org/10.1371/journal.pone.0325310.g015 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 16. Thermal power optimization results. https://doi.org/10.1371/journal.pone.0325310.g016 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 17. Cold power optimization results. https://doi.org/10.1371/journal.pone.0325310.g017 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 18. Gas power optimization results. https://doi.org/10.1371/journal.pone.0325310.g018 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 19. Hydrogen power optimization results. https://doi.org/10.1371/journal.pone.0325310.g019 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 20. Four Pareto frontiers of algorithms. https://doi.org/10.1371/journal.pone.0325310.g020 6.2 Comparative analysis of different scenarios Fig 20 depicts the Pareto frontiers of the proposed integrated energy system model generated by four comparative algorithms. Simulation results indicate that the developed methodology achieves superior Pareto-optimal solutions in multi-objective scheduling optimization for integrated energy systems, outperforming benchmark approaches in simultaneously balancing economic efficiency and carbon emission mitigation. The simulation results for the three scenarios are presented in Table 4. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 4. Results of different scenarios. https://doi.org/10.1371/journal.pone.0325310.t004 Compared to Scenario 1, where the two-stage P2G operation process is absent, the addition of this process in Scenario 2 (replacing the basic P2G equipment) results in a 19.22% increase in the operation and maintenance costs. However, by generating H2 through the electrolyzer, part of which is supplied to the HFC for electricity and heat generation, and part of which is used by the MR to produce natural gas, the system reduces overall energy loss. This leads to a 4.03% reduction in the comprehensive system cost, demonstrating that incorporating the detailed two-stage P2G operation process improves the system’s economic efficiency. In the context of an Integrated Energy System, the introduction of carbon capture equipment in Scenario 3, compared to Scenario 2, led to a 5.11% increase in operation and maintenance costs. However, this addition resulted in a significant reduction of 10,680.27 kg in carbon emissions, corresponding to a 79.13% decrease. Moreover, by utilizing the captured CO₂ as a feedstock for the methanation process, the overall system cost was reduced by 13.83%. This demonstrates that incorporating carbon capture equipment into an IES can yield substantial economic and environmental benefits. 7. Conclusion This paper develops a low-carbon economic optimization scheduling model for IES, integrating a carbon capture unit, a refined two-stage P2G process, and a stepped carbon trading mechanism. The model is optimized using the improved MOAHA, and the performance of three distinct scenarios is evaluated. The main conclusions are as follows: The synergistic integration of power-to-gas facilities, combined cooling-heating-power systems, and advanced energy storage devices establishes a multi-energy complementarity framework that significantly enhances energy utilization efficiency. This configuration minimizes energy cascade losses through coordinated multi-carrier energy conversion while enabling closed-loop carbon recycling via methanation processes. Thus, the IES low-carbon economic operation is improved. The integration of carbon capture facilities demonstrates significant economic and environmental benefits, achieving a system energy procurement cost reduction of 2,944.69 yuan (7.63% decrease rate) while simultaneously reducing carbon emissions by 10,680.27 kg, which consequently lowers carbon trading costs by 3,931.2 yuan.It shows that CCS equipment can better achieve the effect of reducing carbon emissions, so as to realize the low-carbon economic operation of IES. By incorporating carbon capture equipment and refining the two-stage operation process of Power-to-Gas (P2G), the carbon emissions were reduced by 13,664.08 kg compared to the scenario without these considerations. Additionally, the carbon trading costs were lowered by 5,050.13 yuan. As a result, the overall integrated cost was decreased by 8,639.14 yuan, representing a reduction of 17.3%.The refined two-phase power-to-gas operation mechanism, serving as an advanced alternative to conventional P2G systems, demonstrates enhanced performance in renewable energy accommodation while achieving significant reductions in operational costs and carbon emissions. This optimized approach utilizes CO₂ captured through carbon capture and storage technology as the primary carbon source for methanation processes, thus establishing a closed-loop carbon cycle. Furthermore, surplus hydrogen output can be effectively channeled to hydrogen fuel cells, with the subsequent electrical and thermal energy generated from HFC operation supplying a portion of the energy demand for integrated energy systems. This synergistic configuration establishes an energy cascade utilization framework that substantially improves hydrogen energy conversion efficiency across multiple operational phases. The enhanced multi-objective artificial hummingbird algorithm demonstrates superior capability in generating high-quality Pareto-optimal solutions for low-carbon economic dispatch problems in IES.This improved computational framework exhibits significant advantages in achieving optimal equilibrium between economic efficiency and environmental sustainability. The algorithm’s adaptive dynamic search strategy enables comprehensive exploration of the solution space while maintaining computational efficiency, thereby providing decision-makers with diversified optimal operational schemes that satisfy multiple conflicting objectives in IES management. Supporting information S1 Fig. All the pictures in the article. https://doi.org/10.1371/journal.pone.0325310.s001 (ZIP) S1 Table. All the tables in the article. https://doi.org/10.1371/journal.pone.0325310.s002 (DOCX) S1 Data. Original data. https://doi.org/10.1371/journal.pone.0325310.s003 (XLSX)
TI - Research on optimal scheduling of integrated energy system based on improved multi-objective artificial hummingbird algorithm
JF - PLoS ONE
DO - 10.1371/journal.pone.0325310
DA - 2025-06-04
UR - https://www.deepdyve.com/lp/public-library-of-science-plos-journal/research-on-optimal-scheduling-of-integrated-energy-system-based-on-0nBge0tyjw
SP - e0325310
VL - 20
IS - 6
DP - DeepDyve
ER -