A Multi-Objective Chance-Constrained Programming Approach for Uncertainty-Based Optimal Nutrients Load Reduction at the Watershed Scale

A multi-objective chance-constrained programming integrated with Genetic Algorithm and robustness evaluation methods was proposed to weigh the conflict between system investment against risk for watershed load reduction, which was firstly applied to nutrient load reduction in the Lake Qilu watershed of the Yunnan Plateau, China. Eight sets of Pareto solutions were acceptable for both system investment and probability of constraint satisfaction, which were selected from 23 sets of Pareto solutions out of 120 solution sets. Decision-makers can select optimal decisions from the solutions above in accordance with the actual conditions of different sub-watersheds under various engineering measures. The relationship between system investment and risk demonstrated that system investment increased rapidly when the probability level of constraint satisfaction was higher than 0.9, but it reduced significantly if appropriate risk was permitted. Evaluation of robustness of the optimal scheme indicated that the Pareto solution obtained from the model provided the ideal option, since the solutions were always on the Pareto frontier under various distributions and mean values of the random parameters. The application of the multi-objective chance-constrained programming to optimize the reduction of watershed nutrient loads in Lake Qilu indicated that it is also applicable to other environmental problems or study areas that contain uncertainties.


Introduction
Lake eutrophication is one of the greatest environmental challenges, which hampers and destroys ecological functions of water bodies, and leads to deterioration of water quality [1,2].Many irreversible effects on aquatic ecosystems and drinking water supplies have resulted from degeneration of water quality, and this issue has become the principal bottleneck to the sustainable development of watersheds and the protection of human health [3][4][5].Watershed planning and management of water quality is imperative.Load reduction of watershed nutrients (mainly Nitrogen and Phosphorus) is an effective approach to water quality improvement.Load reduction is also an important approach to balancing watershed development and protection of water quality, because nutrient loads to water bodies are mainly due to human production and activities, and they are a primary cause of water quality impairment [6,7].Optimization models are applied widely to provide effective load reduction schemes as a technical basis for watershed management [8][9][10].Optimization models for watershed load reduction provide pollution control schemes that can produce efficient management of watershed load reductions under the constraint of watershed development and human demands.
A multitude of programming methods have been proposed to deal with load reduction for protection of water quality [11][12][13].Among these, much research has focused on representing and dealing with uncertainties, because uncertainties in the components and related impact factors within the system cause optimization problems to be more complicated [14,15].On the other hand, studies on multi-objective models continuously spring up, because taking the diverse demands of social economic progress and environmental protection into account it is necessary, and using multi-objective models can balance the conflict between different optimal objectives, such as economic development and restoration of water quality [16][17][18][19].Uncertainty in the system also resulted from some influence factors having random characteristics, which originates from natural processes in the hydrological cycle, and variation in water pollution over time and space [20][21][22].That is to say, uncertainty and multi-objective problems should be considered simultaneously, however, few studies have been conducted on optimization of load reduction in watersheds while considering both uncertainty and the weighting of multiple objectives.Not considering the balancing of multiple objectives or the handling of uncertainty results in a shortage of information in the optimization process, and schemes that are based on inapposite fundamental assumptions will lack the capacity both to reflect on practical situations and to provide a planning scheme for decision-makers.Thus, a multi-objective programming method that accounts for uncertainty in reducing pollution loads to improve water quality is imperative and practical.
In this study, a multi-objective chance-constrained programming integrated with Genetic Algorithm (GA) and robustness evaluation methods was proposed to optimize load reduction of nutrients in the Lake Qilu watershed.The objective function was composed of minimal system costs and maximization of water qualification rate.Typical impact factors (Total Nitrogen and Total Phosphorus, TN and TP) were incorporated, and the primary pollution sources, such as agricultural non-point pollution, urban non-point pollution, and livestock pollution, were included.From this, constraints in the optimal model were selected, such as reduction of TN and TP from agricultural non-point pollution, domestic sewage treatment level, municipal solid waste disposal, and ecological restoration.Optimization schemes were obtained to find the tradeoff between system investment and risk, and estimation of the model's robustness was conducted.

