A Hybrid Optimization Algorithm for Water Volume Adjustment Problem in District Heating Systems

Nowadays, the winter is getting harsher and harsher in Northern China. Thus, the centralized heating systems (CHSs) are playing even more irreplaceable, essential and critical roles in ensuring general public’s livelihood than never ever before. CHSs are normally composed of one or several combined heat and power (CHP) plants (units) and an extensive vein like district heating networks (DHNs) connecting with chemical plants, paper mills, food processing factories, hospitals, hotels, universities, prisons and residential complexes. A CHP plant in Northern China usually consumes coal to heat the cold water into steam to drive high-pressure turbines and low pressure turbines to generate electricity. Then the low-temperature steam is used to heat up the cold water in a main pipe into hot water travelling through the DHNs to provide heat to each end nodes. The returned water will be heated again for reuse and the surplus steam will be released into air through cooling towers. In 2020, China promised to the world that carbon dioxide will peak in 2030 and net-zero emission will happen in 2060. On the one hand, CHP plants need to guarantee enough hot water flowing within each household’s heating radiator. On the other hand, they should cut down on the consumption on non-renewable resources. Lowering water temperature, adjusting water volume and reducing water pressure will all contribute to energy-saving purpose. Lowering water temperature and reducing water pressure may cause too much heat losses during long-distance transmission in frigid winter. Therefore, a reasonable water volume adjustment becomes an advisable action comparatively. Here, we transfer the hot water supply volume optimization problem (HWSVOP) into a heat exchange station (HES) valve angle adjustment problem (SHWESVAAP). Then, a multi-objective mathematical model is established considering balancing the satisfactory degree of each household in residential quarters and the hot water volume (HWV) in the main pipe. And a hybrid polar bear optimization algorithm integrated with chemical reaction optimization (HA-PBO-CRO) is designed to optimize the valve angle (VA) in each HES. The comparative results between HA-PBO-CRO and non-dominant sorting genetic algorithm (NSGAII) demonstrate HA-PBO-CRO is superior to NSGAII with better Pareto frontiers on one hand and provide a critical reference supporting the management in a CHP plant to make a right decision on what to do to cut energy consumption while satisfying customers’ needs.


Introduction
Centralized heating system (CHS), which is composed of heat sources, heat exchange stations, heat users and district heating networks, is playing even crucial role in escorting northern people through unbearable and annoying cold weather [1]. In a CHS, there are high coupling relationships and great mutual influences among heat sources, networks and users. Normally, heat sources mean CHPs which operate in a combined heat and electricity generation mode producing hot water and electricity together by burning coals and other energy resources [2]. District heating networks mean the underground and built-in-buildings water pipes that reach to each radiator under the window of every northern family. Heat users mean those local residents who rely heavily on those CHPs to deliver hot water into their homes during winter time. The general operational modes are shown in Figs. 1, 2 [3,4].
The overheated hot water from a CHP plant which will go through an extraordinary wide main pipe cannot be directly pumped into each household. Traditionally, there are some HESs standing in between to neutralize the water temperature to an acceptable degree. And thereafter, the hot water in a HES is sent out to warm up each flat in a residential quarter. A simplified process showing how a HES works is shown in Fig. 3.
Right now, what are bothering those CHPs most in northern China is indeed how much hot water they should provide [5]. The reason is the more the hot water is, the more the electricity and energy consumption will be. With more and more greener energies like hydro, wind and solar powers enter into State GRID, more coalfire power means heavier penalties. In Fig. 3, the regular operations on regulating hot water volume are carried out by an experienced operator controlling the gate (disc) angle Ɵ of a butterfly valve (see Fig. 4). As a result, many complaints and grumbles from local people arrive in flocks since a manual adjustment fail to meet the majority's satisfactions. In view of this, if we could set a wise gate (disc) angle Ɵ for valve E, we could help CHPs to make the right decision about how much hot water they need to deliver while pleasing the major customers and avoiding too much penalties from the government. Meanwhile, the energy consumption in CHPs will be further saved to comply with the appeal of China's carbon peak and net-zero carbon discharge.
Here, a mathematical model, which is composed of minimizing users' dissatisfactory degrees and hot water volume (HWV) as two objectives, was proposed after several years of field investigation in CR POWER group in Shenyang. Hundreds of residents' satisfactory degree functions were obtained by collecting and clustering questionnaires from Dadong and Huanggu districts. HA-PBO-CRO was employed to optimize those valve gate angles, which have a direct effect on the users' satisfactory degrees and HWV at each heat exchange station. Through comparing HA-PBO-CRO with NSGAII, it is easy to discover that HA-PBO-CRO is superior to NSGAII in effectiveness and efficiency. And Fig. 1 General operational mode-1 of a CHS [3] the computational result is guidance for CHPs to reach sensible decisions.

