Optimum design of the booster chlorination systems by using hybrid HS-Solver optimization approach Hibrit HS-Solver optimizasyon yaklaşımı kullanılarak ek klorlama istasyonlarının optimum tasarımı

Öz A hybrid simulation-optimization approach is developed in this study for optimally designing the booster chlorination systems in water distribution networks. In the developed approach, chlorine residuals in the demand locations are determined by using the response matrix (RM) approach. The generated RM is then integrated to an optimization model where a hybrid HS–Solver optimization approach is used. HS-Solver is a recently proposed hybrid approach which integrates the harmony search (HS) algorithm and a spreadsheet Solver as the global and local optimizers, respectively. The objective of the HS-Solver in the developed approach is for optimally designing the booster stations by maintaining the chlorine residuals within desired limits. The applicability of the developed approach is tested on a real water distribution network. Identified results indicate that the developed HS-Solver based solution approach determined similar or better results than those obtained by using the different solution approaches in literature. Bu çalışmada, su dağıtım şebekelerindeki ek klorlama istasyonlarının optimum tasarımı amacıyla bir hibrit simülasyon-optimizasyon yaklaşımı geliştirilmiştir. Geliştirilen yaklaşımda, talep noktalarındaki bakiye klor konsantrasyonları tepki matrisi yaklaşımı ile belirlenmiştir. Oluşturulan tepki matrisi ardından hibrit HS-Solver optimizasyon yaklaşımının kullanıldığı bir optimizasyon modeline entegre edilmiştir. HS-Solver, armoni araştırması ile elektronik tablolama programlarında kullanılan Solver eklentisinin sırasıyla global ve lokal optimizasyon aracı olarak kullanıldığı yeni önerilmiş bir optimizasyon yaklaşımıdır. Geliştirilen yaklaşımda HS-Solver’ın amacı bakiye klor konsantrasyonlarını arzu edilen limitlerde tutacak şekilde ek klorlama istasyon sisteminin optimum tasarımını yapmaktır. Geliştirilen yaklaşımın uygulanabilirliği gerçek bir su dağıtım şebekesi üzerinde test edilmiştir. Elde edilen sonuçlar, geliştirilen HS-Solver tabanlı çözüm yaklaşımı ile literatürde verilenlerle uyumlu veya daha iyi sonuçlar elde edilebildiğini göstermiştir.