Study Area
Lake Qilu is one of the nine plateau lakes in Yunnan Province, China, and it is also the natural foundation for development of Tonghai Country, Yuxi City, which has been seriously contaminated through natural processes and human activity in recent years.TN and TP are important indicators of eutrophication for Lake Qilu, and nonpoint-source pollution from agricultural systems and livestock breeding are the primary sources of pollution.TN from agricultural non-point source pollution (NPS) and livestock breeding are 33.04% and 37.56%, respectively, and TP from agricultural, nonpoint-source pollution and livestock breeding are 19.80% and 61.68%, respectively.In addition, household waste is another major component of non-point source pollution, which contributes 23.17% and 14.70% to TN and TP, respectively.Thus, measures to reduce nutrient loads are necessary to control watershed pollution in the Lake Qilu watershed [23].
The Lake Qilu watershed contains five sub-watersheds, the Hongqi River sub-watershed, the Zhewan River sub-watershed, the Zhong River sub-watershed, the Daxin River sub-watershed, and the Northest sub-watershed, which are selected as basic units for optimizing load reductions (Figure 1).According to watershed management requirements, "The 12th Five-Year Plan", and pollution characteristics of Lake Qilu, reduction in TN and TP from agricultural, non-point pollution, Water 2017, 9, 322 3 of 18 sewage treatment of urban areas and villages, garbage disposal, and ecological restoration of riparian buffers were selected as constraint conditions.Correspondingly, municipal sewage treatment, domestic rubbish disposal, agricultural sewage collection and treatment, soil testing and formulated fertilization, constructed wetlands, and ecological restoration were chosen as engineering measures.Among these measures, sewage treatment plants and constructed wetlands can reduce the sewage containment.Agricultural sewage collection and treatment, soil testing and formulated fertilization, and constructed wetlands contribute to the reduction of TN and TP.Domestic rubbish disposal can deal with municipal solid waste, which increases with the growth of the urban population.Ecological restoration is beneficial for restoring the ecosystems of lakes and rivers generally.
Water 2017, 9, 322 3 of 17 treatment of urban areas and villages, garbage disposal, and ecological restoration of riparian buffers were selected as constraint conditions.Correspondingly, municipal sewage treatment, domestic rubbish disposal, agricultural sewage collection and treatment, soil testing and formulated fertilization, constructed wetlands, and ecological restoration were chosen as engineering measures.Among these measures, sewage treatment plants and constructed wetlands can reduce the sewage containment.Agricultural sewage collection and treatment, soil testing and formulated fertilization, and constructed wetlands contribute to the reduction of TN and TP.Domestic rubbish disposal can deal with municipal solid waste, which increases with the growth of the urban population.Ecological restoration is beneficial for restoring the ecosystems of lakes and rivers generally.

Chance-Constrained Programming
Stochastic programming (SP) takes randomness of parameters into consideration.Chance-constrained programming (CCP) is one of the typical methods of SP, which requires that constraints be satisfied at a certain confidence level (αi).The general form of CCP (max) can be expressed as follows [24]:
where x is a decision vector, ξ is a stochastic vector, f (x, ξ) is an objective function, gi (x, ξ) is a set of constraints, and αi and β are the confidence levels.
Generally, there are two methods to solve CCP: (a) Transform CCP into a deterministic model.
The conventional approach for transforming CCP is to translate the chance constraints into the approximated linear forms, and then the linear forms can be solved by traditional mathematical

Chance-Constrained Programming
Stochastic programming (SP) takes randomness of parameters into consideration.Chance-constrained programming (CCP) is one of the typical methods of SP, which requires that constraints be satisfied at a certain confidence level (α i ).The general form of CCP (max) can be expressed as follows [24]: where x is a decision vector, ξ is a stochastic vector, f (x, ξ) is an objective function, g i (x, ξ) is a set of constraints, and α i and β are the confidence levels.
Generally, there are two methods to solve CCP: (a) Transform CCP into a deterministic model.
The conventional approach for transforming CCP is to translate the chance constraints into the approximated linear forms, and then the linear forms can be solved by traditional mathematical models [25][26][27].For objective functions in linear form and some specific nonlinear functions in which parameters follow a normal distribution or an exponential distribution, approximated linear forms can be obtained, and CCP can be transformed into its relevant deterministic programming [28,29].
Stochastic simulation (i.e., Monte Carlo Simulation) is used widely in handling stochastic programming [30].According to the fundamentals of Monte Carlo Simulation, the probability of an event can be estimated by frequency that is gained from numerous tests, and the frequency is probability when the sample size is large enough.For CCP: where ξ is a stochastic vector, and the cumulative probability distribution is Φ (ξ).If we set independent random variables (ξ 1 , ξ 2 , . . ., ξ N ) that are acquired from Φ (ξ) in N experiments, then N is the number of times that the following formula holds: N is the number of random variables that satisfy constraints, and probability can be estimated by N /N.Thus, CCP holds when N /N ≥ α.

