Applying Two-Stage Differential Evolution for Energy Saving in Optimal Chiller Loading

In Taiwan, over 45% of the energy in common buildings is used for the air-conditioning system. In particular, the chiller plant consumes about 70% of the energy in air-conditioning system. The electric energy consumption of air-condition system in a clean room of semiconductor factory is about 5–10 times of that in a common building. Consequently, the optimal chiller loading in energy saving of building is a vital issue. This paper develops a new algorithm to solve optimal chiller loading (OCL) problems. The proposed two-stage differential evolution algorithm integrated the advantages of exploration (global search) in the modified binary differential evolution (MBDE) algorithm and exploitation (local search) in the real-valued differential evolution (DE) algorithm for finding the optimal solution of OCL problems. In order to show the performance of the proposed algorithm, comparison with other optimization methods has been done and analyzed. The result shows that the proposed algorithm can obtain similar or better solution in comparison to previous studies. It is a promising approach for the OCL problem.


Introduction
The air-conditioning systems of large-scale commercial buildings in Taiwan account for about 32%-54% of Taiwan's electrical energy consumption, and chiller plants consume more than 70% of the overall energy consumed by air-conditioning systems [1].In systems where multiple chillers are operated in parallel (multi-chiller systems), each chiller can operate independently; adjusting the chiller operation schedule to provide the venue with a stable refrigeration ton (RT) load and a flexible maintenance schedule [1] is a common practice in large commercial buildings.Because multi-chiller systems are composed of chillers of varying features or even of various types of chillers, the question of how to adjust appropriate numbers of operating chillers and operational control points so that each chiller operates at optimal efficiency is crucial for saving energy in air-conditioning systems.
In recent years, many articles have discussed the optimal chiller loading (OCL) problem.Chang [2,3] proposed using a branch-and-bound method and a Lagrangian multiplier method to solve OCL problems.In addition to traditional numerical methods, numerous heuristic optimization methods have been used to solve OCL problems.Optimal chiller loading that consumed less energy than that of Chang [2,3] was obtained by genetic algorithms (GAs) [4].Chang et al. [5] applied an evolutionary strategy (ES) to OCL problems and found a lower power consumption with higher precision for chillers than found previously.Particle swarm optimization (PSO), which uses the

Introduction to Multi-chiller System
Multi-chiller systems can provide flexible operation, reserve capacity, and less frequent system shutdowns for maintenance.Those systems composed of two or more chillers are widely used in the air-conditioning systems of large buildings.In a multi-chiller system, each chiller can operate independently and provide different refrigeration capabilities; the chillers operate efficiently according to different or similar performance curves to meet a wide range of RT requirements in HVAC (heating, ventilation and air conditioning) system.The architecture of the multi-chiller system is as shown in Figure 1 below.Generally, the maximum capacity of a chiller is designed to meet the maximum peak load demand, but because of actual venue requirements and change of seasons, maximum peak load generally only occurs in summer, and the system operates at low partial load mode during the remaining time.Therefore, if the designed capacity is too large, the system consumes excessive power.The partial load rate (PLR) of the chiller can be expressed as (1): The power consumption of a chiller and its PLR share a certain relationship; according to the chiller power consumption equation of [8] as shown in (2) and (3): Equations ( 2) and (3) are chiller power consumption for Example 1 and Example 2 introduced from literate [8].In (2) and (3), the coefficients ai, bi, ci, and di define the relationship between power consumption and PLR for the chillers in Examples 1 and 2; this analysis was based on [8].The PLR ranged between 0.0 and 1.0.
The overall objective of OCL optimization was to find the optimal partial loading rate of each chiller of the multi-chiller system that satisfied the cooling requirements of the venue but also minimized overall power consumption.Therefore, the objective function of the proposed OCL optimization was defined as shown in (4): In (4), the parameter i signifies the ith chiller; n is the total number of chillers in the multi-chiller system; Pi is defined as the power consumption (kW) of the ith chiller.The object of Equation ( 4) is to find the lowest total power consumption in the multi-chiller system.
OCL problems have two types of restrictions [8] during solving; the first restriction is that the total output of RT must be equal to the RT required by the venue, as shown in (5).Qi signifies the rated RT capacity of the ith chiller; CL is the total RT required by the venue.This is the basic constraint for OCL problems.If the total RT generated by all chillers is smaller or larger than the RT requirement, the people who stay in the venue will feel not very comfortable.