Introduction
Disinfection of water within the water distribution systems is an important problem to protect the public health. A conventional way of disinfection is to inject a large quantity of chlorine from the outlet of the water purification plants to maintain chlorine residuals at the critical locations of the network. Although this is a commonly considered application procedure, it may create excessive chlorine residuals in the vicinity of the injection point. That cannot only cause the taste and odor problems, but it can also result in the formation of carcinogenic disinfection-by-products (DBPs) [1]. On the other hand, reduced chlorine injections at the injection point may not be sufficient to kill bacteria at the remote points of the network [2]. These problems can be alleviated by installing booster chlorination stations in the critical locations where low chlorine residuals are observed.
The advantage of booster chlorination stations is the ease in maintaining the chlorine residuals within desired limits throughout the network. More specifically, the following objectives can be defined for utilizing the booster chlorination stations [3]: i) To minimize the total disinfectant dose, ii) To minimize DBPs, iii) To minimize the investment and operation costs; iv) To maximize the volume of water supplied to consumers within desired concentration limits.
Regarding these objectives, the following questions can be asked: How many booster chlorination stations should be installed and where on the network? What should be the best chlorine injection pattern for maintaining the chlorine residuals within desired limits in the network? Answers to these questions should be evaluated by decision makers when installing the best booster chlorination station network. Therefore, identification of the numbers, locations, and chlorine injection schedule of the booster stations becomes a challenging engineering optimization problem.
The current literature includes various solution approaches to solve the booster station optimization problems by means of the deterministic and heuristic solution approaches. Deterministic approaches include linear programming (LP) [4], mixed integer LP (MILP) [5], quadratic programming (QP) [1], mixed integer QP (MIQP) [6], etc. Similarly, heuristic optimization approaches include genetic algorithms (GA) [2], [7], [8], genetic immune algorithm (GIA) [9], ant colony optimization (ACO) [10], particle swarm optimization (PSO) [11], harmony search (HS) [12], and differential evolution (DE) [13], etc. Among the deterministic approaches, Bocelli et al. [4] first described the problem of booster station optimization by developing an LP formulation. The objective of their LP formulation is defined as the minimization of the chlorine injection dosages at pre-defined injection points by maintaining the chlorine residuals in all demand locations and measurement times. One of the most important outcomes of their study is that chlorine residuals at demand locations linearly changes with the injected chlorine concentrations in case of the first-order wall and bulk reaction kinetics. By using this relationship, the chlorine residuals at demand locations can be calculated by means of the response matrix (RM) approach. As an extension, Tryby et al. [5] used an RM approach to determine both locations and chlorine injection dosages of the booster stations by using an MILP based solution approach. Propato and Uber [1] developed a linear least square (LLS) problem to determine the chlorine injection dosages by means of a QP model. Utilizing the same LLS formulation, Propato and Uber [6] determined the locations and the chlorine injection dosages of the booster stations by using the MIQP based solution approach.
Heuristic optimization approaches were also applied to the solution of booster station optimization problems. These approaches can be used to solve linear or nonlinear optimization problems without taking the special initial solutions to start the search process. Also, they do not require derivative information of the objective function to find a search direction. Therefore, they can handle the optimization problems with both continuous and discrete decision variables.
There are many studies in literature which considered heuristic optimization approaches to solve the booster station optimization problems. Munavalli and Mohan Kumar [2] applied a GA to determine the disinfectant dosages at the pre-defined booster locations. Prasad et al. [3] applied a multiobjective GA model to find the booster locations and the corresponding injection dosages. Ozdemir and Ucaner [8] applied a GA model for determining the booster locations as well as the input concentrations. Ostfelt and Solomons [7] applied a GA model to solve two optimization problems where the first deals with the design/operation of the booster station to satisfy the minimum cost condition, and the second aims to maximize the health protection. Chu et al. [9] developed a solution methodology which applied a GIA-based optimization approach for solving the optimal scheduling problems of the booster disinfection station. Wang and Guo [10] employed an ACO-based optimization model to find the location and injection dosages of the booster stations on a hypothetical water distribution network. Gokce and Ayvaz [12] developed a solution methodology which considered HS to solve the booster station scheduling problems. Similarly, Gokce and Ayvaz [13] applied DE to find the locations as well as the chlorine injection concentrations of the booster stations.
It should be noted that although heuristic optimization approaches can be used to obtain efficient results in terms of the locations as well as the chlorine injection dosages of the booster stations, they usually require a high number of model simulations to find a feasible solution. In order to handle this difficulty, implementation of the hybrid optimization approaches is widely preferred to complex optimization problems. Note that hybrid optimization approaches consist of the mutual integration of the heuristic and deterministic optimization approaches. In this integration, global exploration process starts using the heuristic approach with multiple starting points, and then, a deterministic approach finds the optimal solution by getting the best of the multiple solutions as the starting point [14]. This kind of a solution methodology makes obtaining the global optimum solution easier than both heuristic and deterministic optimization approaches by themselves [15]. Regarding this hybridization scheme, Ayvaz and Kentel [16] integrated a binary GA with a LP based optimization model to determine numbers, locations and chlorine injection dosages of the booster stations. In this integration, a GA is used to determine the numbers and the locations of the booster stations, and LP is used to determine the corresponding chlorine concentrations for each booster location. Although their proposed GA-LP approach can successively solve the booster station optimization problems, implementation of these two algorithms can require advanced programming skills especially for the researchers who are not familiar with the complex coding/decoding tasks in GA.
The main objective of this study is to propose a new simulationoptimization methodology to solve the booster station optimization problems in water distribution networks. In the simulation part of the proposed methodology, water quality response of the given network to a specific chlorination schedule is calculated based on the RM approach which is proposed by Bocelli et al. [4]. In order to use this approach, the required response coefficients are calculated by modeling the given water distribution network on EPANET model [17]. The developed simulation part is then integrated to an optimization model where a hybrid HS-Solver optimization approach is used. HS-Solver is a recently proposed hybrid optimization approach which both integrates the heuristic harmony search algorithm and a spreadsheet Solver add-in as the global and local optimizers, respectively. The objective of the HS-Solver based solution methodology is to determine the locations as well as the chlorine concentrations of the booster stations by maintaining the chlorine residuals within the desired limits at all the demand locations and measurement times. The performance of the proposed methodology is evaluated by solving an existing water distribution network on both HS and HS-Solver. Comparison of the identification results indicated that the proposed HS-Solver based simulation-optimization methodology finds similar or better results than those obtained by using the stand-alone HS and different optimization approaches in literature.

