# Multiobjective Optimized Dispatching for Integrated Energy System Based on Hierarchical Progressive Parallel NSGA-II Algorithm

Authors: Aidong Zeng; Sipeng Hao; Jia Ning; Qingshan Xu; Ling Jiang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Considering the importance of reducing system operating costs and controlling pollutant emissions by optimizing the operation of the integrated energy system, the energy supply structure of the integrated energy system and the joint multiobjective optimization dispatching structure is analyzed in this paper based on a day-ahead economic optimization dispatching model of the integrated energy system. Afterwards, the multiobjective optimization model of the integrated energy system is studied and multiobjective hierarchical progressive parallel algorithm based on improved NSGA-II is proposed according to the characteristics of the model. The algorithm improves the nondominated layer sorting algorithm, changes the convergence judgment condition while introducing the target reaching method to accelerate convergence, and introduces parallel computing technology according to the characteristics of the algorithm. The case shows that the proposed algorithm not only has advantages on the diversity in searching solutions but also can achieve better results in many aspects such as the iteration time and algorithm convergence which are required in practical engineering projects.


---

## Body

## 1. Introduction

Facing the pressure of energy shortage, energy saving, and emission reduction, building a clean, reliable, interactive, and efficient distributed cooling, heating, and power integrated energy system has become an important means to promote the innovation of efficient energy-use models and develop a low-carbon economy.Research and development of optimized dispatching and operation technologies for integrated heating and cooling power systems is a necessary condition for vigorously developing the efficient use of distributed energy systems and also the technical basis for carrying out a multienergy system demonstration project that uses electric power as the core and optimizes the operation of a variety of energy sources, which can provide support for realizing the transformation of energy utilization ways and a prerequisite for promoting changes in production and lifestyle, and becomes an important part of the country's carrying and promoting energy reform and the third industrial revolution.The optimal dispatching goal of the combined heat and power system is not only economical, commonly, the minimum exchange power with the bus and the minimum adjustment of equipment operating power are also available. In [1], an island-type integrated energy system in eastern China was taken as the research object, and the optimization goals were set as economical optimization and minimum power exchange with the tie bus. The operation scheme of the integrated energy system was obtained by solving the optimal scheduling model. In [2], an active distribution network (ADN) with integrated energy system was taken as the research object, and a method for managing the congestion of ADN through reasonable dispatch of integrated energy system was proposed. In the method, a new probability-based risk multiple objective model to achieve the optimal balance between feeder congestion and economic goals was defined. In [3], a comprehensive energy system optimization plan with the goal of optimizing the performance of the combined heat and power (CHP) plant was proposed. The optimization goals in the plan were to maximize power production, minimize fuel consumption, maximize the heat used in the slurry process, and improve the CHP plant’s power production efficiency.In recent years, a number of scholars have taken environmental protection factors into the system’s day-ahead scheduling model. For example, ultralow-energy consumption buildings were used as research objects in [4] to model their energy management and control their operating costs while optimizing their emission of greenhouse gases. In [5], a penalty function of polluting gas emissions was added to the system’s optimization objective function, and a multiobjective opportunity planning model was established by considering factors comprehensively such as the operating costs of the combined heat and power system and overall pollutant emissions. Finally, an improved particle swarm algorithm with local and global memory was adopted to solve the optimization model.In [6], HOMER Pro software was used to simulate the optimization of microgrid based on CHP plants. The microgrid architecture included PV arrays, diesel generators, and batteries (operating in grid-connected and island mode). Multiple operational goals included total net present cost (TNPC), cost of energy (COE), and annual greenhouse gas (GHG) emissions generated, as well as the maximization of annual waste heat recovery (WHR) and annual grid sales (GS). In [7], greenhouse gas emissions and pollutant discharge costs were incorporated into the system’s operating cost function, and a multiobjective optimization model was established with the lowest comprehensive operating cost of the integrated cooling and heating power system as the objective function, and the abovementioned model was solved by particle swarm optimization algorithm. In [8], a multiobjective thermodynamic model of a paper mill’s combined heat and power (CHP) system was established. In the multiobjective optimization model, three optimization goals were considered, including CHP power generation efficiency, total cost of system products, and CO2 of the entire plant emissions, and genetic algorithms are utilized to optimize the objective function.As mentioned above, optimal operation scheduling of an integrated energy system is no longer limited to the optimum economic dispatch; more factors such as sustainable development of environment, prevention of feeder congestion, and energy efficiency improvement need to be taken into account.In terms of optimization model solving methods, there are currently two main types of solving methods: one is to convert multiple optimization targets into a single optimization target and then perform single-objective optimization, and the other is to directly apply intelligent multiobjective algorithms for optimization.In the first method, for example, a multiobjective optimal scheduling model that considers both the economic and environmental objectives of a combined heat and power CHP system was proposed [9]. The optimization problem was transformed into a single-objective optimization problem by maximum fuzzy satisfying method, and then the model was solved by using an improved single-objective genetic algorithm. Finally, the validity of the proposed model and algorithm was verified by a multiobjective optimization example including a wind-solar-storage system and a CHP unit. Energy equivalent factor which converts three kinds of energy (cold, heat, and electricity) into each other equivalently was introduced [10], and a multiobjective energy-saving scheduling model for a combined heat and cold power system that includes system fuel costs, environmental emission costs, and energy coordination costs was established. Target membership function was used to convert multiple targets into single targets, and finally, quadratic programming algorithm was utilized to solve the model.To solve the shortcoming of lacking objective basis for the selection of weight factors among multiple objective functions, the negative ideal point method was proposed in [11], which introduced subjective and objective weighting methods to determine the weight factors among various objective functions. The algorithm was applied to a combined heat and power system with the core of a gas internal combustion engine. An equal performance coefficient from the perspective of pollutant emissions was proposed to reflect the energy consumption of the system in [12], and a multiobjective optimal scheduling model of the integrated cooling, heating, and power system was established. The optimization objectives of the model include power generation costs and emission costs. The two optimization goals were transformed into a single goal for processing by the membership function. Finally, the model was solved using a nonlinear programming algorithm. Distributed energy resources (DER) and microgrids were taken as research objects in [13], and an optimization model for optimal management of MG was proposed based on the combined heat and power (CHP) system considering economic, environmental, and reliability aspects. The exchange market algorithm (EMA) and weighted factors method were used to combine three conflicting targets and treat the multiobjective problem as a single target problem, and the fuzzy satisfaction method was applied to select the best compromise solution.The biggest problem in converting a multiobjective problem into a single-objective optimization is that it cannot fully find a set of solutions that meets the multiobjective optimization problem. The single target solution obtained is often a weight of multiple targets with less effective in balancing trade-offs or preferences for each goal. Therefore, in recent years, more and more researches have shifted the solution method to the direct use of intelligent algorithms such as genetic algorithms and particle swarm optimization to directly solve multiobjective problems and find the appropriate problem solution set.Among the multiobjective genetic algorithms, the widely used algorithms are nondominated sorting genetic algorithm (NSGA) and nondominated sorting genetic algorithm-II (NSGA-II) [14–16]. The sorting method of NSGA is simple and intuitive, that is, to traverse the entire population, compare the individuals one by one and then select individuals that meet the requirements until the layering work is completed in the cycle. In an m-dimensional multiobjective optimization problem, the time complexity of the sorting part is O(mN3) [17]. The NSGA is simple and easy to understand, but there are three major problems [18]: (1) it has high time complexity; (2) it does not adopt elite retention strategy; (3) it needs to specify the shared radius parameter in advance when adjusting the virtual fitness value. As an improved algorithm of NSGA, NSGA-II solved the abovementioned three problems by adopting better accounting and elite strategies which reduced the time complexity of the sorting algorithm from O(mN3) to O(mN2) [19]. NSGA-II reduces the calculation time compared with NSGA, but in some large optimization projects, there is still a long calculation time problem [20].A cold, heat, and power integrated energy system with the core of a small biogas internal combustion engine was taken as the research object in [21]. A multiobjective optimal scheduling model with comprehensive economic benefits, environmental benefits, and energy utilization benefits was established based on energy flow characteristics and operating constraints of the system. The system’s multiobjective optimization scheduling model is finally solved by using the PSO multiobjective algorithm with constraint factors to obtain the system’s Pareto solution set. Time-of-use electricity price was considered and a multiobjective optimal scheduling model of a combined heat and power system with the lowest fuel cost, the lowest pollutant emissions, and the best electricity purchase and sales from public networks was established in [22]. The optimal problem was solved by using a spatially coupled PSO optimization algorithm. An example showed that the algorithm had stronger global search capabilities than traditional PSO algorithms and can obtain more reliable Pareto solutions. Carbon emissions trading was taken as the research point and pollutant emission penalty function was stripped from the system operating costs in [23], then a multiobjective optimal scheduling model for a combined heat and power system that considers carbon emissions trading costs and system overall operating costs was established. A fuzzy self-correcting PSO algorithm was used to solve the model. The example showed that the system can effectively control the carbon emissions under the premise of ensuring normal operation and obtain some transaction benefits, thus improving the system’s economic operation level.A multiobjective optimization model of the CHP integrated energy system was established considering randomness of wind power and load fluctuations in the case of off-grid mode in [24]. The optimization objective considered operating costs, comprehensive benefits, and system power supply reliability. Finally, the model was solved by NSGA algorithm to obtain the optimal solution set. A multiobjective optimization model with the goal of minimizing system operating costs and carbide emissions for a CHP integrated energy system was established in [25]. Considering the nonlinearity and high complexity of the model, modified MOPSO algorithm based on the differential evolution algorithm was used to solve the model. Finally, a laboratory integrated energy system is used as a case to verify the effectiveness of the proposed model and algorithm. Economic and environmental objectives were considered and a dual-objective optimal scheduling model for a CHP system was proposed in [26]. The model was solved using a fuzzy MOPSO algorithm based on chance constraint programming (CCP). The Pareto optimal solution set and the optimal solution under two subobjectives are given in the simulation example.In summary, the research on the multiobjective optimization operation of the integrated energy system is mainly focused on the selection of multiple targets and the solution method of the multiobjective optimization problem. The single-objective transformation method had the disadvantage of a relatively limited solution set, the diversity of the searched solutions was insufficient, and it was easy to fall into a local optimum. However, using intelligent multiobjective algorithm to solve the problem has the problem that the solution process is lengthy and difficult to integrate into actual engineering applications.Therefore, based on the establishment of the optimal operation model of the integrated energy system with multiple optimization objectives, this paper proposes a multiobjective hierarchical progressive parallel NSGA-II algorithm according to the characteristics of the model. The algorithm improves the nondominated layer sorting algorithm, changes the convergence judgment condition to introduce the target reaching method to accelerate convergence, and introduces parallel computing technology according to the characteristics of the algorithm, which accelerates the convergence speed of the algorithm while taking into account the diversity of the solutions searched.The case study shows that the algorithm proposed in this paper not only has advantages on the diversity of solutions searched but also can achieve better results in many aspects such as the iteration time and algorithm convergence which are required in practical engineering projects. In practical engineering applications, the optimal operation scheme can be selected according to the needs of customers to achieve multiobjective optimized operation of the integrated energy system in the smart grid zone.


## 2. Energy Supply Structure of Integrated Energy System