𝑃𝐿𝑅 𝑄 𝐶𝐿
(5) Generally, the maximum capacity of a chiller is designed to meet the maximum peak load demand, but because of actual venue requirements and change of seasons, maximum peak load generally only occurs in summer, and the system operates at low partial load mode during the remaining time.Therefore, if the designed capacity is too large, the system consumes excessive power.The partial load rate (PLR) of the chiller can be expressed as (1): The power consumption of a chiller and its PLR share a certain relationship; according to the chiller power consumption equation of [8] as shown in (2) and (3): Equations ( 2) and (3) are chiller power consumption for Example 1 and Example 2 introduced from literate [8].In (2) and (3), the coefficients a i , b i , c i , and d i define the relationship between power consumption and PLR for the chillers in Examples 1 and 2; this analysis was based on [8].The PLR ranged between 0.0 and 1.0.
The overall objective of OCL optimization was to find the optimal partial loading rate of each chiller of the multi-chiller system that satisfied the cooling requirements of the venue but also minimized overall power consumption.Therefore, the objective function of the proposed OCL optimization was defined as shown in (4): In (4), the parameter i signifies the ith chiller; n is the total number of chillers in the multi-chiller system; P i is defined as the power consumption (kW) of the ith chiller.The object of Equation ( 4) is to find the lowest total power consumption in the multi-chiller system.
OCL problems have two types of restrictions [8] during solving; the first restriction is that the total output of RT must be equal to the RT required by the venue, as shown in (5).Q i signifies the rated RT capacity of the ith chiller; CL is the total RT required by the venue.This is the basic constraint for OCL problems.If the total RT generated by all chillers is smaller or larger than the RT requirement, the people who stay in the venue will feel not very comfortable.
The second restriction is that the partial loading of each chiller cannot be less than 30% [8], as defined in (6): PLR i ≥ 0.3 (6)

Two-Stage Differential Evolution Algorithm
The proposed two-stage DE algorithm is an optimization method that integrates MBDE and DE; it has superior exploration and convergence rates during searches for optimal solutions.

Differential Evolution Algorithm
The DE (differential evolution) algorithm was proposed by Price and Storn in 1995 [18].First, the initial real-number entities were randomly generated; the mutation operator used random numbers to select different entities to generate vector differences to serve as search directions, which were multiplied by a weight to obtain the size of the search, thus forming new search vectors.Then, through crossover, assessment, and selection, the process was iterated until the termination condition was met.The three main steps of the DE algorithm are: mutation, crossover, and selection.

Mutation Operator
The purpose of the DE mutation operator is to generate different mutated vectors.Those vectors are generated by (7) and (8), where G represents the current iteration and G+1 represents the next iteration. V V From the group that ranges from Entity 1 to Entity P, randomly select two entities from whole DE group, the vectors X r2,G and X r3,G , perform subtraction between these two entities to find the difference vector, then multiply it by a given weight F. Finally, add the global best vector X best,G or the self vector X i,G to obtain the next generation of mutated vectors V i,G+1

Crossover Operator
After an individual entity applies the mutation mechanism to generate a mutated vector, the crossover operator is used to generate trial vectors; the probability Cr of the crossover operator can range from 0.0 to 1.0.D represents the total design elements (or variables) in one vector.Regarding the randomly generated R, if the value of R is greater than Cr, then the jth element of original vector X ji,G is selected; if the value of R is smaller than Cr, then the jth element of mutated vector V ji,G+1 is selected as the jth element of trial vector U ji,G+1 for the next generation.The crossover operator is formulated as shown in (9):

Selection Operator
The selection operator is a convergence mechanism of the DE algorithm; it calculates the value of the objective function of the original vector X i,G and the trial vector U i,G+1 , and compares the advantages and disadvantages of the objective function.The vector with the higher (or better) objective Energies 2019, 12, 622 5 of 12 function value becomes the next generation's original vector X i,G+1 , then the vector with the highest objective function value of all original vector is selected as the global best vector X Best,G .

Modified Binary Differential Evolution (MBDE) Algorithm
A modified DE algorithm [15,16] was proposed by Wu and Tseng in 2010.The real number encoding of the original DE algorithm was changed to a binary encoding, and a modified mutation operator was used for evolution.The method was also applied in engineering optimization problems [19]; the mechanism of the MBDE was the same as that of the original DE.