Multi-Objective Programming
The expression of multi-objective programming (MOP) is as follows [31]: Generally, a tradeoff among sub-objects in MOP is necessary to achieve an overall optimization under the conflict of the sub-objectives.The essential difference between MOP and single object optimization is that the optimization solution is not exclusive.The solution set of MOP is composed of a series of Pareto optimal solutions.
The conventional approach for solving MOP is to convert the multi-objective problem to a single-objective problem.Typically, this includes the constraint method, the weighting method, and the distance function method [32].Constraint method converts MOP to a SOP through transferring secondary objectives to constraints that obey a certain requirement.Weighting method establishes a new objective function with original sub-objectives and the weight corresponding with them.Distance function method can search optimal solution that have the shortest distance for each sub-objective in solutions space.Traditional methods inherit the mechanism of single objective optimization, but limitations are obvious at the same time.First, the conventional approach, such as the weighting method, is too sensitive to the optimal frontier of a Pareto solution, and it is strongly subjective with respect to weighted value distributions.Second, traditional methods obtain one solution at a time, thus, they require an enormous amount of time, and it is difficult to obtain an effective decision [33].

Genetic Algorithm
Genetic Algorithm (GA) is a global searching and optimization method developed from simulations of the natural evolutionary process that is based on Darwin's theory of evolution and Mendel's theory of heredity.GA has been used widely in dealing with MOP due to its advantage of being able to handle a large search space and for obtaining a series of Pareto optimal solutions at one time.The basic process of GA is as follows [34]: GA is an efficient algorithm for solving the MOP problem, because it is possible to get multiple, optimal solutions and to acquire an optimal solution set at the same time.A series of new algorithms has been developed from GA, such as Vector Evaluated Genetic Algorithm (VEGA), Multi-Objective Genetic Algorithm (MOGA), Non-dominated Sorting Genetic Algorithm (NSGA), and so on [35,36].VEGA divides the population into M subpopulations and conduct SOP for each subpopulation, and then repeat the process "population partition-optimization searching for each subpopulation-exchange of information between generations" until the optimal solutions are obtained.Although VEGA is convenient in operation, it has the apparent defect that the optimal solutions are generally the extreme points on the Pareto frontier and the Pareto solutions between the extreme points are unable to get.MOGA seeks optimal solution by calculating fitness for each individual, which introduces a fitness-sharing model to apply different types of optimization problems.The deficiency of MOGA is that, if the selection in the Pareto solution is biased resulted from the fitness value, the shape of Pareto frontier will be sensitive to spatial distribution of the solutions.NSGA proposes a hierarchical algorithm for Pareto solution searching, which guarantees the diversity of solution space and the homogeneity of Pareto frontier.
In this study, we selected CE-NSGA-II and applied it to the optimization problem to reduce nutrient loads in a watershed.This algorithm has the advantage on large coverage of the solution set and a fast convergence rate, and has the ability to increase the population diversity under lower fitness value [37].

Robustness of Decision Scheme
Decades of practical experience indicate that not all the information can be described and modeled well by the tools of probability and statistics, and the tools of statistics modeling and analysis may fail to express and handle precisely and properly the uncertainties that are posed by a complex system.The phenomenon above is called "deep uncertainty" [38][39][40].Under deep uncertainty, robust decision-making was proposed to provide policies that are robust across multiple scenarios, or alternative models at low cost [41].Robust decision-making gives priority to robustness of the decision scheme rather than priority to optimization; that is to say, the decision scheme, which is insensitive to uncertainty parameters, will still exhibit good performance under a changing scenario.
Robustness can be expressed as system adaptability: where Λ (d, X) is adaptation function, d is decision scheme, Y T is threshold value, and X is deep uncertainty.In this study, estimating parameter RI is selected to assess robustness of the decision scheme, which is based on the hypothesis that every scenario happens with the same probability.Parameter settings of TN reduction of measure j (ANR j ) and TP reduction of measure j (APR j ) based on the above theory can be expressed as follows: Water 2017, 9, 322 6 of 18 2.3.Multi-Objective Stochastic Programming Model for Load Reduction in the Lake Qilu Watershed