Literature Review
CHSs have a history of several hundred years. There are several components in CHSs. Hear source means CHPs that provide hot water in the main pipe. In Northern China, hot water is the main product and electricity is a byproduct that will infuse into state grid directly if no storage technologies were adopted. Today, the more the electricity from coal, the more the fine from Chinese government will be. Cutting hot water volume will be especially critical since it will indirectly reduce the fine for over supplying electricity. Heat network (HN) means a vast extensive underground water pipes going to each corner of several districts. HES is a node on the heat network linking the main heat network and the secondary network responsible for reducing the hot water temperature in the secondary heat network by adjusting heat exchanger's valve angle to allow more or less hot water in the main pipe to flow in. The main heat network links the main pipes and the secondary network links the residential quarters or complexes. The radiators are end devices placed under each household's windows to heat up the whole apartment. Currently, some European countries like Russia, France and England are among the world's leaders in utilizing CHS. In Russia, more than 70% residential buildings in 95% towns are based on CHS [6]. The Russian heating system has its own unique characteristics: high-quality heating equipment, mature technology for the operation and adjustment, and a high degree of flexibility and automation. Similarly, in Denmark, 70% of the population is benefiting from CHS, where one CHP plant could show off a productivity of up to 90%.
CHS in China has a history of 80 years from 1941. The booming period is after the founding of New China, when economic development continued to develop and the construction of infrastructure was in its full pace. Beijing, Shenyang, Taiyuan and other places with the geographical advantage of being close to the coal mines have built a variety of thermal power plants with different scales. These CHP plants contributed a great deal to supply heat to surrounding neighborhoods. According to statistics, CHS accounts for 62.9 percent of the total heat capacity in China [7].
In CHS related researches, scholars mainly focused their attentions on supplying heat on demand and volume adjustment aspects. The mathematical models were single or multiple objectives and the solving tools were mostly metaheuristics or intelligent optimization algorithms.
Supplying heat on demand has always been a research focus of scholars all over the globe. Woznick et al. [8] applied polar bear algorithm to analyze and calculate the effect of outdoor temperature changes in Germany and Finland. They proposed to use the circulating pump in the heat exchange station to implement on-demand heat supply. Similarly, Turski et al. [9] adopted heat load prediction method to achieve dynamic adjustment in heating stations according to climatic conditions and heat users' demands. Zhong et al. [10] put forward the concept of "smart heating" under the background of "Internet + ". They suggested that with big data, artificial intelligence, simulation and other technologies, CHS can make corresponding self-adjustment and self-adaptation according to actual needs. They applied their theory to a CHS established in Beijing hiring online simulation means and intelligent algorithms to operate realtime multi-source heat load dispatching.
Volume adjustment is a key factor in the operation in CHS. Zhao et al. [11] proposed to install regulating valves in the steam heating network. And they calculated the ideal valves' opening angles with parallel genetic algorithm to change the distribution of steam flow in an industrial park in Shanghai. Coal consumption is an important index to measure the cost of CHS. Wang et al. [12] applied particle swarm optimization algorithm to power load distribution with the optimization goal of minimizing the total coal consumption. They proved their theory by comparing the total daily coal consumption data from before and after optimization. Ikonen et al. [13] established a dynamic linear model minimizing the operating cost of regional heating supply system and verified the effectiveness of their model by simulation computation.
As time went by, the focus shifted from single objective problems to multiple objective ones. Ju [14] established a multi-objective stochastic dispatching model minimizing the operating cost and load fluctuation in virtual thermal power plants and proved the effectiveness and practicability of their method. Wei [15] used a multi-objective interval optimization method to solve the synchronized energy capacity allocation optimization problem in CHS and evaluated Pareto optimal solution set from the perspective of cost, risk and emission with dynamic game theory. Zhang et al. [16] put forward a method combining multi-objective optimization with integrated decision-making to solve economical discharge problem considering different heat demands under changing conditions. Chen et al. [17] established a bi-objective model to maximizing economic benefit and minimizing environmental cost in CHS. An improved particle swarm optimization algorithm was adopted to optimize heat load distribution. The overall benefit of the optimized load distribution plan was better than workers' routine operation. He [18] established a stochastic model with the objective function minimizing the total cost of CHS and the total amount of exhaust emissions.
Meta-heuristics are random search methods mimicking some animal behaviors or physical and chemical phenomena. Although Meta-heuristics may not guarantee to find the optimal solutions, they still possess unique and irreplaceable advantages. For example, they could solve problems without requiring the objective functions and constraints to be continuous and convex.
Since Genetic algorithm was proposed in 1975 [19], many novel meta-heuristics and their hybrid versions emerged in abundance. In 1983, Kirkpatrick et al. [20] invented a simulated annealing (SA) algorithm from the annealing process  [21] put forward an ant colony optimization (ACO) algorithm according to the phenomenon of ants foraging process. In 1995, Eberhart and Kennedy [22] devised particle swarm optimization (PSO) following the social behavior of birds. In 2007, Karaboga et al. [22] created artificial bee colony algorithm (ABC) based on honey collection activities. In 2009, Yang [23] coded firefly algorithm (FA) based on the social behavior of fireflies. In 2010, Albert Lam and Victor Li [24] publicized a chemical reaction optimization algorithm (CRO) in view of the intermolecular interaction in chemical reactions. In 2012, Pan [25] described a Fruit Fly Optimization Algorithm (FFOA) considering the foraging behavior of fruit flies. In 2013, Kaveh and Farhoudi [26] exhibited a dolphin echolocation algorithm (DE) according to the dolphin's ability to locate food in the sea. In 2015, Mirjalili [27] developed a dragonfly algorithm (DA) by observing the static and dynamic clustering behaviors of dragonflies. In 2017, a butterfly algorithm appeared from the butterfly mating mechanism [28]. David and Marchi [29] produced a polar bear optimization (PBO) in 2017. Zeng et al. [30] showed a whale swarm algorithm in 2017. Gravitational search algorithm came into being in 2009 based on the Newtonian gravity and the laws of motion [31]. An antlion optimizer (ALO) was discovered in 2015 displaying hunting mechanism of antlions in nature [32]. Harris Hawks Optimization (HHO), which was introduced by Heidari et al. in 2019, featured the cooperative behaviors and chasing style of hawks [33]. Teaching-Learning (TL) algorithm gained attentions from scholars in 2011 carving out the process of teachers' teaching effects and students' learning efficiency [34].
As is known to all, there is a 'no-free-lunch' theory claiming no algorithm is suitable for all optimization problems. After single versions were known to the world, many hybrid algorithms came into being afterwards to improve the performances of single versions. Garg [35] hybrid PSO and GA together for a series of constrained optimization problems in 2016. Then, Garg [36] integerated GA with GSA to solve constrained optimization problems in 2019. In 2021, Kundu and Garg [37] improved TL algorithm and bound it together with HHO for numerical and engineering optimization problems.
All of the abovementioned hybrid algorithms achieved better results. Therefore, we adopted a hybrid algorithm to solve VAAP.

