Waste load allocation under uncertainty using game theory approach and simulation-optimization process

This article aims to present a new methodology for waste load allocation (WLA) in a riverine system considering the uncertainty and achieve the lowest amount of inequity index, cost, and fuzzy risk of standard violation. To find a surface of undominated solutions, a new modified PAWN method, initially designed for sensitivity analysis, was developed and coupled with a simulation-optimization process using multi-objective particle swarm optimization (MOPSO) algorithm, to consider the uncertainty of all affecting variables and parameters by using their probability distribution. The proposed methodology applied to Sefidrood River in the northern part of Iran. Graph model for conflict resolution (GMCR) as a subset of game theory was implemented to attain a compromise on WLA among the stakeholders of a river system’s quality in Iran: Department of Environment, Municipal Waste Water, and Private Sector. Some undominated solutions were used in GMCR model and modeling the conflict among decision makers reveals that their preferences and the status quo do not lead to a solely stable equilibrium; thus the intervention of a ruler as arbitrator leads them to reach a compromise on a scenario that has a median FRVS and cost. Sensitivity analysis was done using the PAWN method to assess the sensitivity of three intended objectives to all variables and parameters.


GRAPHICAL ABSTRACT INTRODUCTION
Increasing human demand for water has made this vital element noteworthy, and its ownership rights to use have become more challenging. Using water bodies, especially rivers, for industrial purposes such as cooling the industries and assimilating the wastewater, is one of the consequences of the industrialization era. Dwindling the ecological capacity of the rivers to treat wastewater, population growthwhich leads to producing more wastewaterand industries bloom located along the rivers have reduced the quality of the river's water. This reduction is due to the drainage of the wastewater from the riverside dischargers into it. The requirements to maintain water quality in compliance with the approved standards and preserve its capability to receive the maximum flow of wastewater from industries has led to allocating the shares of wastewater discharge to the dischargers as a controversial and problematic issue. In other words, waste load allocation (WLA) in rivers refers to the determination of required pollutant treatment levels to a group of pollution dischargers to ensure that water quality standards are maintained along the river (Ghosh & Mujumdar ). Determining the shares for each of the dischargers, whether industries or municipal wastewater treatment facilities or non-point sources at the river border to discharge the wastewater, can be accomplished using var- Application of game theory in water quality management issues is preferred to conventional decision analysis methods because of two reasons. Firstly, unlike multi-criteria decision making (MCDM), this approach addresses decision-maker behavior and finds solutions that will not be refused by them. Finally, this method, unlike MCDM, is capable of expressing and reflecting the various characteristics of water resource problems, such as engineering, socio-economical, and political issues (Mirchi et al. ).
Optimal waste load allocation implies that the treatment levels not only maintain the water quality at standard level but also result in the best value for the objective function defined for the water quality management problem in the riverine system (Burn & Lence ). Chadderton et al. Burn & Yulianti () presented a methodology that contains three models, two of which are pairwise models (cost versus standard violation, the cost versus equity) to allocate waste in a riverine system. Feizi Ashtiani et al.
() studied an approach to meet environmental, economic, and equity objectives simultaneously. The waste loads were allocated according to the discharge permit market. They used two discrete graphs of cost versus violation and cost versus inequity to select four scenarios based on these two models' Pareto Fronts. Saberi & Niksokhan () proposed a methodology using the decision support system (DSS) for the graph model for conflict resolution II (GMCR II), multi-criteria decision making (MCDM) analysis and the multi-objective particle swarm optimization (MOPSO) algorithm while considering cost and standard violation to allocate waste load among dischargers.
The lack of considering the equity as a third objective in the optimization model leads to the setup of a three objectives simulation-optimization model. The difference between the current work and those mentioned above is that ongoing work leads to a surface of the undominated solution, where all points on that surface are optimum in terms of intended objectives. In contrast, the others are two-axis diagrams where the feasible space of the answers is a curve. This study, instead of using a deterministic model that they used, applied a stochastic model, which uses the distribution function of all variables and parameters that take part in the simulation-optimization process.
In this article, a novel methodology proposed to allocate waste load among dischargers located adjacent to a riverine system equitably and cost-effectively, as well as encompassing the lowest level of fuzzy risk of violating the standard (FRVS) for biochemical oxygen demand (BOD). Achieving a precise allocation for removal percent so that three minimizing objectives attained with a porous ratio, taking into account the uncertainty of variables and parameters affecting the simulation model, is required. Three different criteria (objectives) to allocate the waste load to three dischargers, which have controversial preferences, lead to the use of the game theory instead of conventional optimization approaches.
A stochastic simulation-optimization process was run applying MOPSO optimization algorithms, and quality simulation derived from fate and transport equations to find undominated solutions. In the last section, the GMCRþ model was utilized to find a stable equilibrium for the conflict about selecting a specific scenario by three decision makers: Department of Environment (DOE), Municipal Waste Water (MWW), and Private Sector (PS).
The impossibility to attain one specific equilibrium through conflict resolution caused the entrance of the judiciary; this intervention leads to achieving one particular equilibrium from six selected scenarios from the simulation-optimization process. A sensitivity analysis (SA) by the PAWN method (Pianosi & Wagener ) has been carried out to assess the outputs' sensitivity to the inputs parameters and variables uncertainty in the simulation process.
Using an integrated simulation-optimization which considers all goals simultaneously to find the Pareto Front was not done before since the previous studies usually use two separated simulation-optimization process considering two goals (Burn & Yulianti ; Feizi Ashtiani et al. ). In this article, an integrated model was set up that considers the uncertainty of all parameters.
How to sum up with the obtained output series from stochastic simulation-optimization reliably and flexibly, whether the median, average, or any other statistics should be used, was one of the difficulties which led authors to use a risk-based index: FRVS. The way sensitivity analysis is used for three intended functions output was another challenge, which forced us to run a separate simulation optimization process that uses standard violation function as a substitute function for FRVS to assess the output sensitivity to the inputs. According to existing decision makers and status, the conflict did not lead to a single stable equilibrium; thus, an arbitrator interferes with the conflict resolution process and is forced to omit some of the available scenarios to achieve a unique stable equilibrium.
The current article includes four main sections: introduction, methodology, results, and conclusion. The introduction section comprises a summary of works done in the field of WLA. The methodology section has eight subsections. The proposed model is the first subsection of this part; four different objective functions are presented and discussed. The optimization algorithm, simulation-optimization model, sensitivity analysis are the next subsections.
At the end of this section, the PAWN methodwhich is used for sensitivity analysisconflict resolution and the case study are described in detail. The third section, which is named results, covers the simulation-optimization, conflict resolution, and sensitivity analysis results. At the end of this article, a conclusion will be presented.