Modeling Optimal Load Reduction
We developed an optimal model that used the Lake Qilu watershed as the management objective.We selected reduction of TN, TP, domestic sewage treatment and municipal solid waste; and ecological restoration and the probability of constraint satisfaction as constraint conditions.Among them, reduction of TN and TP were selected as uncertain factors in terms of random parameters, since the pollution load of TN and TP show an obvious random feature depending on season [42,43].We also set the maximization of constraint satisfaction probability and the minimization of system cost on load reduction as optimization objectives.Therefore, measures to reduce the costs of load reduction included the cost of reducing TN and TP from agricultural non-point pollution, the cost of sewage treatment for urban areas and villages, the cost of garbage disposal, and the cost of ecological restoration of riparian buffers.The expression of the multi-objective, stochastic programming model is as follows.
(a) Minimize system cost: (b) Maximize water qualification rate: In objective (b), UC m refers to the constraints containing randomness parameters (ANR and APR).
(a) TN from non-point source pollution in agricultural system: Reduction of TN from agricultural production (j = 3), rural domestic pollution (j = 4) and constructed wetland (j = 5) should meet the requirements of non-point, source pollution (NPS) control.
(b) TP from non-point source pollution in agricultural system: Reduction of TP from agricultural production (j = 3), rural domestic pollution (j = 4) and constructed wetland (j = 5) should meet the requirements of non-point source pollution control.
Water 2017, 9, 322 7 of 18 (c) Domestic sewage: Domestic sewage treatment level should meet sewage treatment requirements, and engineering measures related to domestic sewage is municipal wastewater treatment (j = 1), agricultural sewage collection and treatment (j = 3), and constructed wetland (j = 5).
(d) Municipal solid waste: The level of municipal solid waste disposal should meet the relevant planning requirements.
(e) Ecological restoration: The ecological restoration constraint refers to the percentage of ecological restoration that satisfies the relevant planning requirements.
(f) Technical constraint: Reduction of TN and TP may include random effects that are affected by many uncertainty elements, thus, it is necessary to quantify the random information effectively to obtain precise and proper schemes for decision-making.Assume that reduction of TN and TP obeys a normal distribution, where the mean value is µ and standard variation is σ.µ is obtained from the Manual of Engineering Measure Design, and σ is set as σ = µ × 10% (shown in Table 1).The constraint for CCP can be expressed as follows: In traditional CCP, the confidence level (α i ) is set by decision-makers, nevertheless, setting α i too high will lead to high system cost, because decision-makers cannot predict precisely the relationship between α i and system cost.In this study, the probability of constraint satisfaction (F P ) was set as one of the objective functions, and the single objective model is transformed to multi-objective, chance-constrained programming.Thus, decision-makers can consider system cost and risk comprehensively without the limit of a pre-set confidence α i .
Table 1.Definitions of the parameters used in the model to optimize nutrient loads in the Lake Qilu watershed.

Determination of Parameters
Definitions and values of parameters were obtained from the relevant literature and documents, such as statistical yearbooks, The 11th Five-Year Plan of Lake Qilu, the 12th Five-Year Plan of Lake Qilu, The First National Pollution Census, Manual of Pollution Discharging and Production Coefficient from Domestic Waste in The First National Pollution Census, Design Manual of Constructed Wetland, and so on.Essential data, including water quality target, nutrients load, and engineering measures were acquired from the above materials, which provided a database for determining values for decision variables and parameters in the optimization model.Nomenclature and values for parameters in the model are provided in Tables 1-4.(a) We transformed ANR j and APR j into deterministic parameters and ran the model under the laxest constraints, where confidence level was 99.9%, and we obtained the optimal solution (X min ).
Objective function: Constraints: X i,9 • a i,9 ≥ RFL i (25) (b) We transformed ANR j and APR j into deterministic parameters and ran the model under the strictest constraints, where confidence level was 99.9%, and we obtained the optimal solution (X max ).Inequality ( 20) and ( 22) can transfer to Equations ( 27) and ( 28) respectively: (2) Step II: Stochastic simulation.
Water 2017, 9, 322 10 of 18 We handled uncertainty constraints using stochastic simulation, and we conducted sampling of ANR j and APR j based on the Monte Carlo method, setting test number N (N = 1000).
where t was the number of tests that met the requirement of Formula ( 29), and the objective Function (10) was transformed to Formula (30), and (30) was equivalent to Equation (31). (3) Step III: Solve multi-objective, chance-constrained, programming model.
For solving the high-dimensional linear optimization problem above, the multiple population genetic algorithm is capable to maintain diversity of the solutions and avoid getting local minimum value prematurely.In order to keep diversity of the population, we set three subpopulations to maintain diversity of the population, to ensure the rate of convergence, and to avoid local optimal solutions.Individuals in each subpopulation were set to 40, and genetic algebra was 200.We then added the extreme solutions [X min , X max ] into the population initialization process, and other populations were produced under stochastic constraints.