The integrated energy system is a new type of overall energy supply system on the demand side with many advantages such as various forms of energy supply, numerous energy supply equipment, flexible energy supply schemes, and satisfactory overall economy. Generally, the energy sources of integrated energy systems include power grids, independent distributed power generation equipment, and gas companies. The loads supplied mainly include various industrial customers, residential customers, and commercial customers around the system. Various types of customers have different cold, heat, and electricity needs, and the integrated energy system can formulate corresponding energy supply plans based on these energy needs to meet customers’ load requirements.The main energy-supplying equipment of the integrated energy system studied in this paper includes a CHP unit, which contains microgas turbines, waste heat boilers, and absorption refrigerators. Refrigeration equipment includes electric refrigerator units, refrigeration air conditioners, and cold storage devices. Heating equipment includes gas boilers, heating air conditioners, and heat storage devices. Pure electrical supply equipment includes storage batteries and photovoltaic cells.The integrated energy system can freely trade electricity with the public grid, such as purchase electricity, from the grid when the system power supply is insufficient or the grid power price is lower than the system’s power generation cost and sell electricity to the grid when the system power supply is surplus. At the same time, it is supposed that there is no gas production in the integrated energy system and all the required gas is purchased from external gas companies. The energy supply structure that reflects the energy flow in the system is shown in Figure1.Figure 1
Energy supply structure of the integrated energy system.The integrated energy system performs bilateral energy interaction with the external power grid and unidirectional energy interaction with external natural gas companies. For the CHP unit, the microgas turbine converts the chemical energy of natural gas into electricity, and the waste smoke is utilized by a waste heat boiler to generate high-temperature steam, for heating or cooling through an absorption lithium bromide unit. In the cooling system, the electric refrigerator unit can generate cold water required for cooling, the cold storage device can flexibly store and release the cold water medium according to the needs, and the space cooling load can be met by cooling water through the fan coil cooling or directly through the refrigeration air conditioner. In the thermal system, the gas boiler can generate hot water required for heating and the heat storage device can flexibly store and release the hot water medium according to the needs. The space heat load can be met by hot water through the fan tray or directly by air conditioner. The solar and storage power generation unit in the integrated energy system can exchange power with the external power grid through the public bus.Based on this energy supply structure, this paper establishes a joint multiobjective optimization dispatching model for an integrated energy system. The input conditions for the dispatching model mainly include the next day forecast values of various loads in the system, TOU (time of use) price (electricity and gas) of the energy market, and various types of energy supply equipment models, capacities, and operating conditions in the system.A Pareto solution set is obtained by solving the joint multiobjective optimization model of the integrated energy system, and then the dispatching instructions are issued by one solution of the Pareto solution set according to the goal preference. The input and output of the model are shown in Figure2. The mathematical models of the energy supply equipment in the integrated energy system are established in [27].Figure 2
Joint multiobjective optimization dispatching structure of the integrated energy system.


## 3. Joint Multiobjective Optimization Dispatching Model of Integrated Energy System

The integrated energy system contains a variety of energy equipment; in order to satisfy various loads in the system efficiently, the energy equipment operation status in the integrated energy system must be optimized for dispatching according to the optimization goal.Based on the formulation of the energy supply structure of the integrated energy system and the establishment of mathematical models for various energy supply equipment in the system, a multidimensional optimized operation objective function is established considering the time-of-use energy price and various equipment operation constraints. By solving the multiobjective optimization dispatching model, the purpose of reasonably formulating an operation strategy that meets the requirements of the dispatching objective is achieved.The assumptions of the model are as follows:(1)
It is assumed that the output of the equipment is continuously adjustable within the operating constraints(2)
The equipment has no failures during the optimization operation cycle


### 3.1. Multidimensional Optimized Operation Objective Function

In the process of joint multiobjective operation optimization of the integrated energy system, not only the economic status of the system is considered but also the environmental indicators of the entire integrated energy system. The environmental indicators are also called emission target..(1)
The economic goal of the integrated energy system is to obtain the lowest daily operating cost which mainly includes three aspects: system fuel costs, power interaction costs with the grid, and system operation maintenance costs:(1)f1=minKTotal=minKFuel+KGrid+KMaintain,wheref1 is the economic goal of the system; KTotal is the daily comprehensive operation cost of the system; KFuel is the fuel costs of the system; KGrid is the power interaction costs with the grid; and KMaintain is the operation maintenance costs of the system,wherein the fuel costs function for fuel boilers and microturbines in the integrated energy system is calculated as(2)KFuel=∑t=124∑i=1nCHPcGastfCHPiPit+∑t=124∑i=1nGBicGastFGBit,wherefCHPi is the energy consumption curve function of ith microturbines; Pit is the electric power output of ith microturbines with the unit of kW; cGast is the hourly price of gas with the unit of $/(kW⋅h) after conversion according to the calorific value of natural gas; FGBit is the energy consumption of the ith gas boiler at time t with the unit of kW; and t is the time period number with the unit of hour.The power interaction costs function of the system with the grid is calculated as follows:(3)KGrid=∑t=124cGridtPGridt,wherecGridt is the hourly price of electricity and PGridt is the hourly power exchange value of the integrated energy system with external grid.The operation maintenance costs function of the system is calculated as(4)KMaintain=∑t=124pmWBQWBt+∑t=124pmACPACt+∑t=124pmECPchilt+∑t=124pmstorHint+∑t=124pmstorHoutt+∑t=124pmAirPcondt+∑t=124∑i=1nCHPpmCHPiPit+∑t=124∑i=1ndistripmdistriPdistrit+∑t=124pmGBQGBt,wherepmCHPi is the unit power operation maintenance cost of the micro turbine; pmdistri is the unit power operation maintenance cost of the distributed generation unit; pmstor is the unit power operation maintenance cost of the storage devices; pmGB is the unit power operation maintenance cost of the gas boiler; pmWB is the unit power operation maintenance cost of the waste heat boiler; pmAC is the unit power operation maintenance cost of the absorption chiller; pmEC is the unit power operation maintenance cost of the compression electric chiller; pmAir is the unit power operation maintenance cost of the air conditioner; pi is the electric power output of the ith micro turbine; Pdistrit is the output of the ith distributed generation unit; Hint and Houtt refer to the charging and discharging power of three types of energy storage equipment, respectively; QGBt is the heating power of the gas boiler; QWBt is the heating power of the waste heat boiler; PACt is the refrigeration power of the absorption chiller; Pchilt is the operating power of the compression electric chiller; Pcondt is the operating power of the air-conditioning system; and the abovementioned variable units have been converted into kW.(2)
The environmental protection goal of the integrated energy system is to obtain the lowest daily operating emissions which mainly include three aspects: nitride, carbide, and sulfide emissions.(5)f2=minPGTotal=minδ1×PGC+δ2×PGS+δ3×PGN,wheref2 is the environmental protection goal of the system; PGTotal is the total daily operating emissions of pollutants from the integrated energy system; PGC is the daily carbon emissions of the entire region; PGS is the daily sulfide emissions of the entire region; and PGN is the daily nitride emissions of the entire region and the specific formulas are as follows:(6)PGC=∑t=124∑i=1nGiCO+TiCO+ZiCO×Δt,PGS=∑t=124∑i=1nGiSO+TiSO+ZiSO×Δt,PGN=∑t=124∑i=1nGiNO+TiNO+ZiNO×Δt,whereGiCO, TiCO, and ZiCO are the COX emissions produced from the region’s purchase of electricity, microturbines, and gas boilers with the unit of kg/h; GiNO, TiNO, and ZiNO are the NOX emissions produced from the region’s purchase of electricity, microturbines, and gas boilers with the unit of kg/h; GiSO, TiSO, and ZiSO are the SOX emissions produced from the region’s purchase of electricity, microturbines, and gas boilers with the unit of kg/h; δ1, δ2, and δ3 are the weight coefficients of carbon emissions, sulfur emissions, and nitrogen emissions.


### 3.2. System Constraints

The joint multiobjective operation of the integrated energy system needs to meet three types of constraints: system operation constraints, environmental emission constraints, and economic constraints.(1)
Operational constraints of the system include electric power balance constraints, cooling and heating power supply constraints, equipment capacity constraints, and equipment operation constraints.The constraints of electric power balance are as follows:(7)∑i=1nCHPPit+∑i=1ndistriPdistrit+PGridt+Pstort=PLoadt+Pcondt,wherePGridt is the hourly power exchange value of the integrated energy system with public grid; PLoadt is the hourly electrical load; PStort is the charging or discharging power of battery; Pcondt is the hourly operation power of air conditioner; Pdistrit is the hourly operation power of distributed generation unit; and Pit is the power generated by the ith microturbine.The constraint of heating power supply is as follows:(8)∑i=1nCHPHit+∑i=1nboilerHboilerit+COPcond×Pcondt+Houtt−Hint≥Hspacet+HWatert.The constraint of hot water supply is as follows:(9)∑i=1nCHPHit+∑i=1nboilerHboilerit+Houtt−Hint≥HWatert,whereHit is the retrieved heat value by the ith microturbine through waste heat boiler; Hboilert is the heat output of the ith gas boiler; Hint and Houtt are the input and output of the thermal energy storage equipment; COPcond is the coefficient of performance of air-conditioning equipment; and Hspacet and Hwatert are the hourly space heating load and hot water load of the integrated energy system, respectively.The constraint of cooling power supply is as follows:(10)∑i=1nCHPCit+∑i=1nboilerCchilt+EERcond×Pcondt+Coutt−Cint≥Cspacet+CRefrit.The constraint of refrigeration power supply is as follows:(11)∑i=1nCHPCit+∑i=1nboilerCchilt+Coutt−Cint≥CRefrit,whereCit is the refrigeration power of the ith microturbine manufactured by the absorption refrigeration unit; Cchilt is the refrigeration power generated by the electric chiller; Cint and Coutt are the power input and output of the cold storage equipment, respectively; EERcond is the cooling energy efficiency ratio of air-conditioning equipment; and CSpacet and CRefrit are the hourly space cooling load and refrigeration load of the integrated energy system, respectively.The equipment capacity and operating constraints of the integrated energy system are as follows:for microturbines,(12)Pimin≤Pit≤Pimax,i∈nCHP.For gas boilers,(13)0≤Hboilerit≤Hboilerimax,i∈nboiler.For air-conditioning equipment,(14)0≤Pcondt≤Pcondmax.For heat (cold) energy storage equipment,(15)0≤Hint≤Hinmax,0≤Houtt≤Houtmax,Sstormin≤Sstort≤Sstormax,whereHint and Houtt are the power input and output of the heat (cold) storage equipment at time t; Hinmax and Houtmax are the power input and output limits of the heat (cold) storage equipment; and Sstort is the state of charge of the heat (cold) storage equipment at time t.(2)
During the operation of the integrated energy system, it is necessary to control the amount of pollutants emitted by the entire system within a certain range, including the total pollutant discharge constraints in the region and the hourly discharge constraints of various pollutants.Constraints on total carbon emissions in the region are as follows:(16)0≤PGC≤PGCmax.Constraints on total sulfide emissions in the region are as follows:(17)0≤PGS≤PGSmax.Constraints on total nitride emissions in the region are as follows:(18)0≤PGN≤PGNmax,wherePGCmax is the upper limit of the daily carbon emission value of the entire region; PGSmax is the upper limit of the daily sulfide emission value of the entire region; and PGNmax is the upper limit of the daily nitride emission value of the entire region.The hourly carbon emission constraints in the region are as follows:(19)0≤∑i=1nGiCO+TiCO+ZiCOt≤PGCt.The hourly sulfide emission constraints in the region are as follows:(20)0≤∑i=1nGiSO+TiSO+ZiSOt≤PGSt.The hourly nitride emission constraints in the region are as follows:(21)0≤∑i=1nGiNO+TiNO+ZiNOt≤PGNt,wherePGCt,PGSt, and PGNt are the upper limits of the carbon, sulfide, and nitride emission value of the entire region at time t.(3)
During the operation of the integrated energy system, the operating costs of the entire system needs to be controlled within a certain range. The economic constraint is as follows:(22)0≤KTotal≤Kmax,whereKTotal is the entire daily operating cost of the integrated energy system and Kmax is the upper limit of the daily operating cost.