Problem Descriptions and Mathematical Model
Through an extensive investigation on a cogeneration power plant in Shenyang, we found that the main problems in power plants are low customer satisfaction and excessive power supply to the State Grid, which leads to fines of several millions to dozens of millions of RMB.
In the centralized heating system, HESs connect with the primary pipe and the secondary heat network. HESs can have an indirect effect on the water supply temperature and the backwater temperature in the secondary network and the users' indoor temperature according to the different weather conditions and changes of users' demands by adjusting the volume of hot water entering HESs. Adjusting hot water volume means changing the valve disc angle of each HES. Normally, during the valve disc angle adjustment process in a HES, the operators turn the wheel on the valve to increase or reduce the entering water volume according to their own experiences. This empirical behavior could not satisfy the needs of heat users very well. Thus, the complaints are growing and the satisfaction degrees of heat users keep dropping down in recent years.
A HES system is simplified in Fig. 5. Here, E i is the entering water valve of HES i ; i is the opening angle of E i ; t g1 is the entering water temperature in the primary network; t g2 is the entering water temperature in the secondary network; t h1 is the backwater temperature in the primary network, and t h2 is the backwater temperature in the secondary network.
Assume that the heat loss is zero throughout the centralized heating system. Then, the heat Q in the pipe could be calculated with formula (1) [38][39][40][41][42]. Here, V is circulating water volume; C is the specific heat capacity; t g is entering water temperature and t h is backwater temperature.
From the view of heat balance, if the heat loss of hot water in the pipe is neglected, formula (2) will stand.
Here, V 1 is the water volume entering into HES from the primary network and V 2 is the water volume entering into HES from the secondary network. And, the relationship between t g2 and V 1 can be obtained as follows: When the heating system is in a stable state, there is a balanced relationship between the indoor temperature T n gained from the radiator at the user's end and t g2 (see formula (3)) [40][41][42].
Here, KF r is the heat dispersing coefficient of a radiator; KF b is the heat dispersing coefficient of a building and T w is the outdoor temperature. Suppose KF r , KF b , t h2 and T w are constant. Then, the relationship between T n and t g2 is like T n = c ⋅ t g2 + d . Furthermore, we can get In this paper, we consider how to optimize the water supply volume while satisfying the users' satisfaction with the indoor temperatures. We establish fuzzy trapezoidal functions with ten types according to the users' preferences summarized in the questionnaire. By adjusting the water supply volume from the primary network into the HES, the indoor temperature at a user's house can be indirectly affected and the satisfaction degree of this user could be obtained as well. We assume that the opening angle of each valve in the HES systems is 90 • , and only the opening of the entering water valve disc angle of each HES needs to be optimized. In addition, we assume that the indoor temperature is also affected by the distance, namely T nr ij = T n ij − D ij 100 . Here D ij (meters) is the distance between a user j and a HES i.
A bi-objective optimization model is adopted to represent the real operational problem and x i is used as the optimization variable Here, x i represents the ratio of the opening angle i of the inlet valve at the ith HES to 90 degrees. The smaller the x i , the lower the water supply vo l u m e a n d t h e i n d o o r t e mp e ra t u re . T h e n , (4) indicates that the average value of the users' dissatisfaction degrees should be as small as possible. And Formula (5) exclaims that the average number of total water supply volumes of all HESs should be as small as possible. (4) where N is the total number of HESs; M represents the total number of users and FF ijk is the k-th type satisfaction function value corresponding to the indoor temperature at the j-th user in the i-th HES secondary network system.

