Optimal Chiller Loading for Energy Conservation Using an Improved Fruit Fly Optimization Algorithm

In the multi-chiller of the air conditioning system, the optimal chiller loading (OCL) is an important research topic. This research is to find the appropriate partial load ratio (PLR) for each chiller in order to minimize the total energy consumption of the multi-chiller under the system cooling load (CL) requirements. However, this optimization problem has not been well studied. In this paper, in order to solve the OCL problem, we propose an improved fruit fly optimization algorithm (IFOA). A linear generation mechanism is developed to uniformly generate candidate solutions, and a new dynamic search radius method is employed to balance the local and global search ability of IFOA. To empirically evaluate the performance of the proposed IFOA, a number of comparative experiments are conducted on three well-known cases. The experimental results show that IFOA found 14 optimal values (the optimal values among all algorithms) under a total of 17 CLs in three cases, and the ratio of the optimal values found was 82.4%, which was the highest among all algorithms. In addition, the mean value of all objective functions of IFOA is smaller and the standard deviation is equal to or close to 0, which proves that the algorithm has high stability. It can be concluded that IFOA is an ideal method to solve the OCL problem.


Introduction
To maintain comfort level of an indoor life during hot and humid weather, many people rely on air-conditioning systems. Thus, air-conditioning systems are used in large quantities, and their energy consumption for the supply of CL accounts for more than 30% of the total power generation [1]. In air-conditioning systems, the multi-chiller system is the main energy-consumption equipment. If the chillers are improperly managed, the energy consumption of the chiller will be increased significantly. Because the multi-chiller system is composed of chillers of different design capacities and performance feature, the optimal load distribution of each chiller can obtain minimum energy consumption when meeting CL demand [2]. Therefore, for the part load ratios of all chillers, it is a valuable research topic to find their optimal combination by using optimization methods.

Problem Description
Due to the large demand for cooling in massive buildings, the multi-chiller system is generally used. Figure 1 depicts the structure of a typical multi-chiller decoupled system. It consists of a primary side (chiller side) and a secondary side (load side) [2]. On the primary side, multiple chillers are connected to the distributed system in series or parallel. The system allocates different loads to each chiller by controlling the supply and return water flow because each chiller has different capacities. On the secondary side, the cold water flowing into the cooling coil can be adjusted by a two-way valve according to load changes. For a given wet-bulb temperature, the power consumption (P) of a centrifugal chiller is a convex function of its PLR [6]. Pi is expressed as a polynomial of PLRi, as shown in Equation (3): Here, c1i, c2i, c3i, c4i are constants, representing the coefficients of the ith chiller KW-PLR curve, The multi-chiller system is designed to meet peak loads, but most of the time, the chiller operates at partial load ratio. PLR i can be calculated using Equation (1) [2]. Here, CL i represents the ith chiller CL, and RT i represents its design capacity (RT). Equation (1) is the expression of the meaning of PLR, but we do not use Equation (1) to calculate the value of PLR. The value of PLR is obtained through the IFOA iterative optimization proposed in the paper. Obviously, the value range of the PLR should be in the range [0,1], but when PLR is small, the chiller is prone to surge, thus, the manufacturer recommends that the PLR of each chiller should be greater than or equal to 30%, as shown in Equation (2). Here ON i is the state of the ith chiller. When ON i = 1, it is on; when ON i = 0, it is off. In a multi-chiller system, in order to minimize energy consumption, one or several chillers are allowed to not start up, that is, their PLR is 0. Equation (2) is the constraint condition that the OCL problem must satisfy: For a given wet-bulb temperature, the power consumption (P) of a centrifugal chiller is a convex function of its PLR [6]. P i is expressed as a polynomial of PLR i , as shown in Equation (3): Here, c 1i , c 2i , c 3i , c 4i are constants, representing the coefficients of the ith chiller KW-PLR curve, respectively. Cases 1, 2, and 3 correspond to the case study with six, four, and three chillers in Section 6, respectively. Finding the minimum value of the total energy consumption is the objective function of the OCL problem in the multi-chiller system, as given in Equation (4): where P i represents the ith chiller power consumption and n is the total number of chillers.
The sum of the CLs generated by all the chillers in the multi-chiller system should not be less than the system CL demand. This constraint must be satisfied, and can be denoted using Equation (5): Under the condition of satisfying CL requirements, if the total energy consumption of all chillers is minimized, then the performance of multi-chiller system is the best [15]. The functional objective of the OCL problem is to minimize the total energy consumption of multi-chiller, that is, the value of J obtained in Equation (4) is the minimum. The constraints of the OCL problem are Equation (2) and (5), that is, the PLR of each chiller must be greater than or equal to 0.3, and the CL generated by multi-chiller must satisfy the system CL demand. The solution to the OCL problem is the appropriate PLR value for each chiller in a multi-chiller.