Mutation Operator
The mutation equation for MBDE obtained by the original DE mutation operator [18] uses logical operators to compute binary strings as shown in (10) and (11).X i,G is original vector and X r1,G selects one entities randomly from the whole DE group.X best,G represents the global best vector.This binary mutation uses the logical operator XOR (exclusive or) to divide the binary string into two groups, "0" and "1" strings.The "0" string group represents common characteristics between two entities; the set weight F2 and randomly generated value are compared; when the randomly generated value is smaller than F2, then the bit "0" or "1" is mutated to "1" or "0".Alternatively, when the randomly generated value is greater than F2, then the bit remains unchanged.The "1" string group represents different characteristics between two entities; the given weight F1 and randomly generated value are compared; when the randomly generated value is smaller than F1, then the bit "0" or "1" is converted to "1" or "0"; alternately, when the randomly generated number is greater than F1, the bit remains unchanged.Finally, the two strings are combined to form the next generation of mutated vectors V i,G+1 ; the flowchart of the modified binary mutation operator is depicted in Figure 2.
Weight F1 tends to be greater than F2 because there is a higher probability that the characteristic of the optimal solution is a common characteristic, therefore a lower probability is used to maintain common characteristics, and a higher probability is used to change different characteristics.

Modified Binary Differential Evolution (MBDE) Algorithm
A modified DE algorithm [15,16] was proposed by Wu and Tseng in 2010.The real number encoding of the original DE algorithm was changed to a binary encoding, and a modified mutation operator was used for evolution.The method was also applied in engineering optimization problems [19]; the mechanism of the MBDE was the same as that of the original DE.

Mutation Operator
The mutation equation for MBDE obtained by the original DE mutation operator [18] uses logical operators to compute binary strings as shown in (10) and (11).Xi,G is original vector and Xr1,G selects one entities randomly from the whole DE group.Xbest,G represents the global best vector.This binary mutation uses the logical operator XOR (exclusive or) to divide the binary string into two groups, "0" and "1" strings.The "0" string group represents common characteristics between two entities; the set weight F2 and randomly generated value are compared; when the randomly generated value is smaller than F2, then the bit "0" or "1" is mutated to "1" or "0".Alternatively, when the randomly generated value is greater than F2, then the bit remains unchanged.The "1" string group represents different characteristics between two entities; the given weight F1 and randomly generated value are compared; when the randomly generated value is smaller than F1, then the bit "0" or "1" is converted to "1" or "0"; alternately, when the randomly generated number is greater than F1, the bit remains unchanged.Finally, the two strings are combined to form the next generation of mutated vectors Vi,G+1; the flowchart of the modified binary mutation operator is depicted in Figure 2.
Weight F1 tends to be greater than F2 because there is a higher probability that the characteristic of the optimal solution is a common characteristic, therefore a lower probability is used to maintain common characteristics, and a higher probability is used to change different characteristics.

Crossover Operator
The binary crossover operator of the modified binary DE algorithm is similar to the crossover operator of the DE algorithm.For the mutually equal jth design variable, a randomly generated value and the set crossover rate are compared; if the randomly generated value is greater than the crossover rate, then the original vector is selected as the next-generation trial vector.Alternately, if the randomly generated value is smaller than the crossover rate, then the mutated vector becomes the next-generation trial vector, as shown in Figure 3.

Crossover Operator
The binary crossover operator of the modified binary DE algorithm is similar to the crossover operator of the DE algorithm.For the mutually equal jth design variable, a randomly generated value and the set crossover rate are compared; if the randomly generated value is greater than the crossover rate, then the original vector is selected as the next-generation trial vector.Alternately, if the randomly generated value is smaller than the crossover rate, then the mutated vector becomes the next-generation trial vector, as shown in Figure 3.

Selection Operator
The operation of selection is similar to the assessment and selection of DE.After binary mutation and binary crossover, the objective function F(Xi,G) of the original vector is compared with the objective function F(Ui,G+1) of the trial vector, and the best vector is selected as the next-generation objective vector Xi,G+1.The vector with highest objective function value among all DE group is selected as the global best vector Xbest,G for mutation operator.The evolutionary computational system repeats Sections 3.2.1,3.2.2, and 3.2.3iteratively until the termination condition is met.