Multi-Objective Optimization and NSGA-II
In real life, when several objectives are optimized, conflicts will happen. When the solution of one objective is improved, it may harm the other objectives. Multi-objective optimization is an optimization process of finding solutions optimizing all objectives as much as possible. The solution set of multi-objective optimization problem is called Pareto optimal solution set, which was first put forward by a Italian economist Pareto in 1906 [43]. Finding a set of solutions is the ultimate goal of Pareto optimal process. The multi-objective optimization problem (MOP) is described as follows: There are p inequality constraints and q equality constraints.

Multi-Objective Important Concepts and Definitions [43]
Definition 1(Pareto domination x≺y): Definition 5 (Pareto optimal solution set): A set of all Pareto optimal solutions is a Pareto optimal solution set.
Definition 6 (Pareto frontier): The objective values of all solutions in Pareto solution set in a coordinate space constitute a Pareto frontier. Figure 6 clearly describes the relationship between the dominant solutions and the Pareto frontier for a bi-objective optimization problem.

NSGA-II
With the continuous development of multi-objective research, many multi-objective optimization algorithms emerge. NSGA-II is one of those multi-objective optimization algorithms with ideal results. Thus, it is often used as a benchmark for comparing the calculation results of other algorithms. NSGA-II means non-dominated sorting genetic algorithm with elite mechanism. It was put forward by Deb et al. [44] based on NSGA in 2000. NSGA-II overcomes three shortcomings of NSGA: (1) optimal solutions retention mechanism; (2) adaptive setting for sharing parameters; (3) computation complexity [45][46][47][48]. In NAGA-II, the evolutionary population is divided into several layers by dominance relations. The first layer is the set of non-dominant individuals and the second layer is the set of non-dominant individuals obtained by removing those individuals from the first layer [46]. The advantages of NAGA-II algorithm are as follows: (1) Introducing the elite retention mechanism to expand the population space, and combining the parent and child generations to compete together to produce the next generation population preventing the optimal solution from being missed out; (2) The non-dominated sorting speed is improved by reducing the computational complexity from O(n 3 ) to O(n 2 ); (3) Crowding degree can keep the distribution and diversity of solution population and can avoid sharing parameter selection [43].
The flowchart of NSGA-II is as shown in Fig. 7.