### 3.3. Solving Algorithm

Due to the existence of a variety of energy supply equipment in the optimization model and the complex coupling relationship between cold, heat, and electrical energy, the model has a large number of optimization variables and a large volume of constraint conditions, when using traditional global solution space search algorithms such as NSGA-II or multiobjective particle swarm algorithm will make the solution process unacceptably lengthy which cannot meet the time limit requirements of the rolling optimization process within the dispatching system. When the traditional multiobjective optimization algorithm which transitions from single-objective optimization algorithm such as constraint method or weighting method is used to solve the model, the solution obtained is often limited to the local Pareto optimal solution set and it is difficult to obtain a large range of Pareto optimal solution set;, therefore, the traditional multiobjective optimization algorithm cannot provide the dispatching system with a variety of dispatching schemes for the system to make decisions.Based on the abovementioned analysis, this paper proposes a hierarchical progressive parallel NSGA-II algorithm suitable for solving multiobjective dispatching strategies of integrated energy systems. When solving the optimization model, it can not only make full use of the advantages of the global solution space search algorithm in the balance and diversity of the solution set searched but also achieve fast solution speed for large-scale nonlinear optimization problems which can meet the calculation time limit requirements of the rolling optimization process within the dispatching system. The algorithm has made three improvements based on NSGA-II as follows:(1)
Improve the nondominated layer sorting algorithm to achieve fast layering based on model characteristics.In the process of solving the joint multiobjective optimization operation model of the integrated energy system, considering that the model only includes dual optimization goals, this paper adopts a more efficient and fast nondominated layer sorting method according to the characteristics of the model to meet the requirements of rapid model solving.Firstly, the merge-sorting method is used to sort the individuals in the population by ascending order according to first objective function value. The merge sorting method is a sorting method by using the divide and conquer strategy. In divided stages, the method divides the problem into small questions and solves it recursively, while in the conquer stages, the method combines the answers obtained in divided stages together to get the results [28–30]. Suppose there is a population as shown in Figure 3 and each individual has two objective function values. The overall process of sorting the individuals in the population by ascending order according to first objective function value is shown in Figure 4. In the divided process, this structure is much like a complete binary tree and it can be understood as the process of recursively splitting subsequences. In the conquer process, we need to merge two already ordered subsequences into an ordered sequence, such as the last merge in Figure 4, to merge the already ordered subsequences of [4, 5, 7, 8] and [1–3, 6] into the final sequence [1–8]. The last merge within the conquer process demonstration is shown in Figure 5.If the first objective function values of the two individuals are the same, the individuals are sorted in ascending order according to the second objective function value and the result of the merge sorting is stored in the arrayFi. In this case, the individual with the smallest first objective function value must be an individual in the first layer of the nondominated layer. Then, the adjacent individual of the first individual is searched; if the second objective function value of the adjacent individual is smaller than the previous individual, the individual and the previous individual are belonging to the same layer of the nondominated layer; otherwise, the individual is not belonging to the same layer and will be inserted at the end of the array Fi and so on to find all first-level nondominated individuals in the population. Next, start with the first element at the end of the array Fi to find the second layer of nondominated individuals and so on until the population is completely stratified. In the example described above, the process of stratification occurs which is shown in Figure 6.Since the number of individuals in the next generation is half of that inFi, an incomplete layering strategy is adopted during the layering process to interrupt layering as soon as the number of individuals in the next generation exceeds the next generation population size, which can avoid many unnecessary layering comparison steps. The schematic diagram of the nondominated layer layering algorithm based on improved NSGA-II is shown in Figure 7:In the nondominated layer sorting algorithm, the average time complexity of the merge sorting process is O(N ∗ log2N) and the total times of calculations spent in the layering operation are(23)TTotal=∑i=1kTi=∑i=1kN−i=N−1+N−2+⋯+N−k<kN,whereN is twice the size of the population size; k is the total number of layers after the nondominated sorting is completed; TTotal is the total times of calculations spent in the layering operation after the nondominated sorting is completed; and Ti is the maximum times of comparisons required for generating the ith nondominated layer.The time complexity of the hierarchical (layering) operation is O (kN) after the nondominated sorting is completed and the total time complexity of the algorithm is O (N ∗ log2 N + kN). In general, k≪N, so the total time complexity of the algorithm is O (N ∗ log2 N). The sorting principle of NSGA-II is similar to the bubble sorting which needs to compare the individuals within the population in pairs and the average time complexity of the sorting algorithm is O (mN2), where m is the number of objectives. Compared with NSGA-II, the proposed algorithm reduces the time complexity from square order to linear logarithmic order.(2)
Modify the convergence judgment conditions and combine the target reaching method to accelerate the convergence process.The improved multiobjective genetic algorithm based on NSGA-II can quickly approach the vicinity of the Pareto front set in the process of solving the model. However, in order to achieve the final convergence, it requires a large number of iterations and function comparisons which means longer calculation time. Therefore, a comprehensive hierarchical progressive algorithm is proposed in this paper to speedup the process of finding the Pareto optimal solution set.Firstly, the aforementioned improved genetic algorithm based on NSGA-II is used for searching solutions close to the Pareto front by setting its termination condition to limited generations. Secondly, the target reaching method is utilized for a secondary search which uses the calculation result of the previous step as the initial value and a local search is performed by the target reaching method to obtain the Pareto solution set.The standard form of utilizing the target reaching method to solve the model is(24)minγs.t.fix−wi×γ≤Fi∗,cx≤0,ceqx=0,Ax≤b,Aeqx≤beq,lb≤x≤ub,where optimization variablex is an n-dimensional variable which includes hourly cold, heat, and power output of microturbine-based CHP system, hourly cold, heat, and power output of auxiliary energy supply equipment, hourly cold, heat, and power output of various energy storage equipment, SOC status of energy storage equipment, and hourly interacting power with the grid. fi(x) is the ith optimization objective function; ω is the weight coefficient of the ith objective function; Fi∗is the optimization target value expected to be achieved by the ith target; γ is the target achievement coefficient which is used for representing the achievement of the predetermined goal; c(x) and ceq(x) are the various nonlinear inequality constraints and nonlinear equality constraints met by the optimization variables; A and b are the coefficient matrices of various linear inequality constraints met by the optimization variables; Aeq and beq are the coefficient matrices of various linear equality constraints met by the optimization variables; lb and ub are the lower and upper limits of the optimization variable x, respectively.In the target reaching method, the multiobjective optimization problem is transformed into a single-objective optimization problem by using virtual variableγ and the single-objective optimization model is solved by using the interior point method. In the solving process, the interior point method constructs a Lagrangian auxiliary function which meets the Karush–Kuhn–Tucker conditional:(25)Lx,λ=fx+∑λc,icix+∑λA,iAx−b+∑λceq,iceqix+∑λAeq,iAeqx−beq.The second-order derivative Hessian matrix in its correction equation is(26)H=∇xx2Lx,λ=∇2fx+∑λi∇2cix+∑λi∇2ceqix,whereλc,iλA,iλceq,iλAeq,i are Lagrange multipliers which possess the same dimensions with respective constraints. x, lb, and ub are n × 1 column vectors; b is a p × 1 column vector; beq is a q × 1 column vector; A is a p × n matrix; Aeq is a q × n matrix; Hessian matrix is an n × n matrix.Solutions can be achieved by transformation, LDLT decomposition, iteration, and correction of the matrices in the correction equations. The iteration is automatically stopped when the convergence conditions are met and the solution closest to the expected optimization target in the feasible region is found as a Pareto optimal solution.(3)
Introduce parallel computing technology to accelerate the solution process according to the characteristics of the algorithm.Because there are a large number of repetitive loop iterative processes in the multiobjective hierarchical progressive algorithm based on improved NSGA-II, in addition, each process has no parameter passing to each other and does not depend on the results of other process, parallel computing technology is introduced to accelerate the convergence process of the algorithm. Among them, the part applicable for parallel computing mainly includes the nondominated layer merge sorting calculation and the call of the target reaching method for achieving fast convergence. In the nondominated layer merge sorting section, the fitness of individuals is first calculated on multiple parallel processing units, and then the individuals of the population to be sorted are evenly distributed tom parallel processing units. The complexity of the merge sorting process is reduced to O (N/m ∗ log2N) under parallel processing and the algorithm speeds up with the increase of the parallel processing units. When the target reaching method is utilized for a secondary search, the method starts from different initial values and there is no parameter passing to each other, so it is also suitable for using parallel processing technology to evenly distribute the computing task to multiple parallel processing units of the computer.In addition, because the GPU which is designed for parallel operations has more computing units than the CPU, in order to make full use of the GPU in large-scale array matrix processing, multiple operations without complex logical judgments and repeated calculations with little data correlation, the population crossover and mutation operations in the algorithm are assigned to the GPU. As the GPU is limited by the hardware architecture, single-precision calculations are much faster than double-precision calculations; therefore, this paper uses single-precision (single type variable) calculations for various types of population operations.Through the improvement of the abovementioned three aspects, the hierarchical progressive parallel NSGA-II algorithm can greatly reduce the calculation time in the process of obtaining the Pareto solution set, which can meet the requirements of the rolling optimization process within the dispatching system while ensuring the uniform distribution of the solution and the diversity of the population. The solution set obtained by the algorithm can provide flexible and diverse dispatching schemes for decision makers to choose.The hierarchical progressive parallel NSGA-II algorithm is written in MATLAB to solve the joint multiobjective optimization dispatching model of the integrated energy system in this paper. The algorithm flowchart is shown in Figure8.Figure 3
The individuals with two objective function values in a population.Figure 4
The overall process of the merge-sorting method.Figure 5
The last merge within the conquer process demonstration.Figure 6
The process of stratification in the example.Figure 7
Schematic diagram of nondominated layer layering algorithm based on improved NSGA-IIFigure 8
Flowchart of hierarchical progressive parallel NSGA-II algorithm.


## 3.1. Multidimensional Optimized Operation Objective Function

