Improving the multi-objective evolutionary optimization algorithm for hydropower reservoir operations in the California Oroville–Thermalito complex

This study demonstrates the application of an improved Evolutionary optimization Algorithm (EA), titled Multi-Objective Complex Evolution Global Optimization Method with Principal Component Analysis and Crowding Distance Operator (MOSPD), for the hydropower reservoir operation of the Oroville e Ther- malito Complex (OTC) e a crucial head-water resource for the California State Water Project (SWP). In the OTC's water-hydropower joint management study, the nonlinearity of hydropower generation and the reservoir's water elevation e storage relationship are explicitly formulated by polynomial function in order to closely match realistic situations and reduce linearization approximation errors. Comparison among different curve- ﬁ tting methods is conducted to understand the impact of the simpli ﬁ cation of reservoir topography. In the optimization algorithm development, techniques of crowding distance and principal component analysis are implemented to improve the diversity and convergence of the optimal solutions towards and along the Pareto optimal set in the objective space. A comparative evaluation among the new algorithm MOSPD, the original Multi-Objective Complex Evolution Global Optimization Method (MOCOM), the Multi-Objective Differential Evolution method (MODE), the Multi-Objective Ge- netic Algorithm (MOGA), the Multi-Objective Simulated Annealing approach (MOSA), and the Multi-Objective Particle Swarm Optimization scheme (MOPSO) is conducted using the benchmark functions. The results show that best the MOSPD algorithm demonstrated the best and most consistent performance when compared with other algorithms on the test problems. The newly developed algorithm (MOSPD) is further applied to the OTC reservoir releasing problem during the snow melting season in 1998 (wet year), 2000 (normal year) and 2001 (dry year), in which the more spreading and converged non-dominated solutions of MOSPD provide decision makers with better operational alternatives for effectively and ef ﬁ ciently managing the OTC reservoirs in response to the different climates, especially drought, which has become more and more severe and frequent in California.