Canonical FOA Overview
Canonical FOA is the basic version of FOA proposed by Pan in 2011. In order to show the difference from the various FOA versions improved later, we generally call it "canonical FOA". The canonical FOA adopts a global random search strategy based on swarm. In each generation, there are two phases of osphresis and vision foraging. In the osphresis foraging phase, each individual in the fruit fly swarm searches randomly around swarm, and then finds the maximum smell concentration of the individual after evaluating the new location. In the vision foraging phase, all individuals fly to the location where smell concentration is maximum. The canonical FOA has four steps.
Step 1. Initialization. Set parameters such as the population size (PS), the maximum number of iterations (Iter max ), the random fly direction and distance zone of fruit fly (FR), and the fruit fly swarm location range (LR). The location of fruit fly individual in the swarm is given by its corresponding two-dimensional coordinates (X, Y), whose initial location is defined by Equation (6): Here the function of rand(LR) is to get a number arbitrarily within the range of the positions of the fruit fly swarm.
Step 2.1. Using the osphresis organ, fruit fly individuals search randomly around the swarm. Each individual is given a flight direction and distance randomly. The new location of the individual is computed by Equation (7): Here the function of rand(FR) is to get a number arbitrarily within the fly range of the fruit fly.
Step 2.2. The distance (DIST i ) between the individual and the origin is calculated using Equation (8) due to the exact location of the food being unknown. Then, the smell concentration judgment value (S i ) of a fruit fly individual is computed by Equation (9), which is the reciprocal of the distance: Step 2.3. The smell concentration (Smell i ) of each fruit fly individual in the swarm is calculated by Equation (10), that is, S i is substituted into the fitness function (or smell concentration judgment function): Step 2.4. Find the fruit fly with the best smell concentration in the swarm and record its smell concentration and corresponding location, as shown in Equation (11): Step 3. Vision foraging phase. The best smell concentration value and corresponding fruit fly location are maintained by Equations (12) and (13), respectively, and other fruit flies will use vision to fly toward that location: Energies 2020, 13, 3760 5 of 18 Step 4. Repeat steps 2 and 3 until Iter reaches Iter max .

Disadvantages of Canonical FOA
By analyzing the canonical FOA algorithm, two disadvantages of FOA can be summarized as follows.
(1) Nonuniform generation of candidate solutions Shan et al. [42] first substituted Equation (7) into Equation (8). In Equation (8), X _axis and Y _axis are constants, and rand(FR) is a variable. Then they used the counter-proof method to prove that for the given fruit fly position coordinates X and Y, the taste concentration judgment value S i does not follow a uniform distribution. Usually, S i is the solution to the optimization problem. Obviously, the candidate solutions cannot be generated uniformly, that is, FOA loses the ability to search uniformly in the solution space.
(2) Poor search ability The search radius of the fruit fly is its range of flight (FR). When the FR is large, the search range of fruit fly is large, and the global search ability of FOA is strong, while when the FR is small, the search range of fruit fly is small, and the local search ability of FOA is strong. In the canonical FOA, FR is a fixed value, which cannot better balance the global search ability and local search ability of the algorithm.

Improved FOA Algorithm Based on Dynamic Search Radius
To overcome the first disadvantage summarized above, we should first find the cause of the disadvantage. An analysis of Equations (6)- (9) found that the location of an individual can be represented by two-dimensional coordinates, and the candidate solution is the reciprocal of the distance between individual location and the origin. The candidate solutions generated by FOA are nonlinear, thus, they fail to follow the uniform distribution.
Inspired by the reference literature [42], for FOA, we can change the non-linear of the candidate solution generation mechanism into a linear one, that is, the location of the individual is represented by one-dimensional coordinates, and we no longer calculate the reciprocal of the distance between individuals and the origin. Thus, the original equations have been modified. That is, Equations (6), (7), (9), and (13) are changed into Equations (14)- (17), respectively, and Equation (8) is deleted.
The advantage of using this linear generation mechanism is that candidate solutions can be generated uniformly in the solution space, so that the swarm has the ability to search the global optimal solution: For the second disadvantage mentioned above, this paper proposes a new dynamic search radius method that changes the search radius by the iteration number, as shown in Figure 2. Therefore, the search radius can cover the entire solution space in early iterations, which gives the algorithm a good global search capability. In later iterations, the search radius becomes very small when the swarm Energies 2020, 13, 3760 6 of 18 location is close to the best solution, which guarantees that the algorithm will have a good local search capability. The search radius r can be obtained by Equation (18): where r min is the minimum search radius, r max is the maximum search radius, Iter max is the maximum number of iterations, and Iter is the number of iterations.