In the process of joint multiobjective operation optimization of the integrated energy system, not only the economic status of the system is considered but also the environmental indicators of the entire integrated energy system. The environmental indicators are also called emission target..(1)
The economic goal of the integrated energy system is to obtain the lowest daily operating cost which mainly includes three aspects: system fuel costs, power interaction costs with the grid, and system operation maintenance costs:(1)f1=minKTotal=minKFuel+KGrid+KMaintain,wheref1 is the economic goal of the system; KTotal is the daily comprehensive operation cost of the system; KFuel is the fuel costs of the system; KGrid is the power interaction costs with the grid; and KMaintain is the operation maintenance costs of the system,wherein the fuel costs function for fuel boilers and microturbines in the integrated energy system is calculated as(2)KFuel=∑t=124∑i=1nCHPcGastfCHPiPit+∑t=124∑i=1nGBicGastFGBit,wherefCHPi is the energy consumption curve function of ith microturbines; Pit is the electric power output of ith microturbines with the unit of kW; cGast is the hourly price of gas with the unit of $/(kW⋅h) after conversion according to the calorific value of natural gas; FGBit is the energy consumption of the ith gas boiler at time t with the unit of kW; and t is the time period number with the unit of hour.The power interaction costs function of the system with the grid is calculated as follows:(3)KGrid=∑t=124cGridtPGridt,wherecGridt is the hourly price of electricity and PGridt is the hourly power exchange value of the integrated energy system with external grid.The operation maintenance costs function of the system is calculated as(4)KMaintain=∑t=124pmWBQWBt+∑t=124pmACPACt+∑t=124pmECPchilt+∑t=124pmstorHint+∑t=124pmstorHoutt+∑t=124pmAirPcondt+∑t=124∑i=1nCHPpmCHPiPit+∑t=124∑i=1ndistripmdistriPdistrit+∑t=124pmGBQGBt,wherepmCHPi is the unit power operation maintenance cost of the micro turbine; pmdistri is the unit power operation maintenance cost of the distributed generation unit; pmstor is the unit power operation maintenance cost of the storage devices; pmGB is the unit power operation maintenance cost of the gas boiler; pmWB is the unit power operation maintenance cost of the waste heat boiler; pmAC is the unit power operation maintenance cost of the absorption chiller; pmEC is the unit power operation maintenance cost of the compression electric chiller; pmAir is the unit power operation maintenance cost of the air conditioner; pi is the electric power output of the ith micro turbine; Pdistrit is the output of the ith distributed generation unit; Hint and Houtt refer to the charging and discharging power of three types of energy storage equipment, respectively; QGBt is the heating power of the gas boiler; QWBt is the heating power of the waste heat boiler; PACt is the refrigeration power of the absorption chiller; Pchilt is the operating power of the compression electric chiller; Pcondt is the operating power of the air-conditioning system; and the abovementioned variable units have been converted into kW.(2)
The environmental protection goal of the integrated energy system is to obtain the lowest daily operating emissions which mainly include three aspects: nitride, carbide, and sulfide emissions.(5)f2=minPGTotal=minδ1×PGC+δ2×PGS+δ3×PGN,wheref2 is the environmental protection goal of the system; PGTotal is the total daily operating emissions of pollutants from the integrated energy system; PGC is the daily carbon emissions of the entire region; PGS is the daily sulfide emissions of the entire region; and PGN is the daily nitride emissions of the entire region and the specific formulas are as follows:(6)PGC=∑t=124∑i=1nGiCO+TiCO+ZiCO×Δt,PGS=∑t=124∑i=1nGiSO+TiSO+ZiSO×Δt,PGN=∑t=124∑i=1nGiNO+TiNO+ZiNO×Δt,whereGiCO, TiCO, and ZiCO are the COX emissions produced from the region’s purchase of electricity, microturbines, and gas boilers with the unit of kg/h; GiNO, TiNO, and ZiNO are the NOX emissions produced from the region’s purchase of electricity, microturbines, and gas boilers with the unit of kg/h; GiSO, TiSO, and ZiSO are the SOX emissions produced from the region’s purchase of electricity, microturbines, and gas boilers with the unit of kg/h; δ1, δ2, and δ3 are the weight coefficients of carbon emissions, sulfur emissions, and nitrogen emissions.


## 3.2. System Constraints

The joint multiobjective operation of the integrated energy system needs to meet three types of constraints: system operation constraints, environmental emission constraints, and economic constraints.(1)
Operational constraints of the system include electric power balance constraints, cooling and heating power supply constraints, equipment capacity constraints, and equipment operation constraints.The constraints of electric power balance are as follows:(7)∑i=1nCHPPit+∑i=1ndistriPdistrit+PGridt+Pstort=PLoadt+Pcondt,wherePGridt is the hourly power exchange value of the integrated energy system with public grid; PLoadt is the hourly electrical load; PStort is the charging or discharging power of battery; Pcondt is the hourly operation power of air conditioner; Pdistrit is the hourly operation power of distributed generation unit; and Pit is the power generated by the ith microturbine.The constraint of heating power supply is as follows:(8)∑i=1nCHPHit+∑i=1nboilerHboilerit+COPcond×Pcondt+Houtt−Hint≥Hspacet+HWatert.The constraint of hot water supply is as follows:(9)∑i=1nCHPHit+∑i=1nboilerHboilerit+Houtt−Hint≥HWatert,whereHit is the retrieved heat value by the ith microturbine through waste heat boiler; Hboilert is the heat output of the ith gas boiler; Hint and Houtt are the input and output of the thermal energy storage equipment; COPcond is the coefficient of performance of air-conditioning equipment; and Hspacet and Hwatert are the hourly space heating load and hot water load of the integrated energy system, respectively.The constraint of cooling power supply is as follows:(10)∑i=1nCHPCit+∑i=1nboilerCchilt+EERcond×Pcondt+Coutt−Cint≥Cspacet+CRefrit.The constraint of refrigeration power supply is as follows:(11)∑i=1nCHPCit+∑i=1nboilerCchilt+Coutt−Cint≥CRefrit,whereCit is the refrigeration power of the ith microturbine manufactured by the absorption refrigeration unit; Cchilt is the refrigeration power generated by the electric chiller; Cint and Coutt are the power input and output of the cold storage equipment, respectively; EERcond is the cooling energy efficiency ratio of air-conditioning equipment; and CSpacet and CRefrit are the hourly space cooling load and refrigeration load of the integrated energy system, respectively.The equipment capacity and operating constraints of the integrated energy system are as follows:for microturbines,(12)Pimin≤Pit≤Pimax,i∈nCHP.For gas boilers,(13)0≤Hboilerit≤Hboilerimax,i∈nboiler.For air-conditioning equipment,(14)0≤Pcondt≤Pcondmax.For heat (cold) energy storage equipment,(15)0≤Hint≤Hinmax,0≤Houtt≤Houtmax,Sstormin≤Sstort≤Sstormax,whereHint and Houtt are the power input and output of the heat (cold) storage equipment at time t; Hinmax and Houtmax are the power input and output limits of the heat (cold) storage equipment; and Sstort is the state of charge of the heat (cold) storage equipment at time t.(2)
During the operation of the integrated energy system, it is necessary to control the amount of pollutants emitted by the entire system within a certain range, including the total pollutant discharge constraints in the region and the hourly discharge constraints of various pollutants.Constraints on total carbon emissions in the region are as follows:(16)0≤PGC≤PGCmax.Constraints on total sulfide emissions in the region are as follows:(17)0≤PGS≤PGSmax.Constraints on total nitride emissions in the region are as follows:(18)0≤PGN≤PGNmax,wherePGCmax is the upper limit of the daily carbon emission value of the entire region; PGSmax is the upper limit of the daily sulfide emission value of the entire region; and PGNmax is the upper limit of the daily nitride emission value of the entire region.The hourly carbon emission constraints in the region are as follows:(19)0≤∑i=1nGiCO+TiCO+ZiCOt≤PGCt.The hourly sulfide emission constraints in the region are as follows:(20)0≤∑i=1nGiSO+TiSO+ZiSOt≤PGSt.The hourly nitride emission constraints in the region are as follows:(21)0≤∑i=1nGiNO+TiNO+ZiNOt≤PGNt,wherePGCt,PGSt, and PGNt are the upper limits of the carbon, sulfide, and nitride emission value of the entire region at time t.(3)
During the operation of the integrated energy system, the operating costs of the entire system needs to be controlled within a certain range. The economic constraint is as follows:(22)0≤KTotal≤Kmax,whereKTotal is the entire daily operating cost of the integrated energy system and Kmax is the upper limit of the daily operating cost.


## 3.3. Solving Algorithm