Overview of the Optimization Schemes
An optimization scheme from the multi-objective, chance-constrained, programming model for load reduction in Lake Qilu was obtained using CE-NSGA-II.From 120 sets of solutions, 23 of them were Pareto optimal solutions.Constraint satisfaction probability and system investment constituted the Pareto frontier, which demonstrated the tradeoff between constraint-violation risk and system cost.High system investment led to a high satisfaction level of the constraints, that is, constraint-violation risk of the system was relatively low (Figure 2).

Overview of the Optimization Schemes
An optimization scheme from the multi-objective, chance-constrained, programming model for load reduction in Lake Qilu was obtained using CE-NSGA-II.From 120 sets of solutions, 23 of them were Pareto optimal solutions.Constraint satisfaction probability and system investment constituted the Pareto frontier, which demonstrated the tradeoff between constraint-violation risk and system cost.High system investment led to a high satisfaction level of the constraints, that is, constraint-violation risk of the system was relatively low (Figure 2).To represent and quantify the relationship between system investment and constraint satisfaction probability, the change rate index (CR) was introduced: , , /10 ( ) /10 CR refers to the increment of system investment that corresponded to the constraint satisfaction To represent and quantify the relationship between system investment and constraint satisfaction probability, the change rate index (CR) was introduced: CR refers to the increment of system investment that corresponded to the constraint satisfaction probability from solution m to m + n, after sorting the system investment from small to large.The relationship between CR and F P indicated that there was a dramatic increment on system investment when constraint satisfaction probability was higher than 0.9 (Figure 3).Marginal effectiveness was obvious in schemes 3-7, in which a small investment reduced model risk significantly.Schemes 12-19 were ideal schemes, and decision-makers can consider flexible options regarding their preference for budget constraints or system stability.To represent and quantify the relationship between system investment and constraint satisfaction probability, the change rate index (CR) was introduced: CR refers to the increment of system investment that corresponded to the constraint satisfaction probability from solution m to m + n, after sorting the system investment from small to large.The relationship between CR and FP indicated that there was a dramatic increment on system investment when constraint satisfaction probability was higher than 0.9 (Figure 3).Marginal effectiveness was obvious in schemes 3-7, in which a small investment reduced model risk significantly.Schemes 12-19 were ideal schemes, and decision-makers can consider flexible options regarding their preference for budget constraints or system stability.

Analysis on Optimal Solutions
Schemes 12-19 were optimal solutions for both system investment and constraint satisfaction (Figure 3).The solutions under six engineering measures are shown as Table 5, taking Hongqi sub-watersheds, for example.