Problem formulation
Locations and chlorine injection dosages of the booster stations are determined by using an optimization model. The primary goal of the proposed model is defined as the minimization of the injected chlorine mass to the system by maintaining chlorine residual limits for all demand locations and measurement periods. This problem can be mathematically defined as follows: Let be the number of demand locations where chlorine residuals are monitored, ℎ be the number of observation time periods, be the starting time of the observation, be the volumetric demand [ 3 ] satisfying the chlorine residual limits at location in observation period , be the total volume of water demand [ 3 ] for a given hydraulic period, be the demand be the injected chlorine dosage [ −3 ] from booster station at injection period , and ̃ be the total outflow [ 3 −1 ] at location . Given these definitions, the objective function and the constraints can be formulated as follows: Where is the objective function which represents the chlorine injection mass rate [ −1 ]. The constraint in Equation (2) is used to maintain the non-negativity condition of the injection dosages. Similarly, Equation (3) is used to control the chlorine residuals within the specified lower and upper bounds. For this purpose, the value of should be calculated for all the demand locations and measurement times. Note that this task is usually conducted by directly linking the EPANET model to the optimization approaches to calculate the value of for the generated chlorination schedule. Although this is an effective approach, running the EPANET model for each optimization cycle may require long computational times for the networks having long simulation times and big dimensions. Therefore, the RM approach is used to calculate the chlorine residuals at demand locations. According to Bocelli et al. [4], the value of can be calculated as follows: Where is the composite response coefficient which is determined by using = ⁄ . Note that values of are computed by means of the EPANET. For this purpose, EPANET model is separately executed for each potential booster location by assigning the unit chlorine injections over there.
As indicated previously, the optimization formulation given in Equations (1) to (3) is solved by using the hybrid HS-Solver optimization approach. Since HS is an unconstrained heuristic optimization approach, violation of the constraint in Equation (3) is evaluated by means of the penalty function approach. Based on this approach, the objective function in Equation (1) is transformed to the penalized form as follows: Where ( min ) and ( max ) are the penalty functions for maintaining lower and upper chlorine residual limits, respectively and 1 and 2 are the associated penalty coefficients which are used to adjust the magnitudes of the ( min ) and ( max ). Note that several trial runs are conducted before executing the proposed approach. According to the results of these runs, values of 1 and 2 are taken as 10.000 and 1.000, respectively. Although the penalty functions given above are used to prevent violation of the constraints during the HS solution, they are not required for the Solver since it can solve the constrained optimization problems through its builtin constraint module.

Optimization model
As indicated previously, the optimization problem given above is solved by using the hybrid HS-Solver optimization approach. HS-Solver, proposed by Ayvaz et al. [14], is a powerful and user-friendly optimization approach which integrates the global exploration feature of the HS and strong local search capability of a spreadsheet Solver. Using this integration, various unconstrained, constrained, and structural optimization problems were previously solved in Ayvaz et al. [14]. Also, Ayvaz and Elçi [18], [19] applied HS-Solver to the solutions of two groundwater management problems; one is related to pumping maximization and the other was developed to solve pumping cost minimization problems on the same watershed. Note that HS is a heuristic optimization algorithm inspired from the musical improvisation process. In musical processes, a pleasing harmony can be obtained by following three rules: i) Playing a note randomly, ii) Playing a note from harmony memory, iii) Playing a note which is close to another one stored in harmony memory.
Geem et al [20] first adapted these musical rules to solve engineering optimization problems as: Combinations of these rules in an optimization framework allow obtaining a global optimum solution in HS optimization algorithm because these rules together make a stochastic derivative which shows a direction to the global optimum solution [21]. This stochastic derivative, which is based on collective intelligence of musicians' experiences, is even applicable to discrete variables that cannot have derivative information from traditional differential calculus approach.
Nowadays, electronic spreadsheets have become necessary tools to perform various engineering computations. Almost all the commercially available spreadsheet products such as Excel® include a built-in "Solver" module to solve unconstrained and constrained optimization problems [22]. Solver can solve the linear and nonlinear optimization problems by means of the Simplex and Generalized Reduced Gradient methods, respectively. Since Solver works on Excel®, all its computational features are accessible from the Visual Basic for Applications (VBA) platform. Therefore, HS and Solver processes have been hybridized on VBA platform by developing three separate VBA modules. The first module is for the stand alone HS and can be directly used to solve any optimization problem. The second module is developed for calling the Solver, which is created by using the macro recording feature of Excel®. The last module is used to integrate the HS and Solver modules to generate a hybrid optimization approach. Note that there are two options to integrate the HS and Solver processes in the last module. In the first option, the entire search space is explored by HS, then Solver gets the best result of HS as a starting point to precisely find the global optimum solution. In the second option, both HS and Solver run simultaneously such that all the solutions of HS are subjected to local search by Solver based on a probability of [14]. For the second option, previous experiences state that use of a small probability is sufficient for solving many engineering optimization problems. Therefore, the second option is considered by setting the probability of = 0.10 for all the solutions. The flowchart of the hybrid HS-Solver optimization approach is given in Figure 1. As can be seen, the required solution parameters are given in the first step which are the harmony memory size (HMS), harmony memory considering rate (HMCR), pitch adjusting rate (PAR), fret width (fw), and the probability of . A detailed description of this parameters and HS can be found in [14] and [20].