METHODOLOGY
In this section, four objective functions, optimization algorithm, uncertainty analysis, simulation-optimization, sensitivity analysis, PAWN method, and conflict resolution, are discussed, and finally, the case study used to apply the methodology on is described. As depicted in Figure 1, the proposed model to WLA is made up of three parts: simulation, optimization, and conflict resolution. In the first section, parameters and variables used to solve the problem are defined.
These are collected from official reports and previous studies.
It is vital to determine the acceptable range for data by using a literature review and field study for variables. Although considering the seasonal uncertainty and extreme events is noteworthy and helps to have a cost-effective allocation, the lack of precise data series leads to consider the worst situation that covers all plausible conditions. The worst situation for water quality issues is the lowest flow, defined by the 7Q10 method and the highest temperature.
The river water quality simulation will be run in the deterministic mode. To define sensitive parameters and variables, a sensitivity analysis of water quality simulation was carried out, and sensitive inputs reassessed to be at their correct, acceptable range. In the optimization section, first of all, a stochastic simulation-optimization will be run, and optimal removal percent defined. The undominated surface (Pareto Front) of solutions will be depicted based on the triple objective functions which try to be achieved. The generated Pareto Front presents the feasible space for the problem's answer, and each point on this surface that will be selected has equal worth in comparison to the other points on this surface.
In the next step of this section, six scenarios will be chosen from 100 scenarios, which generate Pareto Front, and each of these points contains removal percent for each discharger.
Finally, a stochastic simulation model will be run against 100 selected scenarios of removal percent to define the value for violating standards from chosen scenarios.
In the last section of the proposed methodology, a conflict resolution approach will be used to select the most acceptable scenario by stakeholders who are called players in this method. After defining players and their preferences based on the six selected scenarios, the impossible states and status quo will be set for the conflict resolution model. Stable equilibriums will be defined by analyzing the game.
If there was no stable equilibrium, the third-party intervenes in the conflict as an arbitrator and imposes its preferences on the game to lead to a sole stable equilibrium and put an end to the conflict.