Two-Stage Differential Evolution Algorithm
The two-stage DE algorithm solves OCL problems in two stages.In stage 1, the solving of MBDE begins after the PLR of each chiller has been encoded to binary, restrictions have been considered, and the objective equation has been defined.The character of binary encoding has low numerical precision but has favorable diversity when searching for an optimal solution [17].After stage 1 has been solved, the resultant real-number encoding is taken into stage 2 of DE for exploitation.The overall architecture of the two-stage DE algorithm is illustrated in Figure 4.

Selection Operator
The operation of selection is similar to the assessment and selection of DE.After binary mutation and binary crossover, the objective function F(X i,G ) of the original vector is compared with the objective function F(U i,G+1 ) of the trial vector, and the best vector is selected as the next-generation objective vector X i,G+1 .The vector with highest objective function value among all DE group is selected as the global best vector X best,G for mutation operator.The evolutionary computational system repeats Sections 3.2.1-3.2.3 iteratively until the termination condition is met.

Two-Stage Differential Evolution Algorithm
The two-stage DE algorithm solves OCL problems in two stages.In stage 1, the solving of MBDE begins after the PLR of each chiller has been encoded to binary, restrictions have been considered, and the objective equation has been defined.The character of binary encoding has low numerical precision but has favorable diversity when searching for an optimal solution [17].After stage 1 has been solved, the resultant real-number encoding is taken into stage 2 of DE for exploitation.The overall architecture of the two-stage DE algorithm is illustrated in Figure 4.

Results, Analyses and Discussions
The OCL problems in [8] can serve as test examples for the proposed two-stage DE algorithm.Examples 1 and 2 are OCL problems for a six-chiller system and a four-chiller system.The parameters were set as follows: the population was 20; the mutation operator used (8) for DE and (11) for MBDE; the DE mutation rate and crossover rate were 0.5; the MBDE crossover rate was 0.5; the mutation rates F1 and F2 were 0.5 and 0.005, respectively.The aforementioned parameters were set in reference to the recommendations of [15] and [18].The total number of iterations was configured to fit the required number of calculations of the example.The number of iterations in of case study 1 was 1000; it was 40 in case study 2.

Case Study 1
This multi-chiller system was composed of six chillers.The objective was to calculate the lowest power consumption combinations of this six-chiller system under different RT requirements.The venue RT requirements in Example 1 were of five types: 6858 (90%); 6477 (85%); 6096 (80%); 5717 (75%); and 5334 (70%).Equation (2) defines the power consumption of this chiller system.The rated RT of each chiller and the power consumption parameters (ai, bi, ci) are listed in Table 1.The present study used a comparison method similar to those of other studies [8]; the total number of function evaluations was 20,000 for each run; 30 runs were executed, after which the mean, standard deviation (SD), maximum, and minimum values were computed in the manner of previously published articles in the literature.

Results, Analyses and Discussions
The OCL problems in [8] can serve as test examples for the proposed two-stage DE algorithm.Examples 1 and 2 are OCL problems for a six-chiller system and a four-chiller system.The parameters were set as follows: the population was 20; the mutation operator used (8) for DE and (11) for MBDE; the DE mutation rate and crossover rate were 0.5; the MBDE crossover rate was 0.5; the mutation rates F1 and F2 were 0.5 and 0.005, respectively.The aforementioned parameters were set in reference to the recommendations of [15,18].The total number of iterations was configured to fit the required number of calculations of the example.The number of iterations in of case study 1 was 1000; it was 40 in case study 2.