Numerical application
In order to evaluate the performance of the proposed simulation-optimization approach, water distribution network of the Cherry Hill-Brushy Plains, CT, USA is used. This network is previously used in different studies in literature to solve the booster station optimization problems [1], [2], [4]- [6], [12], [13], [16]. Figure 2 shows the layout of the water distribution network under consideration. As can be seen from Figure 2, the network consists of a pumping station at the 1 st junction, a storage tank at the 26 th junction, 34 demand locations from 2 nd to 25 th and 27 th to 36 th junctions, and 47 pipes which have the total length of 11.26 km. The network also includes 6 hypothetical junctions to locate booster stations which are also shown in Figure 2 as junctions A-F. The reason of considering these hypothetical junctions in the previous studies was to inject the disinfectant without considering any water demand at the same points. Note that these junctions are also considered in this study in order to obtain comparable results with those obtained in literature. Figure 3 Figure 3(c). Using these data, the network's hydraulic behavior which is given in Figure 4 is determined by modeling the network on EPANET. It is clearly seen that the hydraulic behavior of the network is controlled by the pumping station in the 1 st and the 3 rd six hour time period of a day. For the 2 nd and the 4 th six hour time periods, the tank is used since the pump does not work at those times. This situation results with the change of the flow directions in the network for six hour time periods.
Since the chlorine residuals in the demand locations are determined by utilizing the RM approach, it is required to execute the EPANET model to determine the composite response coefficients of before starting the search process. In these model runs, a unit chlorine concentration is injected from each potential booster location and the corresponding chlorine residuals are collected at demand locations to compute values [23]. Note that this kind of a solution approach requires a periodically changed hydraulic behavior which is observed when chlorine residuals at two successive days are exactly same. Satisfaction of this condition can be possible in cases of long simulation times. Thus, simulation time is taken as 288 hours since the periodically changed hydraulic behavior is observed after 264 th hour in the system. After this process, the coefficients of are calculated by using the chlorine residuals at the last 24 hour of the 288 hour time period. The bulk and wall decay coefficients are assumed to be 0.50 1/d and 0 in these computations.  As can be seen from Figure 5, the optimization process is performed by using the same random number seeds that is why both HS and HS-Solver start from the same initial solutions. When the identified results of HS and HS-Solver are compared, it can be seen that HS-Solver requires less improvisations than HS to minimize the total mass. As an example, for = 1, improvement of the objective function value remains constant after about 2.109 th improvisation in HS. However, the same objective function value is obtained after about the 89 th improvisation in HS-Solver. When the results for different booster station numbers are compared, it is shown that inclusion of the additional booster stations significantly improves the objective function values as an expected behavior. For a different number of booster stations, Table 1 compares the results of HS and HS-Solver in terms of the identified booster locations and injected chlorine concentrations. As can be seen, for = 1, both HS and HS-Solver identified the 2 nd junction as a booster location which is located in just downstream of the pumping station. For this location, both models determined the same chlorine concentrations of 1.780 mg/L. For = 2, both models placed the second booster station next to tank in order to supply chlorine to the network in the second and the fourth 6-hour time period of a day. For other solutions, both HS and HS-Solver determined the same or nearly the same booster stations having the similar injection patterns. Table 2 compares the final identified chlorine injection mass rates for HS and HS-Solver with those obtained by using MIQP and DE solution approaches. As can be seen, for = 1 HS and HS-Solver obtained the same chlorine mass rate (3.010 g/day). For the remaining solutions, although the results are very close to each other, HS-Solver provides slightly lower chlorine mass rates than HS. When the results of the HS-Solver are compared with those obtained by using DE and MIQP, it can be seen that HS-Solver provides identical results with DE and better results than MIQP. For this table, another comparison is conducted in terms of the computation loads of the proposed approaches. As can be seen, the required number of improvisations in HS-Solver is lower than HS. For example, while HS required 2,109 improvisations to find a chlorine mass rate of 3.010 g/day for = 1, the same result is obtained after 89 improvisations in HS-Solver. Same situations are also observed in other solutions such that HS-Solver provided better results in fewer improvisations than the stand-alone HS. Note that since HS-Solver resulted with the identical chlorine mass rates with DE, comparison of their computation loads is a crucial step. Table 2 lists the required number of generations in DE to find the given chlorine mass rates. As can be seen, for = 1, DE required 47 generations which corresponds to 940 simulation runs whereas the identical result is obtained in HS-Solver in 89 improvisations. This difference significantly increases as the complexity of the problem increases.
For both HS and HS-Solver, statistical evaluation of the chlorine residuals in terms of the mean, minimum and maximum values is given in Table 3.   It is clearly seen that the mean of the chlorine residuals is in a decreasing trend with the number of booster stations. Minimum of the chlorine residuals is determined to be equal to the minimum residual limit of 20 mg/L for all the solutions. Note that when the maximum chlorine residuals for both HS and HS-Solver are compared with the injected chlorine concentrations in Table1, it can be seen that almost all the maximum residuals are greater than the identified values in Table 1. This situation is related to the network's hydraulic behavior such that the flow directions in the network change in 6-hour time periods of a day. Regarding this issue, Figure 6 compares the flow directions in the first and the second 6-hour time periods. As can be seen from Figure 6(a), the chlorinated water by the flow paced booster station in the 2 nd junction moves from 2 nd to the 3 rd and the 5 th junctions. Since the flow direction changes in the second 6 hour time period as given in Figure 6(b), the chlorinated water returns from 5 th to 2 nd junction. Since the booster station in the 2 nd junction works continuously through 24 hour time period, there is an additional chlorine injection to the previously chlorinated water that is why maximum chlorine residuals are obtained higher than the injected ones.
Comparisons of the identified booster locations by using HS and HS-Solver and those obtained by using DE and MIQP are given in Figures 7. As can be seen, the identified locations of HS and HS-Solver are all the same except for = 6 where only one booster station is located to a different demand location. The identified locations for HS-Solver and DE are all the same as a result of the chlorine injection patterns given in Table 2. When the model results for MIQP is examined in detail, it can be seen that the MIQP located the 1 st booster station to the junction A whereas the 2 nd junction is selected for all the solutions in HS-Solver. The remaining booster stations are located to the same junctions in = 2 to 4 and close locations in = 5 and 6.