Due to the existence of a variety of energy supply equipment in the optimization model and the complex coupling relationship between cold, heat, and electrical energy, the model has a large number of optimization variables and a large volume of constraint conditions, when using traditional global solution space search algorithms such as NSGA-II or multiobjective particle swarm algorithm will make the solution process unacceptably lengthy which cannot meet the time limit requirements of the rolling optimization process within the dispatching system. When the traditional multiobjective optimization algorithm which transitions from single-objective optimization algorithm such as constraint method or weighting method is used to solve the model, the solution obtained is often limited to the local Pareto optimal solution set and it is difficult to obtain a large range of Pareto optimal solution set;, therefore, the traditional multiobjective optimization algorithm cannot provide the dispatching system with a variety of dispatching schemes for the system to make decisions.Based on the abovementioned analysis, this paper proposes a hierarchical progressive parallel NSGA-II algorithm suitable for solving multiobjective dispatching strategies of integrated energy systems. When solving the optimization model, it can not only make full use of the advantages of the global solution space search algorithm in the balance and diversity of the solution set searched but also achieve fast solution speed for large-scale nonlinear optimization problems which can meet the calculation time limit requirements of the rolling optimization process within the dispatching system. The algorithm has made three improvements based on NSGA-II as follows:(1)
Improve the nondominated layer sorting algorithm to achieve fast layering based on model characteristics.In the process of solving the joint multiobjective optimization operation model of the integrated energy system, considering that the model only includes dual optimization goals, this paper adopts a more efficient and fast nondominated layer sorting method according to the characteristics of the model to meet the requirements of rapid model solving.Firstly, the merge-sorting method is used to sort the individuals in the population by ascending order according to first objective function value. The merge sorting method is a sorting method by using the divide and conquer strategy. In divided stages, the method divides the problem into small questions and solves it recursively, while in the conquer stages, the method combines the answers obtained in divided stages together to get the results [28–30]. Suppose there is a population as shown in Figure 3 and each individual has two objective function values. The overall process of sorting the individuals in the population by ascending order according to first objective function value is shown in Figure 4. In the divided process, this structure is much like a complete binary tree and it can be understood as the process of recursively splitting subsequences. In the conquer process, we need to merge two already ordered subsequences into an ordered sequence, such as the last merge in Figure 4, to merge the already ordered subsequences of [4, 5, 7, 8] and [1–3, 6] into the final sequence [1–8]. The last merge within the conquer process demonstration is shown in Figure 5.If the first objective function values of the two individuals are the same, the individuals are sorted in ascending order according to the second objective function value and the result of the merge sorting is stored in the arrayFi. In this case, the individual with the smallest first objective function value must be an individual in the first layer of the nondominated layer. Then, the adjacent individual of the first individual is searched; if the second objective function value of the adjacent individual is smaller than the previous individual, the individual and the previous individual are belonging to the same layer of the nondominated layer; otherwise, the individual is not belonging to the same layer and will be inserted at the end of the array Fi and so on to find all first-level nondominated individuals in the population. Next, start with the first element at the end of the array Fi to find the second layer of nondominated individuals and so on until the population is completely stratified. In the example described above, the process of stratification occurs which is shown in Figure 6.Since the number of individuals in the next generation is half of that inFi, an incomplete layering strategy is adopted during the layering process to interrupt layering as soon as the number of individuals in the next generation exceeds the next generation population size, which can avoid many unnecessary layering comparison steps. The schematic diagram of the nondominated layer layering algorithm based on improved NSGA-II is shown in Figure 7:In the nondominated layer sorting algorithm, the average time complexity of the merge sorting process is O(N ∗ log2N) and the total times of calculations spent in the layering operation are(23)TTotal=∑i=1kTi=∑i=1kN−i=N−1+N−2+⋯+N−k<kN,whereN is twice the size of the population size; k is the total number of layers after the nondominated sorting is completed; TTotal is the total times of calculations spent in the layering operation after the nondominated sorting is completed; and Ti is the maximum times of comparisons required for generating the ith nondominated layer.The time complexity of the hierarchical (layering) operation is O (kN) after the nondominated sorting is completed and the total time complexity of the algorithm is O (N ∗ log2 N + kN). In general, k≪N, so the total time complexity of the algorithm is O (N ∗ log2 N). The sorting principle of NSGA-II is similar to the bubble sorting which needs to compare the individuals within the population in pairs and the average time complexity of the sorting algorithm is O (mN2), where m is the number of objectives. Compared with NSGA-II, the proposed algorithm reduces the time complexity from square order to linear logarithmic order.(2)
Modify the convergence judgment conditions and combine the target reaching method to accelerate the convergence process.The improved multiobjective genetic algorithm based on NSGA-II can quickly approach the vicinity of the Pareto front set in the process of solving the model. However, in order to achieve the final convergence, it requires a large number of iterations and function comparisons which means longer calculation time. Therefore, a comprehensive hierarchical progressive algorithm is proposed in this paper to speedup the process of finding the Pareto optimal solution set.Firstly, the aforementioned improved genetic algorithm based on NSGA-II is used for searching solutions close to the Pareto front by setting its termination condition to limited generations. Secondly, the target reaching method is utilized for a secondary search which uses the calculation result of the previous step as the initial value and a local search is performed by the target reaching method to obtain the Pareto solution set.The standard form of utilizing the target reaching method to solve the model is(24)minγs.t.fix−wi×γ≤Fi∗,cx≤0,ceqx=0,Ax≤b,Aeqx≤beq,lb≤x≤ub,where optimization variablex is an n-dimensional variable which includes hourly cold, heat, and power output of microturbine-based CHP system, hourly cold, heat, and power output of auxiliary energy supply equipment, hourly cold, heat, and power output of various energy storage equipment, SOC status of energy storage equipment, and hourly interacting power with the grid. fi(x) is the ith optimization objective function; ω is the weight coefficient of the ith objective function; Fi∗is the optimization target value expected to be achieved by the ith target; γ is the target achievement coefficient which is used for representing the achievement of the predetermined goal; c(x) and ceq(x) are the various nonlinear inequality constraints and nonlinear equality constraints met by the optimization variables; A and b are the coefficient matrices of various linear inequality constraints met by the optimization variables; Aeq and beq are the coefficient matrices of various linear equality constraints met by the optimization variables; lb and ub are the lower and upper limits of the optimization variable x, respectively.In the target reaching method, the multiobjective optimization problem is transformed into a single-objective optimization problem by using virtual variableγ and the single-objective optimization model is solved by using the interior point method. In the solving process, the interior point method constructs a Lagrangian auxiliary function which meets the Karush–Kuhn–Tucker conditional:(25)Lx,λ=fx+∑λc,icix+∑λA,iAx−b+∑λceq,iceqix+∑λAeq,iAeqx−beq.The second-order derivative Hessian matrix in its correction equation is(26)H=∇xx2Lx,λ=∇2fx+∑λi∇2cix+∑λi∇2ceqix,whereλc,iλA,iλceq,iλAeq,i are Lagrange multipliers which possess the same dimensions with respective constraints. x, lb, and ub are n × 1 column vectors; b is a p × 1 column vector; beq is a q × 1 column vector; A is a p × n matrix; Aeq is a q × n matrix; Hessian matrix is an n × n matrix.Solutions can be achieved by transformation, LDLT decomposition, iteration, and correction of the matrices in the correction equations. The iteration is automatically stopped when the convergence conditions are met and the solution closest to the expected optimization target in the feasible region is found as a Pareto optimal solution.(3)
Introduce parallel computing technology to accelerate the solution process according to the characteristics of the algorithm.Because there are a large number of repetitive loop iterative processes in the multiobjective hierarchical progressive algorithm based on improved NSGA-II, in addition, each process has no parameter passing to each other and does not depend on the results of other process, parallel computing technology is introduced to accelerate the convergence process of the algorithm. Among them, the part applicable for parallel computing mainly includes the nondominated layer merge sorting calculation and the call of the target reaching method for achieving fast convergence. In the nondominated layer merge sorting section, the fitness of individuals is first calculated on multiple parallel processing units, and then the individuals of the population to be sorted are evenly distributed tom parallel processing units. The complexity of the merge sorting process is reduced to O (N/m ∗ log2N) under parallel processing and the algorithm speeds up with the increase of the parallel processing units. When the target reaching method is utilized for a secondary search, the method starts from different initial values and there is no parameter passing to each other, so it is also suitable for using parallel processing technology to evenly distribute the computing task to multiple parallel processing units of the computer.In addition, because the GPU which is designed for parallel operations has more computing units than the CPU, in order to make full use of the GPU in large-scale array matrix processing, multiple operations without complex logical judgments and repeated calculations with little data correlation, the population crossover and mutation operations in the algorithm are assigned to the GPU. As the GPU is limited by the hardware architecture, single-precision calculations are much faster than double-precision calculations; therefore, this paper uses single-precision (single type variable) calculations for various types of population operations.Through the improvement of the abovementioned three aspects, the hierarchical progressive parallel NSGA-II algorithm can greatly reduce the calculation time in the process of obtaining the Pareto solution set, which can meet the requirements of the rolling optimization process within the dispatching system while ensuring the uniform distribution of the solution and the diversity of the population. The solution set obtained by the algorithm can provide flexible and diverse dispatching schemes for decision makers to choose.The hierarchical progressive parallel NSGA-II algorithm is written in MATLAB to solve the joint multiobjective optimization dispatching model of the integrated energy system in this paper. The algorithm flowchart is shown in Figure8.Figure 3
The individuals with two objective function values in a population.Figure 4
The overall process of the merge-sorting method.Figure 5
The last merge within the conquer process demonstration.Figure 6
The process of stratification in the example.Figure 7
Schematic diagram of nondominated layer layering algorithm based on improved NSGA-IIFigure 8
Flowchart of hierarchical progressive parallel NSGA-II algorithm.


## 4. Case Scenario

### 4.1. Case Situation

Figure9 shows the daily forecasting cooling, heating, and power load curve of a typical integrated energy system in the Tianjin Sino-Singapore Eco-city smart grid zone. The region adopts time-of-use electricity pricing with peak time from 8 to 11 and 18 to 23, flat time from 7 to 8 and 11 to 18, and valley time from 23 to 7. The time-of-use electricity price is shown in Table 1.Figure 9
Curve of the forecasting cooling, heating, and power load of a typical integrated energy system in one day.Table 1
Time-of-use price.
PeriodElectricity price ($/kWh)Peak0.1159Valley0.0491Flat0.0816Assuming the natural gas price is 0.4642 $/cubic meter and the high calorific value is 8571 kcal/cubic meter, thus the gas price is converted to 0.0464 $/kW·h. The weight coefficientsδ1, δ2, and δ3 for carbon, sulfur, and nitrogen emissions are taken as 1,100,100, respectively.The unit of emissions is converted to g/kWh according to [31] and the pollutant emission data of each equipment in the integrated energy system is shown in Table 2.Table 2
Pollutant indicators of major gas and fuel equipment in the integrated energy system.
PollutantMicroturbineInternal combustion engineFuel cellGas boilerPublic gridAbsorption refrigeration unitCOX724543486254922171NOX0.20.150.01360.542.2950.06SOX0.00360.00320.33270.7643.5830.0008Note. The pollutant indicators of microturbine, internal combustion engine, fuel cell, and public grid are based on the production of 1 kWh electricity; the pollutant indicators of gas boilers are based on the production of 1 kWh heat; the pollutant indicators of the absorption refrigeration unit are based on the production of 1 kWh cooling energy.


### 4.2. Equipment Parameters in the Case

The main equipment in the integrated energy system contains microturbines, waste heat boilers, absorption refrigerators, electric refrigerators, gas boilers, batteries, heat (cold) storage devices, household air conditioners, and photovoltaic cells. The integrated energy system exchanges power through a centralized power bus with the public grid, purchases power from the grid when the power supply within the system is insufficient, and sells power to the grid when the power supply within the system is surplus. At the same time, it is assumed that there is no gas production in the integrated energy system and all the gas required is supplied by the external gas company. The parameters of the energy supply equipment and energy storage equipment in the case are shown in Tables3 and 4.Table 3
The parameters of energy supply equipment.
EquipmentParametersValuesCapstone C1000 microturbine systemMaximum power generatedPc1000,max1000 kWRated efficiencyηc10000.33Waste heat boilerMaximum power inputPEB,max2000 kWRated efficiencyηEB0.8Lithium bromide absorption refrigeration unitMaximum power inputPAC,max2000 kWCoefficient of performanceCOPAC1.2Gas boilerMaximum power inputPGB,max1000 kWRated efficiencyηGB0.9Electric refrigeratorMaximum power inputPchil,max500 kWCoefficient of performanceCOPchil4.0Household air conditionerMaximum power inputPcond,max1000 kWEnergy efficiency ratioEERcond2.6Coefficient of performanceCOPcond3.1Public gridMaximum power exchangePBus,max1500 kWPhotovoltaic cellsMaximum power generatedPDG,max187.8 kWTable 4
The parameters of energy storage equipment.
ParametersStorageLead-acid batteriesHeat storage tankCold water storage tankCharging/(heat, cold) efficiency0.970.950.95Discharge/(heat, cold) efficiency0.970.950.95Maximum charge/(heat, cold) rate0.20.20.2Maximum discharge/(heat, cold) rate0.30.20.2Self-discharge/(heat, cold) rate0.020.030.03Maximum state of charge/energy0.90.90.9Minimum state of charge/energy0.20.10.1Capacity200 kW·h1000 kW·h500 kW·hThe calculation results of the case are analyzed from three aspects.


## 4.1. Case Situation

Figure9 shows the daily forecasting cooling, heating, and power load curve of a typical integrated energy system in the Tianjin Sino-Singapore Eco-city smart grid zone. The region adopts time-of-use electricity pricing with peak time from 8 to 11 and 18 to 23, flat time from 7 to 8 and 11 to 18, and valley time from 23 to 7. The time-of-use electricity price is shown in Table 1.Figure 9
Curve of the forecasting cooling, heating, and power load of a typical integrated energy system in one day.Table 1
Time-of-use price.
PeriodElectricity price ($/kWh)Peak0.1159Valley0.0491Flat0.0816Assuming the natural gas price is 0.4642 $/cubic meter and the high calorific value is 8571 kcal/cubic meter, thus the gas price is converted to 0.0464 $/kW·h. The weight coefficientsδ1, δ2, and δ3 for carbon, sulfur, and nitrogen emissions are taken as 1,100,100, respectively.The unit of emissions is converted to g/kWh according to [31] and the pollutant emission data of each equipment in the integrated energy system is shown in Table 2.Table 2
Pollutant indicators of major gas and fuel equipment in the integrated energy system.
PollutantMicroturbineInternal combustion engineFuel cellGas boilerPublic gridAbsorption refrigeration unitCOX724543486254922171NOX0.20.150.01360.542.2950.06SOX0.00360.00320.33270.7643.5830.0008Note. The pollutant indicators of microturbine, internal combustion engine, fuel cell, and public grid are based on the production of 1 kWh electricity; the pollutant indicators of gas boilers are based on the production of 1 kWh heat; the pollutant indicators of the absorption refrigeration unit are based on the production of 1 kWh cooling energy.


## 4.2. Equipment Parameters in the Case