Case Study 1
This multi-chiller system was composed of six chillers.The objective was to calculate the lowest power consumption combinations of this six-chiller system under different RT requirements.The venue RT requirements in Example 1 were of five types: 6858 (90%); 6477 (85%); 6096 (80%); 5717 (75%); and 5334 (70%).Equation (2) defines the power consumption of this chiller system.The rated RT of each chiller and the power consumption parameters (a i , b i , c i ) are listed in Table 1.The present study used a comparison method similar to those of other studies [8]; the total number of function evaluations was 20,000 for each run; 30 runs were executed, after which the mean, standard deviation (SD), maximum, and minimum values were computed in the manner of previously published articles in the literature.
Table 2 summarizes the optimization results of the proposed method with DCSA (differential cuckoo search approach) [8], showing that for both the minimum value and SD, all values of the optimal solution found by the two-stage DE algorithm under different CL conditions that were very close to those of DCSA.When CL was 80%, 75%, and 70%, the proposed method executed 20,000 function evaluations (2000 function evaluations for phase 1 and 18,000 function evaluations for phase 2) to Energies 2019, 12, 622 8 of 12 achieve better results than DCSA.Table 3 compared the optimal solutions of two-stage DE algorithm with those of other methods; when CL was 90% and 80%, a solution similar to that of DCSA was found; the solution was superior to than that of SA [20] and PSO [6].Under other CL conditions (85%, 75%, and 70%), the results of the proposed two-stage DE algorithm were superior to those of other studies [6,8,20].The convergence of two-DE algorithm for case study 1 is shown in Figure 5. phase 2) to achieve better results than DCSA.Table 3 compared the optimal solutions of two-stage DE algorithm with those of other methods; when CL was 90% and 80%, a solution similar to that of DCSA was found; the solution was superior to than that of SA [20] and PSO [6].Under other CL conditions (85%, 75%, and 70%), the results of the proposed two-stage DE algorithm were superior to those of other studies [6,8,20].The convergence of two-DE algorithm for case study 1 is shown in Figure 5.

Case Study 2
This multi-chiller system was composed of four chillers.The objective of this example was the same as that of Example 1, namely, to find the lowest power consumption combination for this four-chiller system at different RT requirements.The RT requirements of the venue were of six types: 2610 (90%); 2320 (80%); 2030 (70%); 1740 (60%); 1450 (50%); and 1160 (40%).Equation (3) defines the power consumption of the chiller, the rated RT of each chiller and their power consumption parameters (a i , b i , c i , d i ) are listed in Table 4.The study referenced the comparison methods of relevant articles [4,[6][7][8], and 800 function evaluations were allowed in each run, and 30 runs were executed.The mean, SD, maximum, and minimum values were calculated for comparison with those of other studies.Table 5 proves that in addition to the power consumption of the proposed two-stage DE algorithm being relatively similar to that of DCSA at CL = 80%, under other CL conditions (70%, 60%, 50%, and 40%), the results of the proposed optimization method using 800 function evaluations (80 function evaluations for phase 1 and 720 function evaluations for phase 2) were also superior to that of DCSA [8].
In particular, at CL = 50% and CL = 40%, the SD results validated that the proposed two-stage DE algorithm was more stable than DCSA.The results of comparisons with other methods [4,[6][7][8] are listed in Table 6, they prove that the proposed two-stage DE algorithm could find a power consumption combination that was similar to or lower than those of other methods (GA [4], PSO [6], DE [7], CSA [8]) under most CL conditions.The convergence of two-DE algorithm for case study 2 is shown in Figure 6.

Conclusions
The proposed two-stage DE algorithm integrates the characteristics of binary and real-valued DE algorithms; it has excellent exploration and convergence rates for optimal solutions.Regarding case studies 1 and 2, the comparisons of the proposed method with other methods indicate that the

Conclusions
The proposed two-stage DE algorithm integrates the characteristics of binary and real-valued DE algorithms; it has excellent exploration and convergence rates for optimal solutions.Regarding case studies 1 and 2, the comparisons of the proposed method with other methods indicate that the proposed two-stage DE algorithm is suitable for optimization of the power consumption configurations of multi-chiller systems.In addition to being able to obtain solutions similar to or better than those of referenced studies, the degree of variance of the optimal solution in each search was better than those of other studies, further validating that the performance and stability of the proposed two-stage DE algorithm are better than those of other methods.The proposed algorithm can be used to solve similar optimization problems.

Figure 5 .
Figure 5. Convergence of two-stage DE for case study 1.

Figure 5 .
Figure 5. Convergence of two-stage DE for case study 1.

Figure 6 .
Figure 6.Convergence of two-stage DE for case study 1.

Figure 6 .
Figure 6.Convergence of two-stage DE for case study 1.

Table 1 .
Power consumption coefficient and rated RT information for case study 1.

Table 3 .
Case study 1: comparison with methods proposed by other studies.

Table 4 .
Case study 2: Power consumption coefficient and rated RT information.

Table 6 .
Case study 2: comparison with the methods proposed by other studies.