Discussions and conclusions
A hybrid simulation-optimization approach is proposed to solve the booster station optimization problems in water distribution networks. In the proposed approach, chlorine residuals in the demand locations are determined by means of the RM approach in the simulation part. The locations as well as the chlorine injection dosages of the booster stations are determined by integrating the RM based simulation part to an optimization model where the hybrid HS-Solver optimization approach is used. The objective of the HS-Solver is to minimize the total injected chlorine mass rates by maintaining the chlorine residual limits for all the consumer nodes and measurement times. The applicability of the proposed approach is evaluated by using an existing water distribution network for different booster station numbers. Identified results indicated the proposed approach does not only efficiently determine the locations and chlorine injection dosages of the booster stations in fewer improvisations, but it also provides identical or slightly better results than those obtained by using different approaches found in literature. The following conclusions and discussions can be drawn from the proposed approach: The main advantage of using the hybrid HS-Solver optimization approach is its easy computational structure compared to other hybrid algorithms in literature. After programming the HS on VBA platform of the MS Excel® spreadsheet, a user can easily integrate the HS and Solver processes to create a hybrid optimization approach without requiring any advanced programming skills. However, this kind of an integration may require high computation times especially for the big problems since default Solver add-in of Excel® has some limitations. However, these limitations can be handled if premium or professional versions of the Solver add-in are used since they can solve the optimization problems 100 times faster than default Solver add-in (Frontline Systems, 2016). Using this improved version of HS-Solver, various engineering optimization problems can be solved without big computation time limitations.
In the proposed approach, chlorine residuals at the demand locations of the network are calculated by means of the RM approach. Although use of this approach significantly reduces the required computation time, it can only be used for the problems where first-order bulk and wall reactions kinetics are observed. However, if these conditions are not satisfied in the system, a water quality simulation model such as EPANET should be directly integrated with the HS-Solver to calculate the chlorine residuals at the demand locations. Such an integration may require high computation times especially for the big networks and/or long simulation times. For such cases, LLS based solution model given in [1] and [6] may be used.
In this study, a single chlorine injection period is considered for all the simulations. Also, all the calculations are conducted by taking the flow paced booster (FPB) into account. Therefore, solution of the problem by considering different booster types and multiple injection periods should be considered as a future research.