Polar Bear Optimization
Polar Bear Optimization (PBO), which is one of the bioinspired meta-heuristic algorithms, was proposed by David et al. in 2017 [29]. The idea of the algorithm came from the foraging process of polar bears (see Fig. 8).
The core of polar bear algorithm is composed of ice floe floating strategy, seal predation strategy and bear colony breeding and death strategy.
Ice floe floating strategy is a global search mechanism which can be represented by formula (9). Here, β is a random number in the interval (0, 1) and γ is a random number in the interval [0,|ω|]. sign(ω) is a piecewise function with a result {-1, 0, 1} changing with ω {negative, 0, positive}. ω, which is depicted in formula (10), is the Euclidean distance between two polar bears.
Seal predation strategy serves as local search mechanism. The hunting mode of each polar bear can be expressed by formula (11). In the searching process, the radius of view is denoted by the following two parameters: a ∈ [0, 0.3] limits the visible distance; � 0 ∈ 0, π 2 indicates the oblique angle at which a polar bear travels around prey.

Fig. 6 Pareto frontier and dominant solutions
If the updated position of a polar bear is better than the current position, " + " sign is applied to indicate a forward move; if the updated position is not ideal, a polar bear will go in the direction of "−". When the results in both directions are not acceptable, the polar bear stays where it is.
Bear colony breeding and death strategy, which is a population updating mechanism, is indicated by formula (13), (14).
In formula (13), k is a random number in the interval [0, 1]. If the death mechanism is triggered, a weakest polar bear will die if the population number is more than 50% of the maximum population capacity n. In formula (14), in the t-th iteration, one of top 10% individuals in the population will reproduce and the population will expand if the population size is less than n-1.
The executive steps of PBO are as follows: Step.1 Setting parameters; Step.2 Randomly generate 75% of population capacity n and record the best solution; Step.3 For each polar bear, randomly set the angle of each dimension; Step.4 Use formula (11) to calculate the search radius r; Step.5 Use " + "to calculate the new position of a polar bear by formula (12) and update the current position if the new position is better than the current position; Step.6 Calculate the new position of a polar bear by formula (12) in the direction of "−" and update the current position if the new position is better; Step.7 According to formula (9) and ω, determine the new position of each polar bear; Step.8 Sort the population and select the first polar bear to judge whether to update the global best solution; Step.9 Randomly select a solution different from the global best solution from top 10% in the current population; Step.10 Randomly generate k. if k < 0.25 and the population size is greater than 50% of the population capacity n, remove the solution that ranked last; Step.11 If the population size is less than n-1, generate a new bear to join the population according to formula (14); Step.12 Judge whether the maximum iteration number has been reached. if yes, go to step.13, if no, go to step.3; Step.13 OuCHPut the global best solution and fitness value.