Figure 2. Variation of r versus iterations.
A good initial swarm location can make the algorithm converge faster. Therefore, in the improved FOA, a population is generated randomly, and then, the location of the optimal individual is selected as the initial swarm location instead of the location generated according to Equation (14). The improved FOA algorithm based on the dynamic search radius is called IFOA. Algorithm 1 shows the main structure of IFOA. A good initial swarm location can make the algorithm converge faster. Therefore, in the improved FOA, a population is generated randomly, and then, the location of the optimal individual is selected as the initial swarm location instead of the location generated according to Equation (14). The improved FOA algorithm based on the dynamic search radius is called IFOA. Algorithm 1 shows the main structure of IFOA.

Algorithm 1. Improved fruit fly optimization algorithm (IFOA) algorithm
Parameters: PS, Iter max , r max , r min Output: Solution X* //Initialization Set PS, Iter max , r max , r min

Implementation of IFOA on OCL Problem
When using IFOA to solve the OCL problem, each individual in the fruit fly swarm corresponds to a solution in the problem. Taking a multi-chiller composed of three chillers as an example, the PLR value of each chiller is represented by a decimal code value, and the encoded values of the three chillers constitute a fruit fly individual, as shown in Figure 3, where each PLR value must satisfy the constraints given by Equation (2).

Implementation of IFOA on OCL Problem
When using IFOA to solve the OCL problem, each individual in the fruit fly swarm corresponds to a solution in the problem. Taking a multi-chiller composed of three chillers as an example, the PLR value of each chiller is represented by a decimal code value, and the encoded values of the three chillers constitute a fruit fly individual, as shown in Figure 3, where each PLR value must satisfy the constraints given by Equation (2). Step 1. Initialization. First determine the center of a good fruit fly swarm. Randomly generate an initial swarm of fruit fly composed of PS fruit fly individuals, calculate the energy consumption of each chiller according to Equation (3), and then obtain the total energy consumption of multi-chiller system by Equation (4). The energy consumption of all PS multi-chiller systems is calculated accordingly. After substituting the PLR values of the three chillers into Equation (5), the total cooling output is obtained and compared with the CL demand of the system. The value of the fitness function of the multi-chiller that satisfy the constraints is equal to its energy consumption value, and the value of the fitness function of the multi-chiller that does not satisfy the constraints is given a larger penalty value. Find the fruit fly with the smallest fitness function value in the swarm, and take the position of this fruit fly as the center position X_axis of the swarm.
Taking Figure 3 as an example, the data of the three chillers are in Ref. [6], and their rated cooling capacity is 800RT. The calculation process is as follows:  Step 1. Initialization. First determine the center of a good fruit fly swarm. Randomly generate an initial swarm of fruit fly composed of PS fruit fly individuals, calculate the energy consumption of each chiller according to Equation (3), and then obtain the total energy consumption of multi-chiller system by Equation (4). The energy consumption of all PS multi-chiller systems is calculated accordingly. After substituting the PLR values of the three chillers into Equation (5), the total cooling output is obtained and compared with the CL demand of the system. The value of the fitness function of the multi-chiller that satisfy the constraints is equal to its energy consumption value, and the value of the fitness function of the multi-chiller that does not satisfy the constraints is given a larger penalty value. Find the fruit fly with the smallest fitness function value in the swarm, and take the position of this fruit fly as the center position X_axis of the swarm. Taking Figure 3 as an example, the data of the three chillers are in Ref. [6], and their rated cooling capacity is 800RT. The calculation process is as follows: Step 2. Osphresis foraging phase. Find the search radius r according to Equation (18), then use the equation X i = X _axis + r·rand() to find the position of the fruit fly, and a swarm of PS fruit fly individuals is generated according to this method. Next, calculate the total energy consumption, total cooling output, and fitness function values of the multi-chiller separately. The calculation method is the same as Step 1.
Find the fruit fly with the smallest fitness function value in the swarm, and take the position of this fruit fly as the center position X_axis of the swarm.
Step 4. Repeat steps 2 and 3 until Iter reaches Iter max . After the iteration is terminated, the coding of the fruit fly individual with the smallest fitness function value is the solution to the OCL problem.

Cases Used in Experiments
In the experiments, three well-known cases are selected to verify the performance of IFOA in solving the OCL problem.

