River water quality management using a fuzzy optimization model and the NSFWQI Index

In this study, a novel multiple-pollutant waste load allocation (WLA) model for a river system is presented based on the National Sanitation Foundation Water Quality Index (NSFWQI). This study aims to determine the value of the quality index as the objective function integrated into the fuzzy set theory so that it could decrease the uncertainties associated with water quality goals as well as specify the river’s water quality status rapidly. The simulation-optimization (S-O) approach is used for solving the proposed model. The QUAL2K model is used for simulating water quality in different parts of the river system and ant colony optimization (ACO) algorithm is applied as an optimizer of the model. The model performance was examined on a hypothetical river system with a length of 30 km and 17 checkpoints. The results show that for a given number of both the simulator model runs and the artificial ants, the maximum objective function will be obtained when the regulatory parameter of the ACO algorithm (i.e., q 0 ) is considered equal to 0.6 and 0.7 (instead of 0.8 and 0.9). Also, the results do not depend on the exponent of the membership function (i.e., γ ). Furthermore, the proposed methodology can find optimum solutions in a shorter time.


INTRODUCTION
Many rivers around the world carry significant amounts of pollutants produced by domestic, agricultural, and industrial activities (Farzadkhoo et al. 2018(Farzadkhoo et al. , 2019Hamidifar et al. 2015). To date, many waste load allocation (WLA) models have been proposed for water quality management in a river system. In general, the purpose of these models is to obtain technically and economically feasible solutions for the pollution control agency (PCA) and dischargers. Organizations set standards and regulations for water use, whereby dischargers are required to remove a percentage of various wastewater pollutants. On the other hand, increasing the fractional removal level of pollutants increases the cost of wastewater treatment. Therefore, trade-offs between conflicting goals between organizations and dischargers have always been a challenge. Many studies have been conducted in the field of water quality management in river systems, each of which has specific characteristics and goals. Mahjouri and Abbasi (2015) attempted to minimize the total cost of wastewater treatment by considering the optimal treatment alternatives and the objectives of the stakeholders. Nikoo et al. (2016) developed a WLA model for a river and applied the policy of allocating total profits from the permitted trade of wastewater discharge among dischargers. Tavakoli et al. (2015) studied the effects of WLA and wastewater load on river water quality management concerning agricultural water return. Estalaki et al. (2015) proposed appropriate penalty functions to comply with the standards imposed by PCA's water quality standard to dischargers. Ashtiani et al. (2015) and Saberi and Niksokhan (2017) optimized the degree of deviation or deviation from water quality standards in different parts of the river, which is one of the water quality management policies in rivers.
Fuzzy rules and concepts are used to reduce uncertainty in WLA models in the river system. There are various techniques in the fuzzy environment for solving water quality problems. Sasikumar and Mujumdar (1998) mathematically quantified the uncertainties related to vagueness in describing the qualitative goals and reducing water pollution and expressed it in a fuzzy framework, and so used the concepts of fuzzy set theory and its relations to the real sets. Karmakar and Mujumdar (2006) used the concept of gray fuzzy to develop their proposed model for water quality management. Qin et al. (2009) proposed an interval-parameter WLA model for river water quality management under uncertainty. Nikoo et al. (2016) simulated a pollutant charge allocation model and used a fuzzy interval solution. Xu et al. (2017) applied the fuzzy random environment approach to assess river pollution under uncertainty conditions. In some studies, the determination of the concentration of water quality parameters (WQPs) (e.g., either BOD or DO concentration or both) has been considered as an indicator of water quality assessment in the river (Takyi and Lence, 1996;Saadatpour and Afshar, 2007;Lei et al., 2015). However, in others, more than two quality parameters have been investigated (Mahjouri and Abbasi, 2015). Furthermore, in some models, some effective quality parameters are expressed as a separate and standard quality index and are considered as a criterion for water quality assessment (Kannel et al., 2007).
In this study, the National Sanitation Foundation Water Quality Index (NSFWQI) is used to determine the water quality in a river system. This popular index determines water quality using 9 effective quality parameters of environmental water. To reduce some uncertainties in the proposed model, the fuzzy set theory has been applied expansively. The objectives of this model are to consider the optimal water quality in the river for the stakeholders so that the costs of the wastewater treatment plants remain acceptable. Another new approach applied here is to optimize the fuzzy objective function in the form of the average NSFWQI quality index at the river checkpoints. The advantage of using the NSFWQI as the membership function in the proposed model is that it considers the uncertainties related to the WLA model and also can generate a set of more accurate solutions.
To simulate and optimize the river water quality model, the QUAL2K simulator and ant colony optimization (ACO) algorithm are used simultaneously to determine the appropriate solutions of the model, and thus to increase the speed and accuracy in the determination of model solutions. The model is run on a hypothetical river and the results are analysed. The model results can contribute to faster and more accurate decisions.

Water quality index (WQI)
One of the known indices that are used for determining water quality is the NSFWQI, which was provided by Brown et al. (1970) and was supported by the National Sanitation Foundation of the United States. For this reason, Brown's index is also referred to as NSFWQI. This index is scaled from 0 to 100 and is considered as a decreasing index because its values reduce with increasing water pollution. Nine WQPs that have a more important role in human health, with different weighting coefficients, are used in calculating the value of this index. The WQPs include DO, faecal coliforms (FC), BOD, pH, nitrates, phosphates, temperature, turbidity, and total solids (TS). For each of these WQPs, several curves have been presented in terms of the amount or concentration of the parameter and index value associated with the same parameter. For example, the DO concentration curve in terms of the index value is shown in Fig. 1.
After specifying the index value of each of the WQPs using such curves, the following equation can be used for obtaining the NSFWQI value: where NSFWQI is the index value, WQI i is the index value related to quality parameter i, w i is the weighting coefficient related to quality parameter i, and n is the number of WQPs (here it has 9 parameters). Then, based on the obtained NSFWQI value, the water quality states are obtained. The water quality state can be very poor (the index value of between 0 and 25), poor (25-30), medium (50-70), good (70-90), and very good (90-100) (Brown et al., 1970;Abbasi and Abbasi, 2012;Tyagi et al., 2013;Nada et al., 2016).

ACO algorithm
In recent decades, several metaheuristic algorithms have been developed that are generally derived from the activities of nature. These algorithms, including genetic algorithm and ACO, are highly accepted and have been widely used and developed (Dorigo and Stützle, 2019;Babaee Tirkolaee et al., 2020;Gholizadeh and Fazlollahtabar, 2020;Nayak et al., 2020). Also, these algorithms are used to solve linear and nonlinear problems and can obtain desirable or optimal solutions to the problem in less time. Inspired by the social behaviour of ants while searching for food through the release or evaporation of pheromone, the ACO can search for shorter paths (Afshar et al. 2009;Kalhor et al. 2011;Afshar et al. 2015). Briefly, for solving the proposed model using the ACO algorithm, the two following steps can be applied: Step 1. In each of the algorithm iterations, each of the ants selects a random path. Therefore, the collection of ants produces different solutions by passing from different paths. The probability of path choice by an ant is equal to: where P ij (t) is the probability of path selection i to j by each of the ants in iteration t, τ ij (t) is the heuristic coefficient of path i to j in iteration t, the exponents α and β parameters specify the importance of pheromone value and heuristic value with respect to each other and S is a set including the nodes of the graph of ants' paths that have the possibility of its selection (l is a component of this set). Suppose q is a random variable that is selected using the normal distribution in the closed interval [0,1] and q 0 (q 0 [0,1]) is a regulatory parameter. Each ant behaves in the following way for selecting the next node in its path (j): where J is a random variable that is selected according to the probability distribution of P ij (t).
Step 2. The solutions obtained in the previous step are evaluated based on the objective function of the problem. Then the updating process of the pheromone trail is done based on the path of the best ant as follows: where τ ij (t) is pheromone concentration deposited in the path i to the path j in iteration t, ρ (ρ[0,1]) is evaporation value of pheromone, the symbol iteration m shows the next iteration and Δτ ij is the updating value of pheromone trail in the path i to the path j. Δτ ij value is obtained by the following equation:
In fuzzy logic, there is a need for a set whose range is flexible and roughly defined. This type of set is called a fuzzy set. The membership of elements in fuzzy sets is expressed by the degree of membership, which is a number between 0 and 1. Also, the function that expresses the degree of membership of a set is called the membership function (Sasikumar and Mujumdar, 1998;Wong and Lai;Zimmermann, 2011).
As noted before, the fuzzy set theory was used for the proposed model. Water quality standards according to PCA's goals were considered based on the average WQI value. The WQI which was applied here (i.e. NSFWQI) was a decreasing index (value decreased with increasing pollution). Also, the water quality standards are usually restricted by bounds (into an interval). Suppose these standards are a fuzzy set. Based on the definition, a membership function of each fuzzy set has to select a continuous and real value between 0 and 1. In accordance with the abovementioned issues, the appropriate membership function (linear or nonlinear) for average WQI is given as follows: where µ(I i ) is the membership function of fuzzy set, I i is average NSFWQI for all checkpoints in the river system for ant i, I imin and I imax are the minimum and maximum of I i , respectively, and γ is the exponent of membership function and non-zero real number. This membership function is plotted in Fig. 2.
To calculate the optimum value of average NSFWQI in the river system, the following model was proposed: Subject to: where I ij is the average NSFWQI value for the ant i in the checkpoint j (the same NSFWQI value in Eq. 1), I ijmin and I ijmax are the minimum and maximum values of I ij , respectively, C k is total wastewater treatment costs for all treatment units in the river system and scenario k, C kl is wastewater treatment cost for scenario k and treatment unit l, C kmin and C kmax are the minimum and maximum C kl values, respectively, NCP is the number of the checkpoints in the river system and ND is the number of treatment units. Other parameters were defined above.
Equation 7 is the objective function that was considered in the proposed model. Equations 8 and 9 are crisp constraints related to the PCA's goal and limit the average amount of NSFWQI and NSFWQI values in each river checkpoint, respectively. Equations 10 and 11 are crisp constraints related to the discharger's goal and limit the total wastewater treatment cost and the amount of treatment cost in each treatment unit, respectively. Equations 12 and 13 address the method of calculating the average amount of NSFWQI and the total treatment cost, respectively (Sasikumar and Mujumdar, 1998;Zimmermann, 2011;Cullis et al., 2018;De Ketele et al., 2018;Namugize and Jewitt, 2018).

Simulation-optimization (S-O) approach
The S-O approach is one of the most practical techniques for solving the quality and quantity problems that determine the required problem parameters with the use of a suitable simulation program or model in the first stage, and determines the optimum problem solutions extracted from the results of the simulation model with the use of a suitable optimization algorithm in the second stage (Saadatpour and Afshar, 2007;Skardi et al., 2015;Mooselu et al., 2019a;Dehghanipour et al. 2020).
As stated in the previous sections, in this research the S-O approach was used to link the WLA model to a water quality simulation model. The simulator model used here was QUAL2K and the optimizer model was the ACO algorithm. The implementation stages of the S-O approach for solving the proposed model are presented in Fig. 3.
The following steps present a perfect description of this approach: Step 1. Before starting, the initial data extracted from Mostafavi (2010) have been entered into the simulator model. Initial data include the hydraulic characteristics of the river, water quality characteristics (such as the concentration of the untreated WQPs in each of the treatment units), and constant coefficients and parameters in the equations imposed on the water quality problems. Then, the optimizer program is implemented.
Step 2. At first, the number of program iterations and the number of required ants for running the program are specified. To do this, an initial pheromone trail is done and the amount equal to each one for all the paths is defined. There are different criteria for specifying the number of appropriate iterations for the optimizer algorithm. One of the criteria is to increase the number of iterations until an appropriate convergence is achieved in the path of all ants.   Step 3. One ant selects various scenarios for each of the treatment units based on a relatively random process. Now, the simulator model is used.
Step 4. Based on the selective scenarios, different levels of pollutants are removed by the units which discharge into the river.
Step 5. The simulator model runs based on the initial data (in Step 1) and the fractional removal levels of pollutants (in Step 4) and then returns the amount of the WQPs in different parts of the river as output data. Furthermore, the optimizer model is run.
Step 6. The index value of each of the WQPs is calculated for all river checkpoints. The relationship between the amount or concentration of WQPs and their respective index values are presented by the special curves (such as the one shown in Fig. 1). In order to use these figures in the optimizer model, after fitting more than 100 effective points of the curves, the equations of the index value based on concentration value for all WQPs with high precision are extracted and the optimizer model is given.
Step 7. The index value is calculated (based on Eq. 1) for each of the river checkpoints. In this step, if the considered ant is not the final ant, there is a return to Step 3; otherwise, it is continued from Step 8. Also, if Constraints 9 and 11 are not established, the next iteration should be followed and Step 3 should be started.
Step 8. At first, the average values of NSFWQI are calculated for all ants (according to Eq. 12). Then, these values are replaced in the fuzzy membership function and their maximum values are selected. If Constraints 8 and 10 are not established, the program should go to the next iteration.
Step 9. In this step, updating of the pheromone concentration takes place based on the maximum value of the objective function (according to Eqs. 4 and 5). All of the above steps will be repeated to the final iteration which was considered for the optimizer program.

Case study
To evaluate the performance of the proposed model, a hypothetical river system (Mostafavi, 2010) with a length of 30 km was used (Fig. 4).
As can be seen in Fig. 4, the river included 6 river reaches and 5 dischargers that discharge wastewater into the river. Here, there were 16 separate elements and, in the midst of these elements, there were the considered checkpoints. The inlet and outlet of the river were considered as the checkpoints; therefore, there were 17 checkpoints in the river system. The river has a wide rectangular cross-section and the average air temperature in the location of the river was 20°C (Mostafavi, 2010). https://doi.org/10.17159/wsa/2021.v47.i1.9444 To determine the re-aeration coefficient, the O'Conner and Dobbins equation was used (Sasikumar and Mujumdar, 1998). The hydraulic characteristics of the river are given in Table 1. The characteristics of the untreated wastewater for each of the treatment units as well as the river upstream are given in Table 2. The concentration of the total suspended solids (TSS) given in this table plays an important role in the calculations of the simulator model. However, it is ineffective in the determination of the WQI value. Fractional removal levels associated with different scenarios and wastewater treatment costs associated with the treatment units are given in Tables 3 and 4, respectively. The number of parameters intended for the proposed and the optimizer model can be seen in  Note et al., 1975;Sasikumar and Mujumdar, 1998;Eheart and Ng, 2004;Afshar et al., 2009;Mostafavi, 2010) and trial-and-error during the model execution.

RESULTS AND DISCUSSION
In this section, the results of the proposed model implementation are investigated and analysed. The model was executed on a hypothetical river (see Fig. 4) with the specifications and data given in Tables 1 to 5. This model was analysed by 30 successful and acceptable implementations of different values of q 0 and γ.
As mentioned before, q 0 is the regulatory parameter of the ACO algorithm and γ is the exponent of the membership function of average NSFWQI. At first, the obtained solution from the QUAL2K was examined with 100 and 150 iterations. As complete stability was not achieved in the solution (i.e., the ant path) in some cases, 200 iterations of parameters q 0 and γ were used. As a result, the QUAL2K model was underwent 200 iterations for each optimal solution. Also, 20 artificial ants were preferred for searching the random paths in the optimizer algorithm. The objective function values obtained with 30 implementations of the optimizer program for each of the various values of q 0 (0.6, 0.7, 0.8, and 0.9) and γ (0.5 and 2) are shown in Figs 5 and 6.
As can be seen, the objective function values was highest for q 0 = 0.6 and lowest for q 0 = 0.9. For q 0 = 0.6 and q 0 = 0.7 in different implementations of the optimizer program, these values were similar and, comparatively, identical. Based on these figures, reduction of q 0 from 0.9 to 0.6 led to a decrease in exploitation and increase in exploration in searching for the next path by ants. In other words, in this model, ants prefer to search more newly discovered paths and there is less likelihood of choosing the paths with a higher concentration of pheromone on which the best ants have passed. Based on these figures, it can be inferred that the objective function values obtained from the model equal to 0.6 and 0.7 were more appropriate, and the exploitation phenomenon was prominent in determining the optimum solutions. Statistical analysis of the data from Figs 5 and 6 and other outputs of the program is presented in Table 6.
In Columns 3 to 7, statistical elements of objective function values for 1 to 30 program implementation iterations are given. Columns 8 to 10 include the data associated with the maximum objective function values (based on data in Column 3). In Columns 11 to 15, the selective scenario number is given for the treatment units according to the maximum values of the objective function.    These data can help in understanding and analysing the model. A brief analysis of the data in the table is offered below: 1. The appropriate solutions in the proposed model can be the same as the maximum values of the objective function which are obtained for q 0 equal to 0.6 and 0.7. Also, in this case, and for each of the γ values, the maximum values of the corresponding objective function had the same values. This is why the data of Columns 9 to 15 demonstrated the same values for values of q 0 equal to 0.6 and 0.7. These q 0 included the same values of the average index, the same total cost, and also the same treatment scenario number.
2. For q 0 equal to 0.6 and γ values equal to 0.5 and 1, the solutions were obtained by 30 runs of the optimizer program. Therefore, it is concluded that obtaining appropriate solutions is possible by implementing a few runs of the optimization model, and this can increase the chances of finding the solutions. Besides, the standard deviation values for γ = 1 were less than that for γ = 0.5. So, the case of q 0 = 0.6 and γ = 1 can be considered as the best possible solution; therefore, the more appropriate objective function of the model was equal to 0.1565 and the total wastewater treatment costs for all the treatment units was 3.98 million USD/yr and the optimal scenarios related to the units of 1 to 5 were 10, 10, 10, 10 and 11, respectively.
3. For q 0 equal to 0.6 and values of γ equal to 0.5 and 1, the maximum value of the objective function was equal to the mode value, i.e., the maximum value of the objective function had the highest iteration obtained from running the optimizer program 30 times, while the maximum value of the objective function was repeated only once according to Column 8 in Table 6. On the other hand, for q 0 = 0.9 and γ = 2, the minimum value of the objective function was equal to the mode value; generally, it can be observed that by increasing the value of q 0 from 0.6 to 0.9 the program deviates from the optimal solutions and searches the nonoptimal solutions. It is recommended that, for high values of q 0 , the number of program iterations is extremely high, although the maximum obtained values might not be very suitable.
4. As is evident from Points 1 and 2, the set of obtained solutions associated with q 0 values equal to 0.6 and 0.7 can be considered as the optimal solutions of the model. Also, because of the negligible difference between the maximum and minimum values of the objective function for two of the q 0 , these show the lowest standard deviation among all q 0 values studied. The same results can be observed for q 0 = 0.9 and γ = 2. The standard deviation values are shown in Column 6 of Table 6. This means that, under these conditions, less time and a lower number of the program's executions (less than 30 times) can lead to near-optimal solutions (not necessarily optimal). 5. According to the contents of Point 3 and the data presented in Table 6, it can be stated that for the values of q 0 equal to 0.6 and 0.7, the standard deviation of the data was lowest and the difference between the values of the model and maximum values of the objective function was less related to other conditions of the problem. Therefore, it can require less time and fewer program iterations (less than 30) to obtain optimal solutions. Also, for q 0 values equal to 0.8 and 0.9, in which the standard deviation of the data was high and the difference between mode value and the maximum value of the objective function was higher than those for other q 0 values, the iteration number of the optimizer model for obtaining the optimal solutions becomes much higher. In this condition, the obtained maximum solutions were not optimal.
6. Figs 7 and 8 show the curves related to the NSFWQI values for 17 checkpoints of the river associated with the maximum values of the objective function for different values of q 0 and for values of γ equal to 0.5 and 2, respectively.
The values of the WQI index between Checkpoints 1 to 3 are the same in all cases because based on the hypothetical river structures (Fig. 4), the first wastewater treatment unit was located on the border of Checkpoints 1 to 3. So these checkpoints had the same qualitative behaviour. There are large jumps at Checkpoints 3, 6, 10, 12, and 14 which showed the establishment of 5 wastewater treatment units before the border checkpoints and drastically altered the behaviour of the river water quality parameters.

CONCLUSION
In this study, a WLA optimization model using the constraints of the wastewater treatment costs along with the objective function of the average NSFWQI and fuzzy approach was proposed. The use of the NSFWQI accelerates the understanding of the water quality in rivers and the fuzzy approach can decrease the uncertainties related to the model goals. Also, the S-O approach with the QUAL2K simulator model and the ACO algorithm were applied for solving the problem. The use of the S-O approach accelerates the program's execution. The model results showed that for the regulatory parameter of the ACO algorithm,