Sewage Treatment
Among the five sub-watersheds, only the Hongqi sub-watershed was suitable for establishing a municipal sewage treatment plant, because the other four sub-watersheds were in rural regions, where sewage disposal and constructed wetlands were used to treat sewage.
For the Hongqi sub-watershed (W 1 ), there were no significant differences between the capacities of municipal wastewater treatment plants between schemes 12-17, although the project scale for sewage collection and disposal increased and constructed wetland areas decreased from schemes 12-17.Schemes 18 and 19 exhibited differences from the others in size decline of sewage treatment plants.In particular, the size of sewage treatment plants and sewage disposal projects in scheme 19 were higher than for other schemes to ensure lower system risk (Figure 4).Among the five sub-watersheds, only the Hongqi sub-watershed was suitable for establishing a municipal sewage treatment plant, because the other four sub-watersheds were in rural regions, where sewage disposal and constructed wetlands were used to treat sewage.
For the Hongqi sub-watershed (W1), there were no significant differences between the capacities of municipal wastewater treatment plants between schemes 12-17, although the project scale for sewage collection and disposal increased and constructed wetland areas decreased from schemes 12-17.Schemes 18 and 19 exhibited differences from the others in size decline of sewage treatment plants.In particular, the size of sewage treatment plants and sewage disposal projects in scheme 19 were higher than for other schemes to ensure lower system risk (Figure 4).In W2 to W5 (Zhewan River sub-watershed, Zhong River sub-watershed, Daxin River sub-watershed, and Northest sub-watershed), the scale for sewage collection and disposal projects showed a small increase, and constructed wetlands increased significantly, because no municipal wastewater treatment plant has been built in these regions.The scale of constructed wetlands varied among different sub-watersheds, and decision-makers can choose alternatives by taking into account any practical issues they might encounter within a particular sub-watershed (sub-watershed 2, for example, shown in Figure 5).In W 2 to W 5 (Zhewan River sub-watershed, Zhong River sub-watershed, Daxin River sub-watershed, and Northest sub-watershed), the scale for sewage collection and disposal projects showed a small increase, and constructed wetlands increased significantly, because no municipal wastewater treatment plant has been built in these regions.The scale of constructed wetlands varied among different sub-watersheds, and decision-makers can choose alternatives by taking into account any practical issues they might encounter within a particular sub-watershed (sub-watershed 2, for example, shown in Figure 5).

Load Reduction of TN and TP
The scope of measures for TN and TP reduction was also provided in the optimal scheme.In W1, constructed wetlands were small-scale in each scheme, but the scale of sewage collection, soil testing, and formulated fertilization increased sharply.The scale of soil testing and formulated fertilization were relatively small in schemes 12-14 and 18 in W1, although in schemes 15-17 and 19, it was much larger.In W5, the project size of rural sewage proposal and soil testing and formulated fertilization showed a linear growth.W1 and W5 were selected as examples (Figure 6).

Load Reduction of TN and TP
The scope of measures for TN and TP reduction was also provided in the optimal scheme.In W 1 , constructed wetlands were small-scale in each scheme, but the scale of sewage collection, soil testing, and formulated fertilization increased sharply.The scale of soil testing and formulated fertilization were relatively small in schemes 12-14 and 18 in W 1 , although in schemes 15-17 and 19, it was much larger.In W 5 , the project size of rural sewage proposal and soil testing and formulated fertilization showed a linear growth.W 1 and W 5 were selected as examples (Figure 6).

Load Reduction of TN and TP
The scope of measures for TN and TP reduction was also provided in the optimal scheme.In W1, constructed wetlands were small-scale in each scheme, but the scale of sewage collection, soil testing, and formulated fertilization increased sharply.The scale of soil testing and formulated fertilization were relatively small in schemes 12-14 and 18 in W1, although in schemes 15-17 and 19, it was much larger.In W5, the project size of rural sewage proposal and soil testing and formulated fertilization showed a linear growth.W1 and W5 were selected as examples (Figure 6).

Estimation of Robustness of the Optimization Schemes
In this study, random parameters were assumed to follow a normal distribution.However, the actual distribution was still unknown; that is, deep uncertainty existed in the model.Thus, estimation of robustness was imperative to ensure that optimal schemes were adapted to various scenarios.

Variation on Distribution
Robustness evaluation was conducted under uniform sampling based on RI (as shown in Formula ( 6)-( 8)). Figure 7 shows that robustness of the optimal scheme was acceptable as before, because the Pareto solution was still on the frontier under parameter variation.It would not be wise to be overly optimistic about system risk, while providing a decision scheme that had a lower theoretical probability of constraint satisfaction is appropriate.Besides, the constraint satisfaction probabilities were lower than 0.7 in schemes 12-18, when the parameters did not follow a normal distribution.
because the Pareto solution was still on the frontier under parameter variation.It would not be wise to be overly optimistic about system risk, while providing a decision scheme that had a lower theoretical probability of constraint satisfaction is appropriate.Besides, the constraint satisfaction probabilities were lower than 0.7 in schemes 12-18, when the parameters did not follow a normal distribution.

Variation in Mean Values of the Parameters
Mean values of parameters were obtained from actual measurements or estimations, which also brought uncertainties.Therefore, robustness evaluations on the optimal schemes based on different mean values evaluation was carried out.New mean values of ANRj and APRj, which decreased by 10%, are shown in Formula ( 33) and ( 34): (1 10%) According to the new mean values, the relationship between system investment and constraint satisfaction probability under the new distribution was obtained (Figure 8).The result indicated that robustness of the optimal scheme was very good.Not only were the solutions still on the Pareto frontier, but also the marginal effect of system investment on the constraint satisfaction probability increased rapidly.