Objective function
Three objectives have been considered to allocate waste load for this river in stochastic mode. In the deterministic mode, one of those objectives will be replaced by a fourth function based on the reasons that will be presented in the following parts: (1) minimizing amount of violation from the national standard of water bodies contamination concentration (standard violation) (Equation (1)), (2) minimizing total cost for wastewater treatment and violating standard penalties (Equations (10)-(13)), and (3) minimizing the inequity index (Equation (14)) are used in a deterministic way, and in addition to the first two, (4) minimizing the FRVS (Equation (15)

Standard violation
In this article, the objective of shielding water bodies' quality is to retain the BOD level under a predefined standardthe Iranian standard of discharge to surface water.
In order to estimate the allowable level of BOD for dischargers, it is essential to use fate and transport equations for BOD to observe predefined criteria under different scenarios. Equations had been applied to estimate the decomposition of BOD and change in the dissolved oxygen (DO) level in the river system. Based on these equations, the standard violation will be estimated by Equations (6) and (9) will be used to calculate DO at the intended location.
Minimize v j (1) Subjected to: In the abovementioned equations, v j is the standard violation at point j that is set at the end of the river path. It is a function of waste load removal percent (r) for each discharger, waste load weight (w), flow (Q), temperature (T ), reaeration and decomposition rates (k) as well as national standards of water quality (wq std ). Cost is a function of the fee paid by dischargers to treat their wastewater and a penalty relating to violating water quality standards.
L is the concentration of oxidizable matter, which here is BOD in mg/L (expressed in oxygen level), L 0 is the initial concentration of oxidizable matter in mg/L at the upstream, k r is the total removal rate in d À1 , which is made of decomposition rate (k d ) and settling rate (k s ). k a is reareation rate (d À1 ), o s is the concen- Saberi & Niksokhan ). They used the following equations to estimate the treatment cost of municipal and industrial wastewater treatment. It has to be considered that for each discharger (PS or MWW), treatment cost and penalty will be calculated and used in the conflict resolution phase. In contrast, the total PS and MWW's costs and penalties have been used in the optimization process.
In the abovementioned equations, c i is the cost pay for r i removal percent treatment at the ith discharger, p i is the penalty applied to this discharger for violating the standard at the checkpoint. The second equation (Equation (11)) is used to estimate the cost for municipal wastewater treatment and the third (Equation (12)) belongs to industries cost to treat their wastewater. In these two equations, c i is the total payment for the treatment of one cubic meter of wastewater in one year (y) in USD. (13)) based on the replacement cost method and assumed five different combinations of BOD level that at each of which, the BOD level at checkpoint exceeded the allowable amount in the regulation, and eventually estimated the penalty function by using regression. In the current paper, the penalty divides equally among all of the dischargers, similar to Saberi & Niksokhan (). In the following equation (Equation (13)) v (mg/L) is a representative for the amount of the violating BOD standard, which based on the Iran national standard equals 5 mg/L in water bodies, and the penalty cost is determined in USD per year.

Inequity index
Burn & Yulianti () used the inequity concept, which is a problematic issue in WLA and is the result of a desire for fairness in the distribution of the treatment level and the associated costs among dischargers. They used the inequity index (Equation (14)) to take fairness into account. A decrease in the amount of this index indicates the improvement of fairness among dischargers, and a rise shows the decline of fairness. In this relation, NS is the number of dischargers, r i and w i are removal percent and input waste load of the ith discharger, respectively. w i can be defined by the multiplying flows (q) of the ith discharger by its BOD concentration. r and w are the average removal level and the average input waste load for all of the NS dischargers.

Fuzzy risk of violating standard
Risk assessment is used extensively in water resources management fields such as water shortage ( (16)) indicates its equation. If the BOD concentration at checkpoint exceeds than 8 mg/L, the fuzzy membership number (μ w (c)) is one. If BOD is less than 5 mg/L, then μ w (c) is zero.
In the case where the BOD concentration at the checkpoint is in the range of 5-8, the μ w (c) will be between zero and one.
The FRVS value declares the risk that BOD or any contamination level proceeds at the predefined level at the checkpoint due to the waste load discharged to the system. FRVS with a low value presents a lower risk for water quality. While zero FRVS is representative of a situation in which the river will not experience any violating from the predefined standard, FRVS equals 1 indicates that the situation at which the BOD level exceeds the allowable level, is inevitable. Therefore, to decrease this risk, either the probability of violating the standard or the level of this violation should be reduced.