The main equipment in the integrated energy system contains microturbines, waste heat boilers, absorption refrigerators, electric refrigerators, gas boilers, batteries, heat (cold) storage devices, household air conditioners, and photovoltaic cells. The integrated energy system exchanges power through a centralized power bus with the public grid, purchases power from the grid when the power supply within the system is insufficient, and sells power to the grid when the power supply within the system is surplus. At the same time, it is assumed that there is no gas production in the integrated energy system and all the gas required is supplied by the external gas company. The parameters of the energy supply equipment and energy storage equipment in the case are shown in Tables3 and 4.Table 3
The parameters of energy supply equipment.
EquipmentParametersValuesCapstone C1000 microturbine systemMaximum power generatedPc1000,max1000 kWRated efficiencyηc10000.33Waste heat boilerMaximum power inputPEB,max2000 kWRated efficiencyηEB0.8Lithium bromide absorption refrigeration unitMaximum power inputPAC,max2000 kWCoefficient of performanceCOPAC1.2Gas boilerMaximum power inputPGB,max1000 kWRated efficiencyηGB0.9Electric refrigeratorMaximum power inputPchil,max500 kWCoefficient of performanceCOPchil4.0Household air conditionerMaximum power inputPcond,max1000 kWEnergy efficiency ratioEERcond2.6Coefficient of performanceCOPcond3.1Public gridMaximum power exchangePBus,max1500 kWPhotovoltaic cellsMaximum power generatedPDG,max187.8 kWTable 4
The parameters of energy storage equipment.
ParametersStorageLead-acid batteriesHeat storage tankCold water storage tankCharging/(heat, cold) efficiency0.970.950.95Discharge/(heat, cold) efficiency0.970.950.95Maximum charge/(heat, cold) rate0.20.20.2Maximum discharge/(heat, cold) rate0.30.20.2Self-discharge/(heat, cold) rate0.020.030.03Maximum state of charge/energy0.90.90.9Minimum state of charge/energy0.20.10.1Capacity200 kW·h1000 kW·h500 kW·hThe calculation results of the case are analyzed from three aspects.


## 5. Case Study

### 5.1. Comparisons of Pareto Boundaries for System Optimal Operation under Different Algorithms

Three different algorithms were used to solve the joint multiobjective optimal operation and dispatching model of the integrated energy system, and the Pareto optimal solution sets under different algorithms were obtained.When using the weighting method, an initial solution is randomly generated within the feasible solution range as the initial iterative value. The weights of the two targets are normalized and 49 sets of weights are taken for iteration. The Pareto optimal solution set boundary obtained is shown in Figure10. As shown in the figure, the Pareto optimal solution set boundary obtained by the weighting method is relatively concentrated in some areas and it is not very continuous on the entire solution set front; in other words, it has disadvantages in the diversity of the solutions searched.Figure 10
Pareto optimal solution boundary obtained by the weighting method.When the NSGA-II algorithm is used, the initial population size is set to 2000 and single-precision real-number coding is used. Each bit of the individual is the electric, heating, and cooling power output values by various types of energy supply equipment. The crossover operator uses two-point crossover and adjacent individuals cross. The crossover rate is 0.8 and the mutation probability is 0.04. The maximum iteration generation is set to 200 and the relative change criterion between the two generations of individuals in the Pareto optimal solution set is set to 1e-1. If the difference between the two generations is less than the criterion, the iteration is terminated. The Pareto optimal solution set boundary obtained is shown in Figure 11.Figure 11
Pareto optimal solution boundary obtained by NSGA-II.As shown in the figure, the Pareto optimal solution set boundary obtained by the NSGA-II algorithm is relatively scattered and continuous on the entire front of the solution set, which has a prominent advantage in the diversity of the solutions searched, that is, the NSGA-II algorithm can find as many Pareto optimal solutions as possible. However, the disadvantage of NSGA-II algorithm is that the calculation time is very long and the optimization process can take several hours at a time, which is unacceptable in practical engineering applications.When the multiobjective hierarchical progressive parallel algorithm based on improved NSGA-II is used, the initial population size is set to 2000 and single-precision real-number coding is used. Each bit of the individual is the electric, heating, and cooling power output values by various types of energy supply equipment. The crossover operator uses two-point crossover and adjacent individuals cross. The crossover rate is 0.8, the mutation probability is 0.04, and the maximum iteration generation is set to 200. If the maximum iteration generation is exceeded, the target reaching method is utilized for a secondary search. The relative change criterion between the two generations of individuals in the Pareto optimal solution set is set to 1e-1. If the difference between the two generations is less than the criterion, the iteration is terminated. The test computer is equipped with Intel’s third-generation Core processor i5-3230M, which can perform dual-core four-thread parallel computing. The main frequency of each core of the CPU is 2.60 GHz; the configured GPU is NVIDIA GeForce 610M. The programming environment is MATLAB R2012a and the result is shown in Figure 12.Figure 12
Pareto optimal solution boundary obtained by hierarchical progressive parallel NSGA-II.As shown in the figure, the Pareto optimal solution set boundary obtained by the multiobjective hierarchical progressive parallel algorithm based on improved NSGA-II is also relatively scattered, and it is passable continuous on the entire solution set front. Because the target reaching method does not pay attention to the diversity of solutions, the diversity of solutions will be sacrificed in the secondary search process in exchange for the convergence time of the calculation. The advantages of the overall algorithm in the diversity of solutions are not as prominent as the NSGA-II algorithm. However, its advantage lies in the combination of the rapid convergence of the target reaching method and the advantages of the NSGA-II algorithm in searching for the diversity of solutions, which makes the algorithm to find as many and more dispersed Pareto optimal solutions as possible in a short time. The optimization process can be reduced to several minutes on the micro-PC, which greatly improves the engineering practicability of the algorithm.


### 5.2. Performance Comparison of Different Algorithms

NSGA-II and the multiobjective hierarchical progressive parallel algorithm based on improved NSGA-II proposed in this paper are used to solve the optimization operation model of the integrated energy system. The average individual comparison times in a single layering process within one solution process, nondominated frontier average number of layers, and the average CPU time consumed by a single layering process of the two algorithms are recorded. The parameter settings of the test are the same as the parameter conditions for obtaining the Pareto boundary using NSGA-II. During the test, the experiment was repeated 20 times for each algorithm and the average value of each test data was taken. The results of the time-consuming comparison tests between NSGA-II and the proposed algorithm on a single nondominated layer construction are shown in Table5.Table 5
Comparison of test results between NSGA-II and the proposed algorithm in every construction of a nondominated frontier set.
Number of testsNSGA-IIThe multiobjective hierarchical progressive parallel algorithmThe average individual comparison timesNondominated frontier average number of layersAverage CPU time consumed by a single layering process(s)The average individual comparison timesNondominated frontier average number of layersAverage CPU time consumed by a single layering process(s)14861898.458.593.362521.929.21.12624665847.750.689.673679.125.31.334585557.947.48861776.723.71.34844448335.541.985.477458.3411.94355255582.674.2100.975764.127.11.37964999614.9649670900.7281.69474983007.963.395.769879.441.71.65585375101.479103.274943.329.51.84895143954.269.898.876288.224.91.899105286829.575.5101.565054.227.71.472115387153.879.5103.445385.829.70.725124806556.656.392.365837.328.11.502134718012.152.790.659488.326.41.261144966130.462.695.448214.931.30.832154987274.163.595.848977.831.70.711164883584.659.393.844668.729.70.69717505473666.297.140401.223.10.535184854897.658.293.26545329.11.487194643834.149.889.270914.424.91.695204998534.463.99650881.1220.933Average value4945322.261.89562424.428.71.302From the table, the algorithm proposed in this paper has very obvious advantages in terms of the average individual comparison times in a single evolutionary layering process, the average number of nondominated frontier layering layers, and the average CPU time consumed by a single layering process. The decrease in the average individual comparison times is mainly due to the use of a nondominated frontier construction method based on dual-objective merge sorting. The decrease in average number of layers of the nondominated frontier is mainly due to the practice of stopping the stratification when the individual is full. The significant reduction in the average CPU time consumed by the layering process is mainly due to the reduction in the number of comparisons and the introduction of CPU dual-core parallel computing.The changes in the number of stratifications in the nondominated layer of the subparent comprehensive population during a typical evolutionary process are as shown in Figure13.Figure 13
Changes in the number of stratifications in the nondominated layer of the subparent comprehensive population in the evolutionary process.As shown in Figure13, the multiobjective hierarchical progressive parallel algorithm based on improved NSGA-II greatly reduces the number of layers in a single evolution process and significantly improved computing efficiency in the process of performing hierarchical selection operations compared with NSGA-II.Although several improvements have been made for the NSGA-II algorithm, due to the large number of evolutionary generations, the large number of comparisons, and the large amount of calculations, the total calculation time is still quite long and the algorithm execution efficiency is lower than the weighting method. Therefore, this paper applies the comprehensive hierarchical progressive algorithm which modifies the convergence condition judgment to solve the model. Taking advantage of the strong scalability of NSGA-II, by setting a certain truncated evolutionary generation and combining with the target reaching method to perform the final step of convergence calculation, the convergence progress of the entire algorithm is accelerated. The specific test calculation time is shown in Table6.Table 6
Time-consuming comparison in whole solution period of algorithms under different convergence conditions.
Number of testsTraditional multiobjective weighting algorithmMultiobjective hierarchical progressive parallel algorithm based on improved NSGA-IIAverage searching time for a single solution(s)Total calculation time(s)Average time taken per iteration(s)Truncated evolutionary generation, 100Truncated evolutionary generation, 200Number of Pareto solutions foundTotal calculation time (s)Number of Pareto solutions foundTotal calculation time(s)16.2306.12.843434.645716.626.23042.943445.946740.136293.23.242467.449788.945.3261.32.941416.146702.555.5268.92.941422.346710.666.5317.73.144470.748782.876.8334.52.945460.646754.686.8331.3344469.247772.995.7279.62.942425.546711.2106.4311.53.243474.749794.5117342.63.245486.949802.6125.5270.73.241454.249773.3136294.73.343475.753804.1146.1299.73.243468.949788.4155.9291.22.942433.446721.4165.7277.63.142449.648761.2176.8333.13.144475.448784.3186.8331.53.244484.449803.3196.7328.7344461.647759206.9336.83.245483.649799.3Average6.2305.73.143.145847.8763.6As shown in the table, when the truncated evolutionary generation is set to 100, the total calculation time of the entire algorithm can be controlled within 10 minutes, which is much faster than NSGA-II (up to several hours, not listed because the magnitudes are not of average significance). Due to the need for specifying the weight and initial position, the fastest traditional algorithm has the characteristics of fast single optimization process but poor overall convergence. It may take many trials to find a satisfactory Pareto front; overall, it is not as good as the algorithm proposed in this paper. The hierarchical progressive algorithm based on the improved convergence judgment accelerates the process of convergence by combining NSGA-II with the target reaching method, but at the expense of a small portion of the Pareto solution set diversity, which is not only reflected in the number of searched solutions but also in the reduction of the dispersion degree of searched solutions. The diversity is slightly weaker than NSGA-II, but since the initial value used by the target reaching method in the convergence calculation process was obtained from the improved NSGA-II algorithm after iteration for several generations and therefore these initial values still have relatively good dispersion, hence the solution diversity of the proposed algorithm is much better than that of traditional weighting methods. In general, the algorithm proposed in this paper has better comprehensive computing capabilities and advantages in terms of computing time and the diversity of solutions.In the multiobjective hierarchical progressive algorithm based on improved NSGA-II, parallel computing acceleration mainly includes two parts: CPU dual-core parallel acceleration and GPU + CPU cooperative processing acceleration. The CPU dual-core parallel acceleration part mainly includes the nondominated layer merge sorting calculation and the final step of calling the target reaching method to achieve rapid convergence. According to the description of introducing the parallel calculation to accelerate the solution process, it can be seen that the ideal speedup ratio is 2 considering the calculation amount is evenly divided into the two cores of the CPU. However, because the calculation amount that needs to be decomposed is not necessarily an even number and data splitting and integration, the speedup ratio is lower than 2 in actual tests. The specific average test results are shown in Table7.Table 7
Algorithmic acceleration test results using parallel computing technology.
Parallel acceleration position in the algorithmIdeal speedup ratioAverage test resultsCPU dual-core parallel accelerationNondominated layer merge sorting section21.834Calling target reaching method to achieve rapid convergence section21.668GPU + CPU cooperative processing accelerationPopulation crossover and mutation operations—1.082Since the GPU has more arithmetic units (multiplication and addition units of integer and floating point number, special arithmetic units, etc.) than the CPU and its design is more oriented to parallel operations, the GPU has a large advantage over the CPU in large-scale array and matrix processing, multiple simple calculations, and repeated operations with little data correlation in theory. At the specific execution level, GPU + CPU cooperative processing acceleration has advantages in terms of population crossover and mutation operation time in a single iteration. However, the elapsed time of GPU-CPU data transmission and communication needs to be added to the total calculation time; thus, in the testing system, the acceleration ratio can only reach an average of about 1.08 in the cooperative processing acceleration process.