Case with Six Chillers
Case 1 is based on the multi-chiller system of a semiconductor plant located in Hsinchu Science Garden (Taiwan), and was originally proposed by Ref. [6]. It consists of six chiller units. Among them, the capacity of four units is 1280RT, and the other two units are 1250RT. Table 1 lists the chiller data for the first case that can be used in Equation (3). Case 2 is based on the multi-chiller system in a hotel located in Taipei, and was originally proposed by Ref. [3]. It consists of four chiller units, among which two units have the capacity of 450RT, and two units have the capacity of 1000RT. Table 2 lists the chiller data for the second case that can be used in Equation (3).  The multi-chiller system in Case 3 is also a semiconductor factory located in Hsinchu Science Garden, which uses three chillers with a design capacity of 800RT [6]. Table 3 lists the chiller data for the third case that can be used in Equation (3).

Results and Analysis
The IFOA algorithm was implemented in the Visual Studio 2010 environment using C++ and run on the Intel Core i5-9300H 2.40 GHz PC. IFOA was run 30 times independently considering that unexpected situations may occur. The maximum, minimum, mean, and standard deviation were calculated from 30 optimal objective function values to comprehensively evaluate the IFOA algorithm's ability to solve the OCL problem.
The optimal results of IFOA (the minimum of 30 optimal objective function values) are compared with those of some algorithms. We bold the optimal value of the algorithms and the reduction in energy consumption of IFOA.
IFOA has four control parameters, and their setting values are shown in Table 4.  Table 5 summarizes the comparison of the optimal values of IFOA and TLBO [21], Two Stage DE [14], and DCEDA [28]. As can be seen from the  The PLR value here retains six digits after the decimal point, and the PLR is a double type variable in the actual program, thus, the accuracy is reduced, resulting in a slight error in the calculation result. In other words, the J value obtained in the actual program is 3838.2079. In addition, we noticed that the value of PLR 3 is 0.000000, that is, the third chiller is off, and its energy consumption should be 0. But Two Stage DE substitutes 0.000000 into the equation and obtains a value of −120.505. That is, the energy consumption of the third chiller is −120.505 KW, which is obviously unreasonable, and its actual optimal value should be 3958.7129 KW (3838.2079 + 120.505 = 3958.7129); (3) when CL is 6477, IFOA is equal to TLBO; when CL is 6858 and 6096, TLBO saves 0.035 KW and 0.066 KW more than IFOA; when CL is 5717 and 5334, IFOA saves 62.147 KW and 96.073 KW more than TLBO. Figure 4 is a histogram representation of the results in Case 1.
The optimal value of IFOA is the minimum value of the objective function obtained by running it 30 times independently. To further evaluate the stability of the algorithm, Table 6 summarizes the maximum, minimum, mean, and standard deviation of the optimal objective function values got by IFOA and DCEDA after 30 independent runs. Measure the stability of the algorithm by first comparing the mean, and if the mean is equal, compare the standard deviation. It can be seen from Table 6 that under all five CLs, the mean and standard deviation of IFOA are smaller compared to DCEDA. Therefore, we can conclude that the stability of IFOA in case 1 is better than DCEDA. In addition, we also give the running time of the algorithm in Table 6, which is the average time of 30 independent runs. It can be seen from the table that CPU time of IFOA and DCEDA is about 0.7 s for all five CLs.  The optimal value of IFOA is the minimum value of the objective function obtained by running it 30 times independently. To further evaluate the stability of the algorithm, Table 6 summarizes the maximum, minimum, mean, and standard deviation of the optimal objective function values got by IFOA and DCEDA after 30 independent runs. Measure the stability of the algorithm by first comparing the mean, and if the mean is equal, compare the standard deviation. It can be seen from Table 6 that under all five CLs, the mean and standard deviation of IFOA are smaller compared to DCEDA. Therefore, we can conclude that the stability of IFOA in case 1 is better than DCEDA. In addition, we also give the running time of the algorithm in Table 6, which is the average time of 30 independent runs. It can be seen from the table that CPU time of IFOA and DCEDA is about 0.7 s for all five CLs.    Table 7 summarizes the optimal values of IFOA and TLBO [21], Two Stage DE [14], and DCEDA [28]. As can be seen from the table, (1) when CL is 1160, IFOA saves energy by 0.02 KW more than DCEDA, otherwise the two are equal; (2) when the CL is 2320, 1740, 1450, and 1160, respectively, IFOA saves energy by 0.073 KW, 10 Figure 5 is a histogram representation of the results in Case 2.   Table 8 summarizes the maximum, minimum, mean, and standard deviation of the optimal objective function values got by IFOA and DCEDA after 30 independent runs. As seen from Table 8, in all six CLs, the mean values of the optimal objective function values of IFOA are equal to the minimum values, and the standard deviations are 0; this result proves that the stability of IFOA in Case 2 outperforms that of DCEDA. In addition, the average running times of the algorithms are given in Table 8. As can be seen from the table, the CPU time of the algorithm decreases with the decrease of CL. The CPU time of IFOA and DCEDA is approximately equal.   Table 8 summarizes the maximum, minimum, mean, and standard deviation of the optimal objective function values got by IFOA and DCEDA after 30 independent runs. As seen from Table 8, in all six CLs, the mean values of the optimal objective function values of IFOA are equal to the minimum values, and the standard deviations are 0; this result proves that the stability of IFOA in Case 2 outperforms that of DCEDA. In addition, the average running times of the algorithms are given in Table 8. As can be seen from the table, the CPU time of the algorithm decreases with the decrease of CL. The CPU time of IFOA and DCEDA is approximately equal.  Table 9 summarizes the comparison of the optimal values of IFOA and TLBO [21], Two Stage DE [14], and DCEDA [28]. It can be seen from the  Figure 6 is a histogram representation of the results in Case 3.  Table 10 summarizes the maximum, minimum, mean, and standard deviation of the optimal objective function values obtained by IFOA and DCEDA after 30 independent runs. As seen from Table 10, in all six CLs, the mean values of the optimal objective function values got by IFOA equal the minimum values, and the standard deviations are 0; this result proves that the stability of IFOA in Case 3 outperforms that of DCEDA. In addition, the average running times of the algorithms are given in Table 10. As can be seen from the table, the CPU time spent by IFOA and DCEDA is around 0.1 s for all six CLs.  Figure 6. Histogram representation of results in case 3. Table 10 summarizes the maximum, minimum, mean, and standard deviation of the optimal objective function values obtained by IFOA and DCEDA after 30 independent runs. As seen from Table 10, in all six CLs, the mean values of the optimal objective function values got by IFOA equal the minimum values, and the standard deviations are 0; this result proves that the stability of IFOA in Case 3 outperforms that of DCEDA. In addition, the average running times of the algorithms are given in Table 10. As can be seen from the table, the CPU time spent by IFOA and DCEDA is around 0.1 s for all six CLs.   Table 11 summarizes the results of the previous three case studies. As can be seen from the table, in Cases 1, 2, and 3, the number of optimal values (the optimal values in all algorithms) found by IFOA are 3, 5, and 6, respectively, and the ratio of the optimal values found is 60%, 83.3%, and 100%, respectively. In three cases, of a total of 17 CLs, IFOA found a total of 14 optimal values; the ratio of the optimal value found is 82.4%, and the ratio of the optimal value found in all algorithms is the highest. In addition, in order to evaluate the stability of IFOA, by comparison with DCEDA, the average value of all the objective function values of IFOA is smaller, and the standard deviation is close to or equal to 0. It can be concluded from the comparison of the above two aspects that IFOA is superior to other comparison algorithms in solving OCL problems. Table 11. Statistics of the optimal results found by IFOA.