Variation in Mean Values of the Parameters
Mean values of parameters were obtained from actual measurements or estimations, which also brought uncertainties.Therefore, robustness evaluations on the optimal schemes based on different mean values evaluation was carried out.New mean values of ANR j and APR j , which decreased by 10%, are shown in Formula ( 33) and ( 34): According to the new mean values, the relationship between system investment and constraint satisfaction probability under the new distribution was obtained (Figure 8).The result indicated that robustness of the optimal scheme was very good.Not only were the solutions still on the Pareto frontier, but also the marginal effect of system investment on the constraint satisfaction probability increased rapidly.The new relationship between system investment and constraint satisfaction probability under the new distribution is shown in Figure 9, and the result shows that solutions were still on the Pareto frontier, that is, the scheme was robust.However, it is inappropriate to be overly optimistic about a constraint satisfaction probability that was obtained in the optimal scheme.New mean values of ANR j and APR j that increased by 10% are shown in Formula ( 35) and ( 36): The new relationship between system investment and constraint satisfaction probability under the new distribution is shown in Figure 9, and the result shows that solutions were still on the Pareto frontier, that is, the scheme was robust.However, it is inappropriate to be overly optimistic about a constraint satisfaction probability that was obtained in the optimal scheme.The new relationship between system investment and constraint satisfaction probability under the new distribution is shown in Figure 9, and the result shows that solutions were still on the Pareto frontier, that is, the scheme was robust.However, it is inappropriate to be overly optimistic about a constraint satisfaction probability that was obtained in the optimal scheme.

Discussion
The multi-objective chance-constrained programming approach integrated with GA and robustness estimation has been provided and applied to a nutrient load reduction problem within the watershed of Plateau Lake.The methodology proposed in this article can present and deal with multi-objective and random information of the system (expressed in minimize system cost and maximize water qualification rate), evaluate robustness of the optimal solutions obtained by GA, and shows the following advantages for making decisions about optimal nutrient load reduction: (a) Multi-objective, chance-constrained, programming has the ability to trade off various optimal objectives in real cases, such as the minimal system cost and maximization of constraint

Discussion
The multi-objective chance-constrained programming approach integrated with GA and robustness estimation has been provided applied to a nutrient load reduction problem within the watershed of Plateau Lake.The methodology proposed in this article can present and deal with multi-objective and random information of the system (expressed in minimize system cost and maximize water qualification rate), evaluate robustness of the optimal solutions obtained by GA, and shows the following advantages for making decisions about optimal nutrient load reduction: (a) Multi-objective, chance-constrained, programming has the ability to trade off various optimal objectives in real cases, such as the minimal system cost and maximization of constraint satisfaction probability, and can represent and address random uncertainty (TN and TP reduction).(b) Robustness evaluation of optimal solutions can verify the robustness of the ideal schemes under variation in random parameters, which can affirm the optimal solution under deep uncertainty in the system.(c) The ideal solutions obtained from the model can provide a reasonable and reliable planning basis for decision-makers, which are screened under robustness estimates for all the optimal schemes that were obtained from multi-objective, chance-constrained, programming by using GA, even when multi-objective tradeoffs and random uncertainty exist simultaneously.
In addition, limitations still exist in the article.For example, the relationship between water quality and nutrient loads has been taken into account, but factors such as phytoplankton, light availability, and nutrients other than N and P, which may play a role in water quality, haven't been Water 2017, 9, 322 16 of 18 considered.In future research, influences from the factors above should be further considered and added into the model, so that the robustness of the methodology can be enhanced.

Conclusions
A multi-objective, chance-constrained, programming approach has been proposed by integrating MOP and CCP with GA and robustness estimation, which can balance the conflict between multiple optimal objectives and deal with random uncertainty within the system simultaneously.GA and robustness estimation were combined in multi-objective, chance-constrained, programming to acquire optimization schemes and to provide robust, optimal decisions.The model was applied to the optimal nutrient load reduction in the Lake Qilu watershed, which performed satisfactorily in balancing minimal system cost with maximization of constraint satisfaction probability, and in handling random uncertainty at the same time.
Ideal schemes obtained from this optimization model reflected the relationship between system investment and risk, and indicated that high system investment led to high satisfaction level with the constraints.Under the ideal schemes, detailed analysis was conducted in five sub-watersheds under six engineering measures on sewage treatment and reduction of TN and TP.Further, an evaluation of the robustness of the scheme was implemented, and the results demonstrated that the ideal schemes have superior robustness under uncertainty.This study is the first attempt to optimize reduction of watershed nutrient loads by using multi-objective, chance-constrained, programming, and it can also be applicable to other areas to manage the reduction of nutrient loads.