### 5.3. Dispatching Strategy of Typical Optimal Solution

In the Pareto optimal solution set, the sum of the pollutant emission values of Pareto solution with the best economic index is 44571.567 kg/day and the operating cost is 3433.16 USD/day. The sum of the pollutant emission values of Pareto solution with the best environmental protection index is 34015.972 kg/day and the operating cost is 4344.79 USD/day. Select a typical Pareto solution that takes into account both economic and environmental protection from the Pareto solutions set; the sum of the pollutant emission values is 42771.233 kg/day and the operating cost is 3667.74 USD/day. The dispatching curve of various load balance and equipment operating status are listed below; the day-ahead optimized dispatching curve of electric load balance is shown in Figure14.Figure 14
Optimized dispatching curve of electric load balance with both economic and environmental advantages.In Figure14, the system purchases electricity from the public grid to meet the demand for electricity loads during the time when the electricity price is low from 23 : 00 to 8 : 00 and the purchase amount is less than the case of optimal economy index and greater than the case of optimal environmental protection. When the peak and flat electricity price periods come, the microturbine enters the rated operating state to meet the demand of the electric load and the total power generation amount is between the economically optimal situation and the environmentally friendly optimal situation.The day-ahead optimized dispatching curve of space heat load balance with both economic and environmental advantages is shown in Figure15.Figure 15
Optimized dispatching curve of space heat load balance with both economic and environmental advantages.In Figure15, the space heat load is jointly supplied by the air-conditioning system, gas boiler, and microgas turbine system. At the time of valley electricity prices, the space heat load is mainly provided by the air conditioner and during the peak electricity price periods, the waste heat boiler in the microgas turbine system and gas boiler are applied for heating supply.The day-ahead optimized dispatching curve of hot water load balance with both economic and environmental advantages is shown in Figure16.Figure 16
Optimized dispatching curve of hot water load balance with both economic and environmental advantages.In Figure16, the hot water load of the integrated energy system in the smart grid zone is met by waste heat boiler and heat storage equipment. The supply is similar to the dispatching strategy with optimal environmental protection indicators and the difference is reflected in the scheduling of heat energy storage.The day-ahead optimized dispatching curve of space cooling load balance with both economic and environmental advantages is shown in Figure17.Figure 17
Optimized dispatching curve of space cooling load balance with both economic and environmental advantages.In Figure17, the space cooling load is jointly supplied by the air-conditioning system and lithium bromide absorption refrigeration unit within the microturbine system. At the time of valley electricity prices, the space cooling load is mainly provided by the air conditioner and the remaining time is mainly provided by the lithium bromide absorption refrigeration unit. When the cooling capacity of the lithium bromide absorption refrigerating unit is insufficient, the household air-conditioning system is used as the peak shaving device of the space cooling load.The day-ahead optimized dispatching curve of refrigeration load balance with both economic and environmental advantages is shown in Figure18.Figure 18
Optimized dispatching curve of refrigeration load balance with both economic and environmental advantages.As shown in Figure18, the refrigeration load of the integrated energy system is jointly met by the lithium bromide absorption refrigeration unit, electric refrigerator, and cold storage equipment, and the output of the three types of energy supply equipment is relatively uniform.The day-ahead optimized dispatching curve of energy storage equipment operating status with both economic and environmental advantages is shown in Figure19.Figure 19
Optimized dispatching curve of energy storage equipment operating status with both economic and environmental advantages.The operating status of the energy storage equipment with both economic and environmental advantages is less than the average state of energy under the dispatching strategies with optimal economic indicators and environmental protection indicators; this is because it is not necessary to reach the limits of economic or environmental protection goals and energy storage equipment can operate with some margin when both goals are considered.The pollutant emissions from various types of energy supply equipment and public grid in the smart grid zone with both economic and environmental advantages are shown in Figure20.Figure 20
Pollutant emissions from various types of energy supply equipment in the optimized dispatching strategy with both economic and environmental advantages.As shown in Figure20, pollutant emissions mainly come from energy production emissions of public grids, followed by emissions from microturbine, CHP auxiliary systems, and gas boiler in the optimized dispatching strategy with both economic and environmental advantages. The strategy utilizes cold and heat energy storage equipment to store the excess energy generated by the microturbine and release it at the required time. When the pollutants emissions are not unacceptably high, the public grid electricity is used for energy supply during nonpeak pricing hours. The emission of pollutants in the scheme is between the best economic situation and the best environmental protection situation. It does not cause a large degree of pollutant emissions and controls the operating cost of the integrated energy system within a certain range, which is rather a balanced dispatching strategy. In practical engineering applications, the optimal operation scheme can be selected according to the needs of customers such as economic-oriented or environmental protection-oriented, thus achieving the multiobjective optimized operation of the integrated energy system in the smart grid zone.


## 5.1. Comparisons of Pareto Boundaries for System Optimal Operation under Different Algorithms

Three different algorithms were used to solve the joint multiobjective optimal operation and dispatching model of the integrated energy system, and the Pareto optimal solution sets under different algorithms were obtained.When using the weighting method, an initial solution is randomly generated within the feasible solution range as the initial iterative value. The weights of the two targets are normalized and 49 sets of weights are taken for iteration. The Pareto optimal solution set boundary obtained is shown in Figure10. As shown in the figure, the Pareto optimal solution set boundary obtained by the weighting method is relatively concentrated in some areas and it is not very continuous on the entire solution set front; in other words, it has disadvantages in the diversity of the solutions searched.Figure 10
Pareto optimal solution boundary obtained by the weighting method.When the NSGA-II algorithm is used, the initial population size is set to 2000 and single-precision real-number coding is used. Each bit of the individual is the electric, heating, and cooling power output values by various types of energy supply equipment. The crossover operator uses two-point crossover and adjacent individuals cross. The crossover rate is 0.8 and the mutation probability is 0.04. The maximum iteration generation is set to 200 and the relative change criterion between the two generations of individuals in the Pareto optimal solution set is set to 1e-1. If the difference between the two generations is less than the criterion, the iteration is terminated. The Pareto optimal solution set boundary obtained is shown in Figure 11.Figure 11
Pareto optimal solution boundary obtained by NSGA-II.As shown in the figure, the Pareto optimal solution set boundary obtained by the NSGA-II algorithm is relatively scattered and continuous on the entire front of the solution set, which has a prominent advantage in the diversity of the solutions searched, that is, the NSGA-II algorithm can find as many Pareto optimal solutions as possible. However, the disadvantage of NSGA-II algorithm is that the calculation time is very long and the optimization process can take several hours at a time, which is unacceptable in practical engineering applications.When the multiobjective hierarchical progressive parallel algorithm based on improved NSGA-II is used, the initial population size is set to 2000 and single-precision real-number coding is used. Each bit of the individual is the electric, heating, and cooling power output values by various types of energy supply equipment. The crossover operator uses two-point crossover and adjacent individuals cross. The crossover rate is 0.8, the mutation probability is 0.04, and the maximum iteration generation is set to 200. If the maximum iteration generation is exceeded, the target reaching method is utilized for a secondary search. The relative change criterion between the two generations of individuals in the Pareto optimal solution set is set to 1e-1. If the difference between the two generations is less than the criterion, the iteration is terminated. The test computer is equipped with Intel’s third-generation Core processor i5-3230M, which can perform dual-core four-thread parallel computing. The main frequency of each core of the CPU is 2.60 GHz; the configured GPU is NVIDIA GeForce 610M. The programming environment is MATLAB R2012a and the result is shown in Figure 12.Figure 12
Pareto optimal solution boundary obtained by hierarchical progressive parallel NSGA-II.As shown in the figure, the Pareto optimal solution set boundary obtained by the multiobjective hierarchical progressive parallel algorithm based on improved NSGA-II is also relatively scattered, and it is passable continuous on the entire solution set front. Because the target reaching method does not pay attention to the diversity of solutions, the diversity of solutions will be sacrificed in the secondary search process in exchange for the convergence time of the calculation. The advantages of the overall algorithm in the diversity of solutions are not as prominent as the NSGA-II algorithm. However, its advantage lies in the combination of the rapid convergence of the target reaching method and the advantages of the NSGA-II algorithm in searching for the diversity of solutions, which makes the algorithm to find as many and more dispersed Pareto optimal solutions as possible in a short time. The optimization process can be reduced to several minutes on the micro-PC, which greatly improves the engineering practicability of the algorithm.


## 5.2. Performance Comparison of Different Algorithms