Chemical Reaction Optimization [24]
Chemical reaction optimization (CRO) was designed in 2010 [24]. CRO mimics the process that molecules react in a closed container and finally reach a low-energy stable state. In CRO, a molecular structure ( ) represents a candidate solution. The potential energy (PE) of a molecule represents the objective function value f (see formula (15)). Kinetic energy (KE) determines the vitality of molecules and represents the tolerance of accepting a molecule worse than the current one. The greater the KE is, the stronger the ability to escape from local minimum will be. The molecular elements are composed of φ, PE, KE, Number of total hits (number of hits), current best solution (minimum structure), current best objective function value (minimum value) and number of hits when best solution is achieved (minimum hit number). In a chemical reaction, kinetic energy and potential energy can be mutually transformed into each other within or between molecules. According to the energy conservation law, if chemical reactions take place in a closed system, energy can only be transferred from one form to another or from one molecule to another. And the total energy will not change. Therefore, in CRO algorithm an energy buffer (EB) is used to store part of the kinetic energy in a molecule. A molecule is allowed to change into ′ if PE ≥ PE ′ or PE + KE ≥ PE �.
In on-wall ineffective collision, Molecules bounce back after hitting the container wall. During collision, the structure of a molecule may remain unchanged or be changed. The molecular structure changes into ′ after the collision, if the relationship between energies satisfies formula (16). The new KE ′ can be obtained with 1] and (1-q) represents the loss of KE when collision happens. KElossrate is a parameter whose value satisfies 0 < KElossrate < 1. The lost KE will be stored in the central EB, which will support decomposition. The more the collision happens, the more KE is stored in the central EB. If formula (16) is not met, a molecule will keep its original structure, PE and KE.
Decomposition starts when a molecule hits the container wall and decomposes into two independent molecules. In CRO, decomposition reaction needs to satisfy formula (17). The KEs of two new molecules are calculated using formula (18), (19) with k ∈[0, 1]. If formula (16) is not met, the KE stored in the central EB can be used to support the decomposition reaction according to formula (20). The KEs of two new molecules are then calculated by formula (21), (22) with (m 1 , m 2 , m 3 , m 4 )∈[0,1]. If formula (20) is not true, the decomposition reaction will fail to happen. Fig. 9 Four basic reactions in CRO [24] The inter-molecular ineffective collision occurs when two molecules collide with each other and then bounce off. The KEs in this reaction will not go into central EB.
Assume that the original molecular structures are 1 and 2 . During the collision, two new molecular structures ′ 1 and ′ 2 are obtained from the neighborhood when formula (23) is satisfied. Two KEs of new molecules are calculated using formula (24) and formula (25) with p∈ [0, 1].
Synthesis reaction forms when two molecules collide and merge into one molecule. Assume that the molecular structures of two molecules are 1 and 2, respectively. During the collision, the existing two molecular structures synthesized into one molecular structure ′ when formula (26)  is satisfied. And the KE of the new molecule is calculated using formula (27).
Looking through all four reactions, they are just criteria judging if new solutions should derive from old ones. However, there are no exact rules about how a solution should be changed into a new one. Therefore, CRO serves as a flexible framework accepting any possible detailed solution transformation operations. The simplified flowchart of CRO is shown in Fig. 10.

A Hybrid Algorithm of Polar Bear Optimization and Chemical Reaction Optimization
In this paper, the adjustment decision of HWV considering customer satisfaction is considered and a hybrid algorithm of polar bear optimization chemical reaction optimization (HA-PBO-CRO) is applied to optimize a multi-objective optimization model.

Modified On-Wall Ineffective Reaction
In this part, ice floe floating strategy in PBO is adopted to generate a new solution.

Modified Decomposition
Here, seals predatory strategy in PBO is applied to produce two other solutions.

Modified Inter-Molecular Ineffective Collision
In GA, crossover operation is often performed to turn two solutions into two different ones. Therefore, a crossover operator should fit well with this type of collision reaction.