Figure 1 .
Figure 1.Location and watershed division of the study area.

Figure 1 .
Figure 1.Location and watershed division of the study area.
Generate an initial population under some coding scheme; (b) Convert the encoded individuals to decision variables under the corresponding decoding method, and calculate the adaptive value; (c) Set up a mating pool by selecting larger individuals from the population; (d) Develop a new generation after crossover and mutation; (e) Repeat steps 2-4 until the convergence criterion is satisfied.

Implication Interpretation i sub-watershed i = 1 :
Hongqi River sub-watershed i = 2: Zhewan River sub-watershed i = 3: Daxin River sub-watershed i = 4: Zhong River sub-watershed i = 5: Northest River sub-watershed j engineering measure j = 1: municipal wastewater treatment j = 2: municipal solid waste treatment j = 3: agricultural sewage collection and treatment j = 4: soil testing and formulated fertilization j = 5: constructed wetland j = 6: ecological restoration a i,j applicable factor for engineering measure j in sub-watershed i a = 1: appropriate a = 0: inappropriate ASC j Unit cost of measure j From "The 11th Five-Year Plan" and "The 12th Five-Year Plan" ANR j TN reduction of measure j Random parameter Water 2017, 9, 322 8 of 18

Figure 2 .
Figure 2. Relationship between system investment and constraint satisfaction probability.

Figure 2 .
Figure 2. Relationship between system investment and constraint satisfaction probability.

Figure 2 .
Figure 2. Relationship between system investment and constraint satisfaction probability.

Figure 3 .
Figure 3. Relationship between CR and constraint satisfaction probability for 19 Pareto solutions.Figure 3. Relationship between CR and constraint satisfaction probability for 19 Pareto solutions.

Figure 3 .
Figure 3. Relationship between CR and constraint satisfaction probability for 19 Pareto solutions.Figure 3. Relationship between CR and constraint satisfaction probability for 19 Pareto solutions.

Figure 4 .
Figure 4. Scale of sewage treatment measures in W1.

Figure 4 .
Figure 4. Scale of sewage treatment measures in W 1 .

Figure 5 .
Figure 5. Scale of sewage treatment measures in W2.

Figure 5 .
Figure 5. Scale of sewage treatment measures in W 2 .

Figure 6 .
Figure 6.Scale of sewage treatment measures in W2 and W5.

Figure 6 .
Figure 6.Scale of sewage treatment measures in W 2 and W 5 .

Figure 7 .
Figure 7. Relationship between system investment and constraint satisfaction probability when stochastic parameters obey uniform distribution.

Figure 7 .
Figure 7. Relationship between system investment and constraint satisfaction probability when stochastic parameters obey uniform distribution.

Figure 8 .
Figure 8. Relationship between system investment and constraint satisfaction probability when the mean values of stochastic parameters decrease by 10%.

Figure 8 .
Figure 8. Relationship between system investment and constraint satisfaction probability when the mean values of stochastic parameters decrease by 10%.

Figure 8 .
Figure 8. Relationship between system investment and constraint satisfaction probability when the mean values of stochastic parameters decrease by 10%.

Figure 9 .
Figure 9. Relationship between system investment and constraint satisfaction probability when the mean values of stochastic parameters increase by 10%.

Figure 9 .
Figure 9. Relationship between system investment and constraint satisfaction probability when the mean values of stochastic parameters increase by 10%.
i Ecological restoration target of sub-watershed i From "The 12th Five-Year Plan"

Table 2 .
Unit cost of the projects used in the model to optimize nutrient loads in the Lake Qilu watershed.

Table 3 .
Unit reduction of TN/TP for projects used in the model to optimize nutrient loads in the Lake Qilu watershed.

Table 4 .
Load and reduction target of the sub-watersheds in the model to optimize nutrient loads in the Lake Qilu watershed.

Table 5 .
Scale of measures for optimal schemes 12 to 19.