Introduction
In the last few decades, studies in the field of Evolutionary Algorithms (EAs) in water resources have focused mainly on incorporating new strategies in development of new algorithms or on verifying the suitability of a particular algorithm in solving different kinds of conceptual water-resources management problems (Maier et al., 2014). Tremendous efforts have been made to develop EAs with new methodologies, which are inspired by various natural phenomena, such as foraging behavior of ants (Dorigo et al., 1996) and bees (Pham et al., 2011) and social behavior of fish colonies (Kennedy et al., 2001). These innovations lead to a large number of studies that apply EAs to many fields of water-resources management problems (Nicklow et al., 2010), including ground-water calibration, water treatment and reservoir operation. For reservoir operation problems in particular, EAs are recognized as good decision-making support tools because of their multiple advantages. In general, according to Abraham and Jain (2005), there are many advantages of using EAs because they require little knowledge about the problem being solved, are less susceptible to the shape or continuity of the Pareto front, are easy to implement, robust, and could be implemented in a parallel environment. The performance of EAs on human-designed test functions with various mathematical properties, such as discontinuous, nondifferentiable, non-convex, and multimodal Kumar 2006, Reed et al., 2013), have been reported in many studies. The convenience of using EAs to find the Pareto optimums, and their effectiveness for highly complex test functions provide some confidence that EAs will also be successful with complex real-world problems. In addition, EAs have proved to be effective and suitable, especially for solving multi-objective problems. In multiobjective optimization problems, EAs consider all of the objectives simultaneously without any user defined weights to each objective, and the population-based searching mechanism enables EAs to generate a set of equally important solutions, termed nondominated solutions (Deb, 2001), in a single run instead of performing a series of separate runs (Abraham and Jain, 2005). The non-dominated solution set forms a Pareto front which is able to provide decision makers with trade-off information between conflicting objectives. These benefits mentioned above allow EAs to provide decision makers with a reliable set of solutions to realworld problems, as well as the confidence about the use of these solutions that are able to consider the often multi-objective nature of decision making in reservoir operations. Therefore, as Maier et al. (2014) concluded, EAs have bolstered our ability to solve problems that are more relevant to real-world systems. Many studies have applied EAs to realistic reservoir management cases, especially to reservoir release scheduling and hydropower scheduling problems (Reddy and Kumar 2006, Reddy and Nagesh Kumar 2007, Afshar et al., 2007Chen et al., 2007;Farmani et al., 2005;Chang and Chang, 2009;Wardlaw and Sharif, 1999;Cheng et al., 2008).
Although EAs have gained much popularity in the field of algorithm development and theoretical application to reservoir operations, there are few reports about the actual uses of EAs on realistic reservoir-operation problems. As Jones et al. (2002) and Oliveira and Loucks (1997) concluded, there are few documented applications in which decision makers actually use EAs. Shepherd and Ortolano (1996) reported on personal communications with system operators, stating that they "don't like being told what to do" and "want to do it in his own way". Yeh and Becker (1982), Yeh (1985), Wurbs (1991Wurbs ( , 1993 summarized that there is a gap between theoretical development and actual real-world implementation. In the state-of-the art review of optimal operation of multi-reservoir systems, Labadie (2004) concluded that one reason for the gap is that many reservoir operators lack the confidence in simplistic problem formulation, which purports decision makers to replace their judgment and prescribe solution strategies. To fill the gap between theoretical development and realistic problem implementations, Labadie (2004) and Goulter (1992) suggested more interactions and involvement of decision makers in the problem formulation as well as better suitable packaging of the optimization approaches along with designed problems. U.S. Army corps of engineering (USACE, 1990) stated that, among all of the applications of optimization techniques, the combined use of simulation models and optimization has been found to be an effective analysis strategy for reservoir operation problems. Maier et al. (2014) also outlined the current status and future research challenges and directions in the development of a fundamental understanding of both realworld problem and algorithms. These suggestions indicate that both problem formulation and algorithm development are equally important to promote the uptakes of EA by decision makers.
In order to enable EAs to be more broadly applied to realistic water resources problems, this paper will focus on applying an improved EA to a conceptual reservoir operation model that is built on realistic objectives and constraints for the head-water supply region in California. The simplification of reservoir topography using different methods is analyzed so that the impact of simplifications and assumptions can be better understood.
First, we focus on developing a water-hydropower joint reservoir management framework of the OrovilleeThermalito Complex (OTC) in the California's State Water Project (SWP). We show that, by explicitly formulating the non-linear relationship between reservoir storage volume and water elevation, the use of the optimal polynomial method to fit the reservoir topography is able to extensively reduce the approximation errors compared to the traditional piecewise linearization method. The piecewise linearization method is commonly used to simplify these non-linear problems in reservoir-system planning and operation models, such as that in CALSIM (Draper et al., 2004), CALSIMii (Ferreira et al., 2005) from DWR and HEC-5 from the U.S. Army Corps of Engineers (USACE, 1998). These models implement piecewise linearization for its computation efficiency and programming simplicity. However, Labadie (2004) concluded that the non-linear challenges in optimal reservoir-system management, especially with inclusion of hydropower generation as objectives or constraints, should be addressed directly by non-linear programming as well as non-linear optimization techniques. Bay on et al. (2009) suggested that linear simplification of the storageeelevation (SeE) curve can produce serious errors in the optimal solution. With respect to these concerns, in this paper, we experiment with a non-linear approach to OTC's problem formulation, in which electrical generation and consumption are formulated as one of the objectives, while the reservoir topography is approximated with optimal order polynomial function. Two SeE curve-fitting methods, including linear simplification and optimal polynomial are compared with successive parabolic interpolation (Siddall, 1982;Kharab and Guenther, 2011), which is a technique for reconstructing a continuous unimodal function by successively using the parabola function (polynomials of degree of two) to fit the measurement points. Because the true SeE relationship is unknown, we assume that the SeE curve generated by successive parabolic interpolation is the "true" value.
In addition, we present the enhancements made to the Multi-Objective Shuffled Complex Evolution Global Optimization Algorithm (MOCOM-UA) using crowding distance and principal component analysis. The MOCOM-UA algorithm  is an extension of the successful single-objective Shuffled Complex Evolution (SCE-UA) global optimization algorithm (Duan et al., 1992). It combines the strength of competitive evolution (Holland, 1975), the Nelder-Mead (simplex) method (Nelder and Mead, 1965), and Pareto ranking (Goldberg and Holland, 1988). The MOCOM uses a triangular possibility distribution to assign the possibilities to the members in the population so that parents that have better objective function values are more likely to be chosen to produce offspring. This strategy is able to generate a fairly uniformed non-dominated solution set. Various studies have tested its usefulness in hydrologic models Gupta et al., 1998Gupta et al., , 1999Boyle et al., 2000;Xia et al., 2002;Leplastrier et al., 2002). However, Gupta et al. (2003) pointed out several weaknesses of MOCOM-UA, in which two major issues are its clustering tendency of non-dominated solutions and premature phenomena in certain cases. Vrugt et al. (2003) concluded that the failing is due to its fitness assignment of the Pareto ranking, in which members having an identical rank are not distinguished when assigning a selection possibility. The second issue is due primarily to the high dimensionality of the optimization problem (the so-called "curse of dimensionality"), which prevents the evolution of the population from exploiting the entire decision variable space. In response to these concerns, we introduce two modules into the MOCOM. The first module revises the selection possibilities of the members with identical Pareto ranking so that the generated non-dominated solutions can form a more uniformed distribution along the Pareto front. The second module monitors the diversity of the population during evolution based on principal component analysis, which has been shown to prevent the population from degenerating.
Finally, the algorithm is applied to the OTC's reservoir-operation model to generate better operation alternatives than the real operation in response to the drought climates and other watersupply conditions as compared to the MOCOM-UA. The generated non-dominated solutions from MOSPD are able to discover the limitation of the system, determine the optimal operations strategies, and detail the sensitive days, when reservoir operators are making decisions on reservoir-release amounts to maximize their specific objectives (e.g., hydropower generation or storage volume).
Because climate variability and its impacts on regional water resources is the leading cause for uncertainty in the waterdistribution system, our goal is to provide improved solutions, based on the reservoir operation's objectives under varying historical climate conditions. The application of MOSPD on OTC's reservoir operation problem provides valuable guidance for making improved and enhanced decisions with this climate variability. Specifically, certain solutions we demonstrate show an increase of the system's storage volume, which promotes the OTC's capabilities to prepare certain special operations, such as "water supply loan" and "water consolidations" strategies (Kelly, 1986), in response to dry conditions. The other solutions that increase hydropower generation would be better operational alternatives for wet year and average year because the extra clean energy from hydro-power plants supplements the annual energy consumption from the SWP's water transferring and pumping sectors.
The contribution of this paper focuses on narrowing the gap between theoretical development and actual real-world implementation of EAs. In detail, we build a conceptual model based on the survey of realistic reservoir functionalities, system configuration and the decision maker's goals in operating the OTC's facilities. As Maier et al. (2014) summarized, it is important to understand the impacts of simplifications and assumptions in the way real problems are represented mathematically. The comparison study of the different curve-fitting methods helps to quantify the errors introduced by simplifying the physical relationship of reservoir topography and its impact on optimal reservoir operation strategies. In addition, we apply two popular techniques to the original MOCOM algorithm to promote the effectiveness and accuracy of the algorithm, as well as to strengthen the EAs capability in solving realistic problems. Comparison studies among several state-of-art EAs on human-designed test problems are carried out, and the application to the OTC's problem is provided for assisting the decision making for OTC. This paper is organized in six sections. Section 2 explains the OTC's configuration, current operation, and the explicit formulation to reduce simulation errors. Section 3 introduces the optimization algorithm and the two enhancement modules. Verification results are also included. In Section 4, we implement the new optimization algorithm into OTC's problem. Section 5 presents the conclusion, and some limitations and future works are listed in Section 6.

Water and hydropower's sustainability in OTC
The main functionalities of OTC are to supply fresh water and generate hydro-electricity for the California's SWP. The SWP (Fig. 1) is the nation's largest state-built water and power development and conveyance system (DWR, 2013a) operated by the California Department of Water Resources (DWR), with its main purpose to store and supply water from the precipitation-concentrated northern area of to the water-scarce central and southern regions. Currently, SWP has conveyed an annual average of 3.577 Â 10 9 cubic meters (m 3 ) of water with the potential of providing 5.181 Â 10 9 cubic meters (m 3 ) designated water allocation to its users. The main concern for the future SWP project lies in its ability to make efficient water release to meet increasing water demands. As DWR (1993) projected, a net, statewide water-supply deficit of 4.07 Â 10 9 to 5.181 Â10 9 m 3 by 2020 is expected, which implies that SWP is under a significantly severe burden to supply water under potential drought condition. Calendar year 2013 closed as the driest year in recorded history for many areas of California, and current conditions suggest no change in sight for 2014. On January 31, 2014, the DWR announced several actions to protect the health and safety of Californians from the effects of more severe water shortages. Those actions include lowering the anticipated allocation of water to customers of the SWP from 5% to zero, which marks the first zero allocation announcement for all customers of the SWP in its 54year history. In order to better meet the water shortage and response to varying water-supply conditions, a more efficient operation of SWP's water storage facilities is required so that the SWP's water supply can be more stable and sustainable.
The OrovilleeThermalito Complex (OTC), which is located in northern California in the foothill of the Sierra Nevada Mountains, consists of several reservoirs and hydropower plants in and around the city of Oroville in Butte County. As the most vital fresh headwater supply source and power development for SWP, OTC delivers about 4.317 Â 10 9 m 3 of water at maximum capacity and generates more than 2.8 billion kilowatt-hours of power annually. When water is needed in SWP, the OTC releases water into the Feather River through the Oroville Dam, Thermalito Diversion Dam, and the Thermalito Afterbay. The released water travels down the river to where the river converges with the Sacramento River and continues to be pumped or diverted to southern and central California for various demands. These processes are carried out by jointly operating OTC's 10 major facilities: (1) three cascade reservoirs: Lake Oroville, the Thermalito Forebay, and the Thermalito Afterbay; (2) three hydro-electrical powerplants: the Hyatt Powerplant, the Thermalito Diversion Powerplant, and the Thermalito Pumping-Generating Plant; and (3) five other facilities, including the Thermalito Diversion Dam, the Lake Oroville Dam, the Feather River Fish Hatchery, the Fish Barrier Pool, and the Thermalito Power Canal.
To increase the water's sustainability in OTC, we consider the accumulated daily-storage volume during one month as the first objective. This objective indicates the regional "water supply loaning" availability and capability to consolidate storages in response to emergent drought conditions. The concept of "water supply loaning" is a special reservoir operation scheme during drought conditions. According to Kelly (1986), "water supply loaning" is recognized as an efficient special reservoir operation during drought conditions, which temporarily transfers large amounts of water from one water-supply system to another waterscarce system in a short period of time (days or weeks) to mitigate the drought condition. This action requires the loading system to have sufficient and continuous water storage so that the transferring of water does not jeopardize the loading system's engineering constraints and demand constraints. Historically, the "water supply loaning" operation has been carried out between SWP and California's Central Valley Project (CVP) since August 1977, during which 9.25 m 3 of water has been transferred from the OTC and the SWP's reservoirs to the San Luis Reservoir to temporarily meet the CVP irrigation demand. This action is completed in 4 days. During that time the SWP stores relatively high levels of water storage in its reservoir system, while the CVP is facing a forecasted water shortage caused by the 1976e1977 California droughts. Consolidation of water storage is another drought-response operation scheme that merges the water storage of several reservoirs into one as quickly as possible to decrease the evaporation, seepage and in-stream losses. The consolidation of water storage is extended through the Orland Project, a project within the CVP, during the same period of 1976e1977. A necessity of these two special operations is that water must to be transferred from the region with relatively higher storage volumes in a short period of time (days or weeks) to the region with low storage volumes. Therefore, the accumulated daily storage volume during a particular month is able to measure the capability to quickly enable these special operations.
In order to increase OTC's power-supply stability, we set the net electricity generation during a given month as the second objective. California's SWP is the largest single electricity user in the state, accounting for 3% of all electricity consumed statewide (DWR, 2013b). Most of the energy is used to deliver water to the southern California region, where pumping 1 m 3 water through the Tehachapi Mountain consumes about 2.432 kWh electricity. Annually, about two-thirds of SWP's power comes from its hydro-electrical plants, and the remainder is supplemented by the coalfired Reid Gardner plant in Nevada and power purchases and exchange programs (DWR, 2011(DWR, , 2012(DWR, , 2013c. Thus, increasing OTC's total net hydro-power generation could help SWP build its electrical power development towards cleaner and more self-sufficient levels. More important, from the Governor of California's perspective, these two objectives fit into the future focus of managing water and energy. California, along with other places around the world, has already recognized the strongly connected interaction between water and energy, termed as water and energy nexus (DWR, 2013b;CEC, 2005;Cohen et al., 2004, U.S.D.E., 2006, Wilkinson, 2000. The two objectives we considered in the OTC complex are based on the functionalities and the purposes of the reservoirs in the complex. In other literature, there are significant ways to include hydropower generation in reservoir operation, such as considering short-term hydro scheduling, annual production, and trade-offs with agricultural water uses (Cheng et al., 2014;Hassanzadeh et al., 2014;Li et al., 2013;Gil et al., 2003) In these ways, hydropower generation is either considered as a single objective in complex reservoir systems or measured by economic value to other water uses. However, we set the objectives based on the realistic functionalities of the regional reservoir system. The total storage in the complex indicates the water supply potential for California's SWP, while the hydropower generation ensures that the water released from this complex will be delivered through the whole California's water distribution system. The two objectives are essential to the operation of OTC and SWP.
Besides the two major purposes (water supply and power generation), current operation in OTC has other goals, such as flood control, temperature control for fish habitat, and water-quality monitoring. The priorities for these goals vary from month to month. First, OTC's facilities are operated under flood-control requirements specified by the U.S. Army Corps of Engineers (USACE), which requires that the OTC's reservoirs to keep a flood-control space during the peak flood season (October 15eMarch 31). The required flood-control space varies between 462.55 Â 10 6 and 925.11 Â 10 6 cubic meters, depending on the accumulated precipitation parameter prescribed in the flood-control manual (DWR, 2013d). In addition, OTC is required to maintain a specific outlet temperature for salmon and steelhead trout spawning conditions during the whole year. A temperature range of plus or minus 4 F is allowed (DWR, 2006) only from April through November. Moreover, water-quality standards are designed to meet several water quality objectives such as salinity, Delta outflow, river flows, and export limits (DWR, 2006).

Mathematical formulation
In Fig. 2, the facilities and key infrastructures for the OCT's reservoir system are presented. In order to mathematically formulate OTC's reservoir system, the main components are described separately. Reservoirs, dams and power plants are represented by node numbers 1e6. Water flows are illustrated by blue arrows (in the web version), and the interactions between power facilities and electricity grids are shown with red arrows. OTC has two major outlets to the Feather River. One is located at the Thermalito Diversion Dam, which belongs to the Thermalito Forebay area, and the other is connected to the Thermalito Afterbay. Other water flows are upstream inflows and the regulated daily flows, such as deliveries to nearby cities, power flows, and pumpback flows. In this system, we consider the Feather River daily releases to be optimized because of the regulations as mentioned above. Detailed information for the fixed water flows is listed in the Appendix (Fig. A1 and Table A1).
The two main types of constraints that relate to mean sea-level water elevation are included in Fig. 2 and are the reservoir waterelevation constraints and the normal static water head constraints of power plants (DWR, 2013a). The reservoir waterelevation constraints are described as the upper and lower bounds, which represent the reservoirs' physical water capacities based on the engineering design of the dams. The normal static head upper and lower bounds ensure the safety operation requirements for the power plant's generator and other components. The OTC's problem can be expressed as: Optimizing decision variables U t in order to Max 8 > > > < > > > : which is subject to: G t;min G t G t;max (8) P t;min P t P t;max (9) In Eqs.
(1)e (10), T is the number of days in a month; N is the number of storage facilities; M is the number of power facilities; S t is the storage at time step t; I t is the reservoir inflow at time step t; O t is the reservoir outflow which consists of regulated power discharge term Q t , Feather River release term U t , evaporation and other losses E t , and water deliveries to local urban areas D t at each time step t; G t is the power-generation capacity at time step t; P t is the pump-back electricity capacity at time step t; h and m are the electricity generation and pumping efficiency, respectively; Q 0 t is the regulated pump-back flow; and H t is the normal static waterhead difference between the upper and lower reservoirs. In Eqs. (5)e(10), lower bounds (upper bounds) of the constraints are noted with "min" ("max") as subscripts.
Note that in Eq. (4), normal static water-head difference is determined by the upper and lower reservoir water elevation h upper and h lower , each of which is a function of the reservoir's storage volume. The relationship between reservoir's storage volume and elevation is non-linear due to the irregular shape of reservoir topography. Here, we use function f and g to represent this relationship.
In Eqs.
(1)e(10), there are two crucial non-linear aspects that could influence decision maker's choice and the optimized solution. First, according to Eq. (1), the net electricity generation changes non-linearly with the discharge. The electricity generation is calculated by the discharge Q t multiplied by the water-level difference H t . More discharge generates more hydro-electricity and also results in less storage volume in a given reservoir. The less storage volume leads to a lower water level, according to the nonlinear functions f and g. Eventually, the increasing discharge also decreases the net electricity generation amount. Second, in this process, reservoir topography functions f and g play an important role in determining the actual changing value of normal water static head difference H t . In the model framework, we choose the time step to be daily for the mid-range (month or seasonal) optimal reservoir operation. By using a shorter the time step, such as hourly, the formulation is more realistic to the dynamic hydropower generating and pumping mechanics, because the turbine efficiency curve and other dynamic process are able to be included (Diniz et al., 2007). However, the hourly operation requires quick opening and closing the releasing vaults or gates, and the failure to ensure this could cause operation difficulties and quicken the aging of facilities. Moreover, with regard to the optimal reservoir watersupply objective, the operation on a daily scale is a relatively good tuning scale for mid-range optimal water-supply planning. Therefore, we use daily as the time step and follow the similar manner that used by Li et al. (2013) to estimate hydropower generation and pumping in cascade reservoir system for long term operation.
To present the relationship between storage volume and water elevation, a piece-wise storageeelevation curve is typically used. As mentioned earlier, Bay on et al. (2009) argued that linear approximation of the reservoir storageeelevation (SeE) curve could result in serious errors in calculating hydro power. Here, we use a nonlinear approach, in which the SeE curves for Oroville Lake, Thermalito Forebay, and Thermalito Afterbay are fitted with 8-, 13-and 14-order polynomial functions (Fig. 3(a)e(c)). The orders of fitting functions are chosen based on error variability between observations and fitted values. As shown in Fig. 3(d)e(f), with the increasing polynomial orders, the sum of the squared residual is decreasing. The polynomial function can generate a near-perfect fit by increasing the number of orders. However, higher order of polynomial could result in requirements for more complex computation and difficulties for optimization algorithm to find local minimums. Here, we choose the order of fitting function based on the fact that the sum of squared residual is relatively small, but it does not decrease dramatically with increasing orders, which are 8, 13, and 14 for the Oroville Lake, the Thermalito Forebay and the Thermalito Afterbay, respectively. The sample SeE measurements are retrieved from the SWP construction reports (DWR, 1974) and the monthly reports published by SWP operations control office.

Building an enhanced multi-objective optimization algorithm
In this section, we improve a Multi-Objective Complex Evolution (MOCOM-UA) global-optimization method by introducing two new modules in order to enhance the diversity and convergence performance of the non-dominated solutions.

Diversity and convergence
Diversity and convergence (Deb, 2001) are the two major concerns to evaluate whether a set of non-dominated solutions is beyond another set. The diversity refers to the non-dominated solutions set's coverage along the Pareto front. And the convergence measures the non-dominated solutions set's closeness towards the Pareto front. The ideal case is when the non-dominated solutions are identical (i.e., exact overlap) to the global Pareto front.
A solution set with more diverse members is preferred by decision makers than the one in which its members are relatively similar. The diversity of the operation alternatives is able to give decision makers various options in response to the changing climate conditions. In dry conditions, various operation alternatives could help decision makers to revise their operations to conserve certain amounts of reservoir storage for future water supply. In wet scenario, diversified solutions are able to provide decision makers with efficiently water releasing strategies for maximizing other non-water-supply objectives, such as hydropower generation.
Equally important, convergence of a non-dominated solution set helps decision makers understand the limitations of the system as well as the range of potential benefits. A non-dominated solution set which fails to approach the higher values of objective function values is not likely to be used by decision makers. In the decision makers' perspectives, greater benefits can be gained by using a more converged non-dominated solution set. In dry scenarios, the failure to reach a converged solution, which maximizes the sustainable water supply, means the water is not used efficiently in response to potential drought. In addition, the solution that is closer to the global Pareto front could provide greater benefit for both water and power supply in wet scenarios, when the supply is not as critical as in drought conditions.

MOSPD algorithm
To address these issues, we update the MOCOM algorithm  with two enhancement techniques, and entitle the new version Multi-Objective Shuffled Complex Evolution with Principal Component Analysis and Crowding Distance Operator (MOSPD). We add two distinct modules to the original MOCOM algorithm (Fig. 4). The first module is called the "possibility-adjustment" module (Fig. 5) while the second is called the "dimension monitoring and restoring" module (Fig. 6).
The general steps of MOCOM can be summarized as follows according to the flowchart in Fig. 4: (1) a total of m Â p points are randomly sampled in the parameter space to form the initial population, where m is the number of complexes and p is the total number of individuals in a complex; (2) the functions are evaluated for each individual; (3) the entire population is shuffled and split into m groups (complexes). In each of the groups (complexes), the p individuals form the sub-population; (4) the Pareto ranks (Goldberg and Holland, 1988) are calculated for the entire  population; (5) a triangular possibility function is used to assign a selection possibility to each individual according to its Pareto ranks; (6) a simplex is constructed by selecting n þ 1 individuals according to the possibility distribution of the sub-population derived from the previous step; (7) the Nelder-mead evolution strategy (Nelder and Mead, 1965) is implemented to obtain a new individual, and the population is updated; (8) the steps from (3) to (7) are repeated until the maximum of the Pareto ranks in step (4) becomes 1, which means the individuals in the population are all non-dominated in relation to each other. The details of the "possibility-adjustment" module are exhibited in Fig. 5. As Figs. 4 and 5 show, this module is embedded in the main routine of the MOCOM. The steps for this module can be summarized as follow: (1) the individuals with a Pareto rank of 1 are selected and stored in a temporary set. For better illustration, the total number of individuals in this temporary set is t; (2) the crowding distance vector D(d1, d2, …, dt) is calculated according to the Euler distance between an individual and its neighbors in the objective function space (Deb, 2001); (3) the selection possibility of each individual in this temporary set is calculated as where the P i is the selection possibility from the main routine that is calculated using the triangular possibility density function; (4) the selection possibility is updated from P i to P l for the individuals in the temporary set; (5) the process in steps (2)e(4) is repeated for the individuals with Pareto ranks that equal to 2, 3 … and so on, until the maximum Pareto rank is reached. This module is intended to ensure that the members with the same Pareto ranking are put into the same group as well as their selection possibilities. In each group, the crowding distance (Deb, 2001) is calculated for all members according to the distance between two of the closest neighboring members. The crowding distance technique has been successfully tested and applied as an enhancement to evolutionary algorithms Nagesh Kumar 2007, Azadnia andZahraie, 2010). A detailed example of the crowding distance calculation is in the Appendix (Fig. A2 and Table A2). Then, a new selection possibility for each member in this group is computed, which equals the total selection possibility, which is assigned by the main routine of MOCOM multiplied by a distance coefficient. The distance coefficient is calculated as each member's crowding distance divided by the total crowding distance of all members in this group. The new selection possibility replaces the original one for each member in this group. This adjustment is looped from the group with the lowest Pareto ranking to the one with the highest Pareto ranking until all members in the population are assigned with a new selection possibility. Different from the MOCOM, in this new strategy, members with identical Pareto rankings are assigned different possibilities, with the criterion that a member which locates remotely from its neighboring members is more likely to be selected than those closely clustered in a limited objective space.
The second module, "dimension monitoring and restoring", is shown in Fig. 6. The aim of this module is to capture and restore the lost dimension during evolution. The module uses the Principal Component Analysis (PCA), which was invented by Pearson (1901) and further developed by Hotelling (1933Hotelling ( , 1936. The PCA is a statistical procedure that transforms a given dataset into an orthogonal coordinate system, in which the first coordinate e termed first principal component of PC has the largest variance of the projection from the data set and other coordinates have smaller values of variances in descending order. Sometime the lower-ranked PCs have negligible variance, which means that the data set has a dimensionality reduction. In the MOSPD algorithm, once a lost dimension is discovered by PCA, a new point is sampled along the corresponding PC, and the member that results in the dimension lost is replaced with this new point. The lost dimension here is defined as a dimension has an eigenvalue less than 1% of the summation of all eigenvalues of the covariance matrix. This module can be generalized in terms of two steps. The first step is to check the dimensionality of the space spanned by all individuals in the population using the following procedures: (a) Let C ¼ ½c ij ¼ ½x 1 …x mp be the matrix with the coordinates of each point as its columns. Then, C has dimensions of n Â (m Â p), where, x i ; i ¼ 1; 2…mp are the points in the population, n is the dimensionality of the problem, m is the number of complexes, and p is the number of individuals in a complex, as mentioned above; (b) Transform the original coordinate system to a normalized coordinate system by centering and normalizing each row of C and obtain C 0     (c) Calculate the covariance matrix of C ' and denote it as R.
Obtain the eigenvector and eigenvalues of R. Each eigenvector is a principal component (PC) of the population, and its corresponding eigenvalue measures the variance of the individuals along the direction of that PC. (d) If the variance along one PC is too small, it means the individuals are not spanning well over that direction, and we define it as a lost dimension. Mathematically, a threshold of 1% of the total variance as a lost dimension is used to measure whether a variance is too small.
Once a lost dimension is detected, the second step is to restore the lost dimension by randomly sampling a new point along the PC: (a) Sample a point from the side of centroid of C ' along the PC: where c 0 is the centroid of C 0 , a is random number from a normal distribution with mean ¼ 2 and variance ¼ 1, r is the radius of the entire population in C 0 , and l is the unit vector along the PC.
(b) Transform x 0 back to the original coordinates and evaluate the objective function. Then, update the population and terminate the module.
This method has already proved efficient and effective in solving population degeneration problem in high-dimension, singleobjective optimization problems (Chu et al., 2010(Chu et al., , 2011 and reducing the number of objectives in multi-objective optimization problems . However, the implementation of this method in the MOCOM algorithm is rarely reported.

Population size
The population sizes for the selected algorithms are identical for each test function. For SCH, POL, and FON, the population size is set to 16, while for KUR, ZDT1, ZDT2, ZDT3, and ZDT4, a population size of 124 is used. The results, along with simulated global Pareto front are shown in Fig. 7. For MOSPD and MOCOM, the population at selected number of function evaluation during the evolution is shown in Fig. 8. Detailed information on these benchmark functions is included in the Appendix (Table A4). Note that it is possible that algorithm performance could vary as the population size changes; nevertheless, we use an identical population size for each algorithm in order to conduct a fair comparison. Detailed population sizes for each test problem can be found in the Appendix Table A4, last column.

Other settings
The key settings for each algorithm are listed below: For the MODE, the crossover constant is set to 0.9 and the scalar factor is 0.45; For the MOGA, the crossover probability is 0.5, and the mutation probability is 0.1; For MOSA, the cooling factor is set to 0.87; reheating temperature is set to 5. For the MOPSO, the cognitive learning factor is 2, the social learning factor is 2, the movement velocity is 0.5, the inertial constant is 0.5, and the maximum number of individuals in each particle is set to 5. All test runs are set a maximum of 50,000 function evaluations.
To evaluate the performance of a multi-objective optimization algorithm, Deb (2001) suggested at least two metrics should be used. One metric measures how close is the non-dominated solution set towards the global Pareto front, and the other evaluates the spreading extent of the non-dominated solution set along the global Pareto front. Here we choose diversity metric D and generational distance (GD), which are well-established performancemeasurement indices from Zitzler and Thiele (1999) and Deb (2001. The formula to calculate D and GD and are listed in the Appendix in Eqs. (A1) and (A2).
The range of GD is greater or equal to zero, and D lies within the range of [0, 1]. A smaller value of GD indicates that the nondominated solutions have better convergence towards the global Pareto front. Similarly, when the D value is closer to zero, a more diverse spread of the non-dominated solutions along the global Pareto front is denoted. We calculate these two metrics to show the difference between the MOSPD and the original MOCOM algorithm.
Tables 1 and 2 list the spread metric D and GD values for each of the algorithms on the six test problems.
According to the test results, MOSPD shows better capability of expanding non-dominated solutions along the Pareto front. According to Table 1,  Based on the result of limited sensitivity test shown in this study, MOSPD has exhibited superior performance over MOCOM in all cases and better diversity measurements over the other four algorithms, especially for high dimension problems. In some of the test functions (Fig. 7(d)e(g)), the non-dominated solutions with MOCOM tend to cluster in a fairly small region in the objective space, while the solutions with MOSPD have a likely uniformed spread. As Table 2 shows, excluding the SCH and KUR cases, the smaller values of GD indicate that the MOSPD algorithm is also able to generate more converged non-dominated solutions than those generated from other algorithms. However, the GD values for SCH and KUR with MOSPD are still very competitive when compared to others. For the functions POL, ZDT1, ZDT2, and ZDT4 ( Fig. 7(b), (e),(f), and (h)), MODE, MOGA, MOSA, MOPSO, and MOCOM fail to discover the global Pareto front in the objective function space and the search stops at a local optimal with higher values for both of the objectives. In contrast, the MOSPD is able to escape from the local attractions and reach the objective function values that are very close to the global optimums. When compared with the MOCOM and MOSPD evolution processes (Figs. 8 and 9), MOCOM fails to maintain the population diversity during the evolution (Figs. 8(g), 9(a) and (c), (e), and (g)).
Although MOCOM converges more quickly toward lower objective function values for some cases during the evolution ( Fig. 9(a) and (c)), and the expansion of the population is poor. For most of the cases, the population evolution of MOCOM eventually stops with higher objective function values due to the large local minimum attraction, while the population evolution of MOSPD in the same tests (Figs. 8(h),9(b) and (d), (f), and (h)), are likely to maintain the population diversity, and the final population is able to reach a more expanded location with lower objective function values.
The improved performance of MOSPD is due to the two newly introduced modules. In the original MOCOM, the evolution process is guaranteed to be competitive based on the selection criteria that the possibility that "better" parents contribute to the generation of offspring is higher than that of "worse" parents. However, this strategy does not guarantee that the offspring can uniformly locate along a certain Pareto front. In MOCOM, according to the identical Pareto ranking, the parents, which can contribute to generating an offspring with a better location that is towards a sparse surrounding area, are equally treated as the parents that are only able to produce the offspring with a worse location. The new strategy used in MOSPD remedies the equal treatment criteria by adjusting the selection possibility using the crowding distance measurement. The crowding distance-based possibility selection strategy ensures that the parents that produce offspring with a better location are assigned with a higher chance to be selected. This is a more robust strategy because more diversified offspring help to form the final non-dominated solutions towards a uniform distribution along the global Pareto front. Similarly, for some of the cases (Fig. 8(d)e(f)), the original MOCOM cannot escape the local attractions due to the fact that the decision variables in certain dimensions happen to be the same or very close values, which causes the population to lose the ability to continue searching decision variable space in this dimension. Eventually, the members with lost dimensions will be stuck at the local Pareto optimal with higher values on both objectives than the global Pareto optimal. MOSPD overcomes this problem by repeatedly restoring the lost dimensions and preserving the population's vitality of searching larger parameter spaces.
The newly developed MOSPD algorithm combines (1) the strengths of the MOCOM algorithm  (2) the concept of the crowding distance-based offspring selection probability strategy (Deb, 2001); and (3) the tool of PCA (Hotelling, 1933(Hotelling, , 1936) that restores and maintains the population diversity during searching. According to the test results, MOSPD is a more efficient and effective algorithm over MOCOM regarding convergence and diversity of the non-dominated solutions. The improved algorithm (MOSPD) theoretically is able to provide more diversified and accurate alternatives as a better decision-making tool on the OTC's reservoir daily release-operation problem.

Application
In this section, we will implement the improved MOSPD into the OTC's reservoir daily release problem to generate operational alternatives based on different climate conditions. We also analyze the difference between two extreme solutions, in which one maximizes storage and one maximizes net electricity generation. Their potential benefits in responding to these wet, average, and dry conditions are analyzed as well.

Settings
In the OTC's problem, we conduct the simulation with both the MOSPD and MOCOM algorithms with identical settings: (1) The tunable parameters are the Thermalito Forebay and Afterbay daily Feather River release amounts and the remaining internal flow amounts within the OTC system are set to the realistic operation values. As mentioned in Section 2, there is a constraint on the monthly total release amount for the Thermalito Forebay and the Afterbay. Thus, the tunable parameters have a dimension of 2 Â (number of days in one month À 1). Other constraints are (1) the storage capacity constraints for the upper and lower limits as shown in Eq. (5); (2) power generation capacity constraints as shown in Eq. (8); (3) pump-back electricity capacity constraints as shown in Eq. (9); and (4) the static water-head constraints represented in Eq. (10).
(2) The boundary-handling method (referred to as the reflecting method) is intended to reset any newly generated offspring that violates its respective constraint. During the evolution, the boundary acts as a mirror and reflects the projection of the displacement. Then, the displacement adjusts the offspring's location in the parameter space.
(3) Objective functions are the daily storage volume total and net electricity generation as shown in Eq. (1). As mentioned in Section 2, the first objective is an important factor for initiating special operations during drought conditions, and the second objective supplements the energy shortage for transferring and pumping water in the SWP.
(4) The simulations are carried out for the period of AprileJune, which is snow-melt season for the Sierra Nevada Mountain. Three different years (1998, 2000 and 2001) are included because, according to the SWP water supply office, these three years are officially recognized as the typical wet, average and dry year respectively.

Results
To demonstrate the accuracy of the joint model, we carry out an objective function value comparison between the real operations scenarios, and the model calculated scenarios with randomly sampled initial parameters for April, May, and June in 1998, 2000. The comparison result is shown in Fig. 10. The colored solid circles represent the real operation scenarios for each month. The hollow star symbols are the objective function values for 25 independent initial parameter sets. For better illustration, the symbols for the model simulated values for each month are all plotted with the same color (black). Nevertheless, the model simulated points of each month are closely clustered around the corresponding real operation points. The results indicate the model is able to give reasonably accurate simulations with randomly sampled parameters with the initial settings mentioned above.
To apply the proposed optimization scheme to the model, we set the population size to 128 with 64 individuals in each complex, and the maximum of function evaluation to 10,000 as one of the stopping criteria. The simulation results for the accumulated daily storage and net electricity generation in each month of OTC are shown in Fig. 11, in which each of the solution points represents one feasible Feather River release strategy during a given month. Different SeE curve-fitting methods (polynomial, piecewise linearization and successive parabolic interpolation) are compared. The comparison results for different months using MOSPD are shown in Fig. 12.
We also present the non-dominated solution set for May-2000 in Fig. 13(a)e(b). Among the non-dominated solutions, we choose two extreme examples. One of the solutions maximizes total daily storage volume and is termed the storage optimal operation. The other one maximizes the total net electricity generation and is titled the electricity optimal operation. Their corresponding daily storage volume is plotted in Fig. 13(c)e(d).

Discussion
MOSPD is a better operation-support tool than the original MOCOM. According to Fig. 11, both MOCOM and MOSPD are able to generate daily operation strategies in the feasible space with two differences: (1) the non-dominated solutions with MOSPD are located towards higher objective function values and (2) the generated solutions from MOSPD are more uniformly and diversely distributed along the Pareto front. These differences imply that MOSPD is capable of generating better non-dominated solutions in the OTC's problem. The simulation results from Fig. 11 also show that, when a monthly total release volume is regulated by DWR, an efficient daily-release adjustment can still benefit the system's output regarding the potential total daily storage and net electricity generation during the month.
The SeE curve-fitting methods comparison results (Fig. 12) show that the polynomial fitting is a better method for reducing the residuals. The results generated using the polynomial fitting are closer to the assumed "true" (successive parabolic interpolation) compared to that using piecewise linearization. Here, the successive parabolic interpolation is used as a reference because the true Pareto front of the realistic system is unknown for most of the case. By successive fitting of the SeE measurement points with the parabolic function, the global Pareto front is approximated by the "near real" solution front. The distance, which measures how close one solution front is to the "near real" solution front, represents the errors caused by different fitting methods. The closer of the two fronts indicates a better explanation of the "near real" situation and vice versa.
Even though the daily releases are similar among different strategies ( Fig. 13 (a)e(b)), the daily storage volumes could be dramatically different, as shown in Fig. 13 (c)e(d). Therefore, the changes or adjustments are the crucial factors to influence the entire system regarding the storage volume and other targets. This fact is in agreement with our consultation with DWR staff saying that the Feather River releases are very important operating objectives with regard to the management of OTC's facilities and SWP. The reason for the large change in storage volume is that the daily release contributes a daily carryover of storage, which forces the storage to either drain or fill the reservoirs. The carryover of storage accumulates so quickly that in several days the level of reservoir storage reaches its maximum, such as that shown in Fig. 13(c) on the 8th, 13th and 19th days. Similarly, the carryover storage forces the storage volume to reach the lower bounds of the reservoirs on the 21st and 26th day in the optimal-electricity operation (Fig. 13(d)). Larger static water-head difference arises, and more electrical power is generated, when the water level reaches the lower boundary on these two days.
Moreover, the optimized solutions are able to provide better daily operation alternatives in response to dry and wet water supply conditions. If drought conditions are foreseen in the near future, the storage-dominated solutions can increase the capability of the system to deal with emergent water needs or potential water loan. The reservoirs in the system are holding a higher storage volume by reducing the release amount for certain days but supplementing on other days. The higher short-term storage volume is Fig. 10. The comparison of the objective function values between the real operation scenarios and 25 independent runs of the model generated scenarios using the randomly sampled initial parameters for April, May June in 1998, 2000, and 2001. able to execute emergent drought-response operations, such as water supply loaning and storage consolidation. Both these special operations require high-storage volume of one region so that the short-term (days or weeks) water transfer does not drain the water lenders. If the water-supply condition is above normal and no shortages are expected, electricity generation-dominated solutions result in increased hydropower generation in order for the SWP to pump and deliver the water to its users. The increased clean-energy production from hydropower sectors also helps to mitigate the greenhouse gas emissions because the power supply from coalfired and other forms of energy is replaced by the extra hydropower generation. For the average water supply condition, whether one objective surpasses the other one depends on decision maker's preference and consideration, which are difficult to be generalized at this point. Nevertheless, the compromised solutions are recommended because both strategies for drought mitigation and increased hydropower generation become important.

Conclusions
In this study, we aim to enhance the capability and strengthen the application of a Multi-Objective Complex Evolution (MOCOM-UA) global-optimization method on the OTC's reservoir system. We built the optimization model based on a realistic reservoir system configuration, engineering constraints and decision makers' goals, but with several simplifications and assumptions. The impact of the simplifications and assumptions on reservoir topography are analyzed by comparing different curve-fitting methods. Different from the traditional procedure that separately considers algorithm development and modeling of the real-world problem, we try to narrow the gap between theoretical development and real world implementation by improving the original MOCOM algorithm, and build a two objective reservoir model for the OTC problem. This study provides an integrated platform to exhibit choices in a more transparent and clear format to decision makers in OTC. In detail, OTC's reservoir operation problem is studied, in which reservoir topography and hydro power generations are explicitly formulated, which dramatically reduces the errors when estimating the SeE relationship. Different curve-fitting methods are compared. The goal is to make the problem formulation as realistic as possible so that decision makers gain more confidence about the simplifications and assumptions in modeling the real world problem. On the algorithm development side, an improved algorithm, titled Multiobjective Shuffled Complex Evolution with Principal Component Analysis and Crowding Distance Operator (MOSPD), is developed to meet decision makers' requirement of accuracy and diversity on the obtained solutions. Test results give the following conclusions: (1) Test results show that the newly developed MOSPD algorithm is able to generate better non-dominated solutions, based on the diversity and the convergence performance criteria, as compared to the MOCOM algorithm. Both the Fig. 11. Simulation results for MOCOM and MOSPD for April, May, and June in 1998, 2000 diversity and convergence criteria are important in the decision making process allowing water managers to choose the most appropriate reservoir-operation options. We credit the resulting improvements to the effectiveness of the newly introduced modules, namely, the "possibility-adjustment" and the "dimension monitoring and restoring" modules. (2) The comparison among the optimal order of polynomial fitting, linearization and successive parabola fitting helps us to understand the impact of simplifications and assumptions in the way the real reservoir topography is mathematically represented. The results again confirm the claim by Labadie (2004) that non-linear challenges in optimal reservoirsystem management should be addressed directly by nonlinear programming, as well as the conclusion by Bay on et al. (2009) that linear simplification of the storageeelevation (SeE) curve can produce serious errors. (3) The optimal solutions derived by the proposed algorithm (MOSPD) are able to provide operation alternatives in response to different water supply conditions, as well as various preferences from decision makers. For the case study provided in this paper, the following overarching recommendations emerge. i. During dry conditions, the storage maximizing solutions are recommended in order to better respond to any special operating scheme triggered by drought.
ii. During wet conditions, the electricity maximizing solutions are recommended in order to mitigate power shortages and allow production of more clean energy. iii. The compromised solutions (in the middle ranges of the Pareto front) might be preferred by decision makers, based on their consensus preferences.
Finally, it is the authors' belief that the proposed approach, which combines the capabilities of advanced multi-objective optimization algorithms with more realistic (i.e., considering the nonlinearity and complexity) formulations of the system, can provide decision makers with the better picture of the range of options to choose from.

Limitations and future works
Regretfully, there are still several other non-linear aspects in modeling the OTC's problem (i.e., water rights (DWR, 2013b), environmental requirements (DWR, 2006), and engineeringoptimal design including heterogeneous hydropower units (Li et al., 2013), and non-stable short-term turbine efficiency influence (Diniz et al., 2007)), which are not fully considered in this study. These issues currently are either simplified or omitted for further study. In addition, more interactions between decision makers and algorithm developers are needed in order to allow for better and more realistic formulation of the real-world problem and greater appreciation by the algorithm developers about the complexity of issues facing decision makers. Other potential future work involves the adaptive changes of the constraints in the optimization process to obtain better Pareto optima (Piscopo et al., 2015).