Simulation
Here  Tables 2, 3.  According to Tables 2, 3, we can see that there is always an imbalance between heat supply and actual need. Insufficient supply incurs dissatisfactions and overheating leads to a dramatic waste of heat and energy resources.
China's adjustment on indoor temperature in winter is between 18 and 22 °C. We collected 50 indoor temperatures in a building on a winter day in 2018-2019 winter and listed them in (Table 4).
Seeing from Table 4, we can know that 26% users are having temperatures below 18 °C resulting in an ever-growing complaint rate each year.
Here, we consider a primary network with N HESs and each HES serves M heat users. We optimized the inlet valve disc angles of HESs with a bi-objective hybrid algorithm based on polar bear optimization and chemical reaction optimization. As we have mentioned, adjusting valve dis angle will affect the indoor temperature and indirectly have an    [24,25] 0 t i ∉ [18,25] effect on energy consumption and online electricity volume. Therefore, a handsome save on operational cost in a CHP could be made.   Hypervolume (HV) is a widely used criterion for evaluating convergence and solution distribution of a multi-objective algorithm in recent years (see Fig. 12) [49]. According to Fig. 12, we can draw a conclusion that the higher the HV, the better the Pareto frontier.
Some parameters for HA-PBO-CRO are listed in Table 5. And other parameters used in HA-PBO-CRO are set according to simple initial tests and the experiences reached in literature [24] and [29]. In HA-PBO-CRO, PE = (f1+f2) 2 . Parameters of NSGA-II are set according to literature [46]. The data about the distances between heat users and their HES are shown in Appendix.
Furthermore, we concluded ten types of satisfaction fuzzy trapezoidal functions based on investigations on 3000 local heat users (see Table 6).
The best computational results of HA-PBO-CRO and NSGA-II after 100, 200, 300 and 500 iterations in 20 runs with 10 HESs (each has 30 users) are presented below in Fig. 13.
From the comparison of Pareto frontiers in Fig. 13, we can see that with the increase of iterations, the Pareto frontier integrity of HA-PBO-CRO is obviously better than that of NSGA-II algorithm. When the iterations are small, NSGA-II algorithm can not even give a complete Pareto frontier. Therefore, HA-PBO-CRO has an obvious advantage over NSGA-II algorithm in dealing with heat supply problem. After 20 independent runs, the average HV index is shown in Table 7.
According to the comparison results of hypervolume index value (HV), it can be seen that the solution set quality produced by HA-PBO-CRO is better than that of NSGA-II algorithm.
In order to further evaluate the performance of HA-PBO-CRO, we increase the number of HES to 20 and the number of heat users to 600 with each HES serving 60 users. The Pareto frontier comparison results are presented in Fig. 14 and the HV values are listed in Table 8.
Two computational experiments show that HA-PBO-CRO is clearly superior to NSGA-II in solution diversity, solution distribution and convergence. And HA-POB-CRO only takes two-thirds of the computing time of NSGA-II to produce a result.

Conclusions
In the 8-month-long cold environment, CHP plants in northern China are vital factors for local residents' livelihoods. Through on-the-spot investigation, field research and communication with the management level of a large coal-fired CHP plant in Shenyang, we found that what is bugging them the most is the excessive electricity generated during the heat supplying process. It seems that the CHP plants should be profitable, but as a fact the fines from the government have caused serious losses in revenue. Thus, the management enquired about whether it was possible to optimize the water supplying volume with operations optimization technology to reduce the electricity entering the State Grid. Therefore, we carried out theoretical research by analyzing and examining the DHS carefully. First, we transformed the optimization on water supply volume into the optimization on inlet valve disc angle adjustment of HESs in DHS. Then, we established a bi-objective optimization model minimizing customer dissatisfaction degree and water supply cost. And we adopted a hybrid meta-heuristic algorithm to solve the mathematical model. By comparing the computational result with that of a NSGAII, we exhibited that the proposed algorithm in this paper has advantages over NSGAII in solution quality and computation efficiency. Moreover, the computational result provided significant insights for CHP management.
The research limitation in this paper is that the feasibilities of assumptions and objectives in the mathematical model need to be further verified in real-world settings. And, it is necessary to carry out more extensive and large amount of calculation tests to verify the performance and parameters of the proposed algorithm in this paper.
Future research should focus on developing online realtime heat network monitoring and operational optimization software. Therefore, our algorithm could be integrated into optimization module in the software easily as one of the optimization tools. In addition, some other interesting topics revolving CHPs are worth caring about as well. For example, research on applying Ant Lion Optimizer to work out reasonable supplier Side Bidding Strategies at Day-ahead Electricity Market. And, study on hiring whale optimization algorithm for robots' patrolling path perception and optimization in a highly automated CHP.