NSGA-II and the multiobjective hierarchical progressive parallel algorithm based on improved NSGA-II proposed in this paper are used to solve the optimization operation model of the integrated energy system. The average individual comparison times in a single layering process within one solution process, nondominated frontier average number of layers, and the average CPU time consumed by a single layering process of the two algorithms are recorded. The parameter settings of the test are the same as the parameter conditions for obtaining the Pareto boundary using NSGA-II. During the test, the experiment was repeated 20 times for each algorithm and the average value of each test data was taken. The results of the time-consuming comparison tests between NSGA-II and the proposed algorithm on a single nondominated layer construction are shown in Table5.Table 5
Comparison of test results between NSGA-II and the proposed algorithm in every construction of a nondominated frontier set.
Number of testsNSGA-IIThe multiobjective hierarchical progressive parallel algorithmThe average individual comparison timesNondominated frontier average number of layersAverage CPU time consumed by a single layering process(s)The average individual comparison timesNondominated frontier average number of layersAverage CPU time consumed by a single layering process(s)14861898.458.593.362521.929.21.12624665847.750.689.673679.125.31.334585557.947.48861776.723.71.34844448335.541.985.477458.3411.94355255582.674.2100.975764.127.11.37964999614.9649670900.7281.69474983007.963.395.769879.441.71.65585375101.479103.274943.329.51.84895143954.269.898.876288.224.91.899105286829.575.5101.565054.227.71.472115387153.879.5103.445385.829.70.725124806556.656.392.365837.328.11.502134718012.152.790.659488.326.41.261144966130.462.695.448214.931.30.832154987274.163.595.848977.831.70.711164883584.659.393.844668.729.70.69717505473666.297.140401.223.10.535184854897.658.293.26545329.11.487194643834.149.889.270914.424.91.695204998534.463.99650881.1220.933Average value4945322.261.89562424.428.71.302From the table, the algorithm proposed in this paper has very obvious advantages in terms of the average individual comparison times in a single evolutionary layering process, the average number of nondominated frontier layering layers, and the average CPU time consumed by a single layering process. The decrease in the average individual comparison times is mainly due to the use of a nondominated frontier construction method based on dual-objective merge sorting. The decrease in average number of layers of the nondominated frontier is mainly due to the practice of stopping the stratification when the individual is full. The significant reduction in the average CPU time consumed by the layering process is mainly due to the reduction in the number of comparisons and the introduction of CPU dual-core parallel computing.The changes in the number of stratifications in the nondominated layer of the subparent comprehensive population during a typical evolutionary process are as shown in Figure13.Figure 13
Changes in the number of stratifications in the nondominated layer of the subparent comprehensive population in the evolutionary process.As shown in Figure13, the multiobjective hierarchical progressive parallel algorithm based on improved NSGA-II greatly reduces the number of layers in a single evolution process and significantly improved computing efficiency in the process of performing hierarchical selection operations compared with NSGA-II.Although several improvements have been made for the NSGA-II algorithm, due to the large number of evolutionary generations, the large number of comparisons, and the large amount of calculations, the total calculation time is still quite long and the algorithm execution efficiency is lower than the weighting method. Therefore, this paper applies the comprehensive hierarchical progressive algorithm which modifies the convergence condition judgment to solve the model. Taking advantage of the strong scalability of NSGA-II, by setting a certain truncated evolutionary generation and combining with the target reaching method to perform the final step of convergence calculation, the convergence progress of the entire algorithm is accelerated. The specific test calculation time is shown in Table6.Table 6
Time-consuming comparison in whole solution period of algorithms under different convergence conditions.
Number of testsTraditional multiobjective weighting algorithmMultiobjective hierarchical progressive parallel algorithm based on improved NSGA-IIAverage searching time for a single solution(s)Total calculation time(s)Average time taken per iteration(s)Truncated evolutionary generation, 100Truncated evolutionary generation, 200Number of Pareto solutions foundTotal calculation time (s)Number of Pareto solutions foundTotal calculation time(s)16.2306.12.843434.645716.626.23042.943445.946740.136293.23.242467.449788.945.3261.32.941416.146702.555.5268.92.941422.346710.666.5317.73.144470.748782.876.8334.52.945460.646754.686.8331.3344469.247772.995.7279.62.942425.546711.2106.4311.53.243474.749794.5117342.63.245486.949802.6125.5270.73.241454.249773.3136294.73.343475.753804.1146.1299.73.243468.949788.4155.9291.22.942433.446721.4165.7277.63.142449.648761.2176.8333.13.144475.448784.3186.8331.53.244484.449803.3196.7328.7344461.647759206.9336.83.245483.649799.3Average6.2305.73.143.145847.8763.6As shown in the table, when the truncated evolutionary generation is set to 100, the total calculation time of the entire algorithm can be controlled within 10 minutes, which is much faster than NSGA-II (up to several hours, not listed because the magnitudes are not of average significance). Due to the need for specifying the weight and initial position, the fastest traditional algorithm has the characteristics of fast single optimization process but poor overall convergence. It may take many trials to find a satisfactory Pareto front; overall, it is not as good as the algorithm proposed in this paper. The hierarchical progressive algorithm based on the improved convergence judgment accelerates the process of convergence by combining NSGA-II with the target reaching method, but at the expense of a small portion of the Pareto solution set diversity, which is not only reflected in the number of searched solutions but also in the reduction of the dispersion degree of searched solutions. The diversity is slightly weaker than NSGA-II, but since the initial value used by the target reaching method in the convergence calculation process was obtained from the improved NSGA-II algorithm after iteration for several generations and therefore these initial values still have relatively good dispersion, hence the solution diversity of the proposed algorithm is much better than that of traditional weighting methods. In general, the algorithm proposed in this paper has better comprehensive computing capabilities and advantages in terms of computing time and the diversity of solutions.In the multiobjective hierarchical progressive algorithm based on improved NSGA-II, parallel computing acceleration mainly includes two parts: CPU dual-core parallel acceleration and GPU + CPU cooperative processing acceleration. The CPU dual-core parallel acceleration part mainly includes the nondominated layer merge sorting calculation and the final step of calling the target reaching method to achieve rapid convergence. According to the description of introducing the parallel calculation to accelerate the solution process, it can be seen that the ideal speedup ratio is 2 considering the calculation amount is evenly divided into the two cores of the CPU. However, because the calculation amount that needs to be decomposed is not necessarily an even number and data splitting and integration, the speedup ratio is lower than 2 in actual tests. The specific average test results are shown in Table7.Table 7
Algorithmic acceleration test results using parallel computing technology.
Parallel acceleration position in the algorithmIdeal speedup ratioAverage test resultsCPU dual-core parallel accelerationNondominated layer merge sorting section21.834Calling target reaching method to achieve rapid convergence section21.668GPU + CPU cooperative processing accelerationPopulation crossover and mutation operations—1.082Since the GPU has more arithmetic units (multiplication and addition units of integer and floating point number, special arithmetic units, etc.) than the CPU and its design is more oriented to parallel operations, the GPU has a large advantage over the CPU in large-scale array and matrix processing, multiple simple calculations, and repeated operations with little data correlation in theory. At the specific execution level, GPU + CPU cooperative processing acceleration has advantages in terms of population crossover and mutation operation time in a single iteration. However, the elapsed time of GPU-CPU data transmission and communication needs to be added to the total calculation time; thus, in the testing system, the acceleration ratio can only reach an average of about 1.08 in the cooperative processing acceleration process.


## 5.3. Dispatching Strategy of Typical Optimal Solution

In the Pareto optimal solution set, the sum of the pollutant emission values of Pareto solution with the best economic index is 44571.567 kg/day and the operating cost is 3433.16 USD/day. The sum of the pollutant emission values of Pareto solution with the best environmental protection index is 34015.972 kg/day and the operating cost is 4344.79 USD/day. Select a typical Pareto solution that takes into account both economic and environmental protection from the Pareto solutions set; the sum of the pollutant emission values is 42771.233 kg/day and the operating cost is 3667.74 USD/day. The dispatching curve of various load balance and equipment operating status are listed below; the day-ahead optimized dispatching curve of electric load balance is shown in Figure14.Figure 14
Optimized dispatching curve of electric load balance with both economic and environmental advantages.In Figure14, the system purchases electricity from the public grid to meet the demand for electricity loads during the time when the electricity price is low from 23 : 00 to 8 : 00 and the purchase amount is less than the case of optimal economy index and greater than the case of optimal environmental protection. When the peak and flat electricity price periods come, the microturbine enters the rated operating state to meet the demand of the electric load and the total power generation amount is between the economically optimal situation and the environmentally friendly optimal situation.The day-ahead optimized dispatching curve of space heat load balance with both economic and environmental advantages is shown in Figure15.Figure 15
Optimized dispatching curve of space heat load balance with both economic and environmental advantages.In Figure15, the space heat load is jointly supplied by the air-conditioning system, gas boiler, and microgas turbine system. At the time of valley electricity prices, the space heat load is mainly provided by the air conditioner and during the peak electricity price periods, the waste heat boiler in the microgas turbine system and gas boiler are applied for heating supply.The day-ahead optimized dispatching curve of hot water load balance with both economic and environmental advantages is shown in Figure16.Figure 16
Optimized dispatching curve of hot water load balance with both economic and environmental advantages.In Figure16, the hot water load of the integrated energy system in the smart grid zone is met by waste heat boiler and heat storage equipment. The supply is similar to the dispatching strategy with optimal environmental protection indicators and the difference is reflected in the scheduling of heat energy storage.The day-ahead optimized dispatching curve of space cooling load balance with both economic and environmental advantages is shown in Figure17.Figure 17
Optimized dispatching curve of space cooling load balance with both economic and environmental advantages.In Figure17, the space cooling load is jointly supplied by the air-conditioning system and lithium bromide absorption refrigeration unit within the microturbine system. At the time of valley electricity prices, the space cooling load is mainly provided by the air conditioner and the remaining time is mainly provided by the lithium bromide absorption refrigeration unit. When the cooling capacity of the lithium bromide absorption refrigerating unit is insufficient, the household air-conditioning system is used as the peak shaving device of the space cooling load.The day-ahead optimized dispatching curve of refrigeration load balance with both economic and environmental advantages is shown in Figure18.Figure 18
Optimized dispatching curve of refrigeration load balance with both economic and environmental advantages.As shown in Figure18, the refrigeration load of the integrated energy system is jointly met by the lithium bromide absorption refrigeration unit, electric refrigerator, and cold storage equipment, and the output of the three types of energy supply equipment is relatively uniform.The day-ahead optimized dispatching curve of energy storage equipment operating status with both economic and environmental advantages is shown in Figure19.Figure 19
Optimized dispatching curve of energy storage equipment operating status with both economic and environmental advantages.The operating status of the energy storage equipment with both economic and environmental advantages is less than the average state of energy under the dispatching strategies with optimal economic indicators and environmental protection indicators; this is because it is not necessary to reach the limits of economic or environmental protection goals and energy storage equipment can operate with some margin when both goals are considered.The pollutant emissions from various types of energy supply equipment and public grid in the smart grid zone with both economic and environmental advantages are shown in Figure20.Figure 20
Pollutant emissions from various types of energy supply equipment in the optimized dispatching strategy with both economic and environmental advantages.As shown in Figure20, pollutant emissions mainly come from energy production emissions of public grids, followed by emissions from microturbine, CHP auxiliary systems, and gas boiler in the optimized dispatching strategy with both economic and environmental advantages. The strategy utilizes cold and heat energy storage equipment to store the excess energy generated by the microturbine and release it at the required time. When the pollutants emissions are not unacceptably high, the public grid electricity is used for energy supply during nonpeak pricing hours. The emission of pollutants in the scheme is between the best economic situation and the best environmental protection situation. It does not cause a large degree of pollutant emissions and controls the operating cost of the integrated energy system within a certain range, which is rather a balanced dispatching strategy. In practical engineering applications, the optimal operation scheme can be selected according to the needs of customers such as economic-oriented or environmental protection-oriented, thus achieving the multiobjective optimized operation of the integrated energy system in the smart grid zone.


## 6. Conclusion

This paper analyzes the typical energy supply structure of the integrated energy system and the joint multiobjective optimization dispatching architecture, researches on the optimization dispatching model of the integrated energy system under multiple objectives, and proposes multiobjective hierarchical progressive parallel algorithm based on improved NSGA-II according to the characteristics of the model.The algorithm improves the nondominated layer sorting algorithm, changes the convergence judgment condition while introducing the target reaching method to accelerate convergence, and introduces parallel computing technology according to the characteristics of the algorithm. The algorithm accelerates the convergence speed of the solution process while taking into account the diversity of the solution searched.By analyzing the application scenarios of the optimized operation of the integrated energy system under multiobjectives, the typical daily operation optimization model of integrated energy systems in the smart grid zone is solved, and system’s joint multiobjective optimized operation Pareto boundary is obtained. Finally, the algorithm proposed in this paper is compared with traditional weighting algorithm and NSGA-II algorithm in terms of solution diversity and solution time. The case study shows that the algorithm proposed in this paper not only has advantages on the diversity of solutions searched but also can achieve better results in many aspects such as the iteration time and algorithm convergence which are required in practical engineering projects. In practical engineering applications, the optimal operation scheme can be selected according to the needs of customers to achieve multiobjective optimized operation of the integrated energy system in the smart grid zone.



---

*Source: 6541782-2020-05-31.xml*