Optimization algorithm
The use of evolutionary algorithms for multi-objective optimization problems has been increased and has led to Global, quantitative, model-independent, unconditional on any assumed input value, easy to interpret, easy to compute, stable, and moment independent are those features which a sound sensitivity analysis index should satisfy (Liu & Homma ). While variance-based sensitivity indices satisfy all the above features, but the last one, the PAWN method, meets all features (Pianosi & Wagener ).

PAWN method
In this article, the authors use the PAWN method ( PAWN measures the distances between conditional and unconditional CDF of outputs for a specific parameter in which this method is trying to estimate its sensitivity index. In the first phase of estimating this distance, selecting all parameter and variables to calculate CDF of the output y (F y (y)) will be done in a stochastic manner, and every value in the region of the allowable space of inputs has the chance to be selected and used; this kind of CDF is called unconditional CDF. To create a conditional CDF, parameter x, the parameter which the model tries to measure output sensitivity relative to it must be set to a specific value and selection of the other parameters and variables will be stochastic (F yjx i (y)). As one of the fundamental properties that a 'good' global sensitivity index should satisfy is to be unconditional about any assumed input value (Liu & Homma ), PAWN uses n different values for x to estimate n conditional CDFs and address this issue (Pianosi & Wagener ).
As is evident in the third box of Figure 3, there is an infinite amount of distance between an unconditional and a conditional CDF; therefore, the PAWN uses the Kolmogorov-Smirnov (Kolmogorov ; Smirnov ) statistic (KS) to calculate this length (Equation (17)). Since there are n different conditional CDFs and each one has a specific distance to the conditional CDF, this method uses the maximum of the KS statistic (T i ) set to declare the final result, which is the x sensitivity index (Equation (18)). Steps in the numerical implementation of the PAWN sensitivity index can be found in Figure 3:

Conflict resolution
The existence   Iran could lead to a stable equilibrium.

Case study
The proposed methodology was applied to the Sefidrood river in Gilan province, in northern Iran, to assess its applicability and efficiency. This river is 114.95 km length from the Sefidrood Dam outlet to the river meadows at the Caspian Sea. The river is divided into nine reaches according to its hydraulic characteristics. Twelve industries and municipal wastewater are discharging to the river. In this study it is supposed that each discharger uses its facilities to treat wastewater. Figure 5 depicts the location of Sefidrood and three tributaries that join the mainstream of this river: Tootakbon, Tarikrood, and Zilakirood.

RESULTS
This section provides the results of different processes and models that are used in the proposed methodology in three subsections. In the first subsection, simulation-optimization results are presented, then the results of sensitivity analysis are described, and ultimately, the final solution for the WLA is presented through a conflict resolution subsection.

Simulation-optimization
One hundred undominated solutions, demonstrated in simulation-optimization model, which used MOPSO as an evolutionary optimization algorithm. Figure 6(b)-6(d) display the undominated solution in a 2D view from different perspectives. Six scenarios (Table 1)  scenario has the lowest value in its acceptable range for FRVS and a low value for the inequity index. The S4 scenario, which is representative of a highly risky situation, has a low inequity index and a median cost. Nevertheless, the decline of FRVS will lead to an increase in the system's cost.

Sensitivity analysis
In the current study, all 76 parameters and variables used in the stochastic simulation optimization process were used to be assessed in the sensitivity analysis. A critical point that should be considered is that instead of analyzing the FRVS, authors used standard violation as an objective function to assess its sensitivity; the reason is related to the nature of these two different functions, despite their similarity. The second phase of sensitivity analysis leads to the generation of a series of N values that could be used to make a CDF or PDF, hence comparing the different CDFs  of numbers. One has to use the model output's PDF to calculate the fuzzy risk. In contrast, in the second phase of sensitivity analysis one has to reach a PDF or CDF, and this is not possible, hence instead of using FRVS, which is the final outcome of uncertainty analysis and is a scalar value, one has to use a function whose outcome of uncertainty analysis will be a series of scalar values instead of a single scalar value. Due to this point, a standard violation function, which has a type the same as p(c) in the FRVS (Equation (7)    Based on the sensitivity index driven from the PAWN method, which is depicted in Figure 10, the inequity index function has the highest sensitivity to the BOD concentration and flow of the ninth discharger's (BOD9 and Q9), the seventh discharger's flow (Q7), and the ninth discharger's flow (Q9), respectively. The order of the parameters sensitivity index indicates that after the previous parameters, the function has the most sensitivity to the different dischargers' Q and BOD concentrations in the next place. It should be noticed that the upstream discharge (Q Upstream ) has the least sensitivity index, among other parameters and variables. Since the definition of inequity (Equation (14)) is just dependent on the amount of waste load (w), and removal percent (r), such ordering of the sensitivity index can be justified.
According to the fact that the sensitivity index (the average of Ti index) is a dimensionless value, if the average of this index is calculated, it will reveal that the proposed methodology has the most sensitivity to upstream flow (Q Upstream ), the ninth discharger's BOD concentration and flow (BOD9 and Q9), and the seventh discharger's flow (Q7). It has also revealed that the methodology has less sensitivity to the reareation rate at the second reach of the river.
An overall view of the sensitivity analysis for all 76 inputs is displayed in Figure 11.
To standard; the CDF for this conditional output is drawn up in black color at the end of the graph.

Conflict resolution
Two of the intended objective functions, cost, and FRVS, are in contrast to each other, so that achieving better performance for one of them requires getting away from the best performance for another one. Reducing the FRVS needs higher treatment percentage, and hence increases dischargers' costs, therefore reducing the system's total cost will lead to a higher FRVS amount.  On the other hand, DOE only considers declining waste load discharged from municipal and private wastewater treatment facilities to the water bodies to increase water quality, specifically in checkpoints.
The third player in this conflict is the private sector (PS), which is representative of a coalition of industries surrounding the river and owned by the non-governmental sector.
This player or DM aims to maximize its profit or minimize its cost; thus, it has specific preferences, which would be in contrast to other players' preferences. also try to reduce their costs; thus, they will sort their preferences based on the cost that they must pay for scenarios.   thus, the current status quo alongside the players' preferences will not lead to a common option selection by two other players, and a stable equilibrium will not be achieved.
For instance, PS comprehends that selecting state number 58 is in its favor; thus, select this state, but two other players try to choose other states, which has more payoff for them in comparison to status number 58.
The intervention of a third party was assessed to comprehend whether it is useful to resolve the conflict.
Hence the judiciary entered the conflict as an arbitrator who can force players not to select specific options, including the options with the highest FRVS and the most significant inequity index. The judiciary will ban players from selecting S2 and S4 as the options, with the highest amount of FRVS and S1 as the option with the highest inequity index.
Removing the afore-mentioned triple options causes a change in the status quo, and selected options will move to the nearest ones. In this state, which is representative of the worst situation and numbered 23 by GMCR þ , MWT still select S5, DOE selected option will move from S4 to S5, and finally, PS will select S6.
Analyzing the new conflict results from the entrance of arbitrator to previous conflicts shows that all possible stable equilibriums will decline to nine items. By considering the new status quo, GMCRþ outputs prove PS will move from S6 to S5, and consequently, two other players comprehending moving from S5 to another option will not increase their payoffs, and thus they do not change their selection. This state is named number 14 by the model. Therefore, one movement from status number 23 to status number 14 by PS will lead to a stable equilibrium, which is the final resolution ( Figure 15). Table 4 shows the twelve dischargers' removal percent for stable equilibrium, which is the final WLA based on the predefined objective functions and triple decision makers' preferences.
S5 will be the option in which the entrance of an arbitrator causes it to be accepted by all parties of the conflict. This option leads to a total cost of 1,298,112 USD per year and MWW should pay 1,137,180 USD, and PS should pay 160,932 USD. Selecting this option also leads to an inequity value of 14.65 and FRVS of 11.29.

CONCLUSION
Three objectives were considered in this article simultaneously to allocate waste load quotes among dischargers located adjacent to a river. Minimizing dischargers' cost for their wastewater treatment, as well as the penalty for violating the standard, was one of those objectives. Fuzzy risk of violating the standard (FRVS) was the second function, which represents the amount of the risk that may exist in the riverine system to infract the predefined BOD level at the checkpoint due to the wastewater discharged. The