Conclusions and Future Work
OCL is a crucial optimization problem in many realistic applications. However, this problem has not been solved well. Therefore, an efficient optimization algorithm is urgently needed to solve the OCL problem. In this study, in order to solve the OCL problem, in view of the two disadvantages of FOA, we propose an improved FOA algorithm (IFOA). The main contributions of the research are summarized as follows: (1) a mechanism for generating candidate solutions is developed. This mechanism changes the generation of candidate solutions from nonlinear to linear, so that IFOA has the ability to search uniformly in the solution space; (2) a new dynamic search radius method is proposed, which makes the search radius change smoothly with the number of iterations. In the early iterations, the large search radius gave IFOA global search ability. Then, as the number of iterations increases, the search radius gradually decreases. In the later stage of the iteration, the position of the fruit fly swarm is close to the optimal solution. At this time, a small search radius can enhance the local search ability of IFOA.
In order to verify the ability of IFOA to solve OCL problems, we selected three well-known cases to our study. The experiment results show that IFOA has two advantages: (1) strong optimization ability. The probability of IFOA finding the optimal solution is 82.4%, which is the highest among all algorithms; (2) high stability. The mean value of all objective functions of IFOA is smaller and the standard deviation is equal or close to 0. From the above two aspects, we can draw a conclusion that the IFOA algorithm is a highly recommended effective algorithm for solving OCL problems, and it can also be used in other optimization fields.
There are several opportunities for future research on optimization algorithms for the OCL problem. First, the multi-chiller load balancing problem will be considered. Then, some new multi-objective evolutionary optimization algorithms [43][44][45][46][47][48][49][50] can be used to solve the above OCL problems.