Determination of cost coefficients of a priority-based water allocation linear programming model – a network flow approach

This paper presents a method to establish the objective function of a network flow programming model for simulating river–reservoir system operations and associated water allocation, with an emphasis on situations when the links other than demand or storage have to be assigned with nonzero cost coefficients. The method preserves the priorities defined by rule curves of reservoir, operational preferences for conveying water, allocation of storage among multiple reservoirs, and transbasin water diversions. Path enumeration analysis transforms these water allocation rules into linear constraints that can be solved to determine link cost coefficients. An approach to prune the original system into a reduced network is proposed to establish the precise constraints of nonzero cost coefficients, which can then be efficiently solved. The cost coefficients for the water allocation in the Feitsui and Shihmen reservoirs’ joint operating system of northern Taiwan was adequately assigned by the proposed method. This case study demonstrates how practitioners can correctly utilize network-flow-based models to allocate water supply throughout complex systems that are subject to strict operating rules.


Introduction
The allocation of water in river-reservoir systems usually involves a number of priority-based decisions, which include water rights, reservoir operating rules, commitments and negotiation between stakeholders, preferences for the conveyance of water and other requirements. Such systems usually comprise reservoirs, weirs, river channels, canals, diversion tunnels, pipelines and treatment plants as well as the demands of different purposes. The configuration of a regional system may extend to include multiple reservoirs, transbasin diversion and in-stream flow requirements at different reaches. Such modeling is further complicated by the need to determine the ideal means of regulating flow, such that demands are satisfied according to assigned priorities, while minimizing the residual water flowing into the receiving water body to ensure the efficient utilization of water resources. The means by which water is moved must also conform to the associated conveyance capacity.
Solving the above problem requires a clear identification and proper modeling of the allocating rules that account for every possible combination of supply and demand conditions (Ilich, 2008). A common approach is to utilize optimization methods (Yeh, 1985;Labadie, 2004;Rani and Moreira, 2010), among which the most widely applied is the linear programming (LP). This approach relies on LP to find the optimal feasible way of routing water in a regional system, given that the allocation objective, governing equations of physical water movement and operational constraints are appropriately linearly formulated. This formulating process requires sufficient knowledge of the optimization method as well as the under-analysis problem to transform the physical and operational features into mathematical representation. Moreover, satisfying the allocating rules usually requires a trial-and-error process to determine the most appropriate set of weighting factors or cost coefficients, which multiplying with respective allocated water constitutes the objective function. The lack of a systematic and precise way to establish and interpret the objective function may prevent the model from being entrusted or accepted by all involved stakeholders. For example, Juízo and Lidén (2010) reported the experiences of implementing an optimization-based model on transboundary water allocation in South Africa. They found that "the results from the system analysis tool are not easily understood by the stakeholders, and government representatives of different countries bear some suspicion about the results". In order to resolve this problem, two other nonoptimization-based models were evaluated and compared with the original one. Nevertheless, the authors still could not conclude on which model was more adequate for their case due to the structurally differences of simulating water allocation priorities in different models.
As a specialization of LP, network flow programming (NFP) only focuses on solving a specific subset of general LP problems that can be formulated in a more restrictive format. This loss of generality allows the resources allocation problem to be visually and precisely displayed by the network structure, and gains in return higher computational efficiency and easier comprehension of the priority-based allocation mechanism. These characteristics have prompted model developers to incorporate NFP into many general models (Evenson and Moseley, 1970;Sigvaldason, 1976;Labadie et al., 1986;Martin, 1987;Kuczera and Diment, 1988;Brendecke, 1989;Chung et al., 1989;Andrews et al., 1992;Wurbs, 1993;Andreu et al., 1996;Yerrameddy and Wurbs, 1996;Fredericks et al., 1998;Ilich et al., 2000;Dai and Labadie, 2001;Chou and Wu, 2010). The NFP represents the physical aspect of a water resources system as a directed network G(N, L), where N is the set of n nodes and L is the set of m links. The formulation of a minimum cost NFP problem can be expressed as (Ahuja et al., 1993) where i and j are the indices of node; (i, j ) is the link from the tail node i to the head node j ; x ij represents the amount of flow on link (i, j ); c ij is the unit shipping cost along link (i, j ); l ij and u ij are the lower and upper limits on flow in link (i, j ).
In a NFP-based water allocation model, nodes can represent storage or nonstorage points of confluence or divergence, and links represent reservoir outlet works, channels or pipes, water consumption, and carryover storage. Equation (2) indicates the continuity and availability of water at a node, for it states that the flow out of the node should equal to all incoming water. The upper and lower limits of a link represent its physical flow capacity, thus Eq. (3) states the transportability of water conveyance. The cost coefficient promotes flow routes that minimize net cost, thus determining the most preferable allocation of water supply with respect to a given allocating rule. Thus, correct assignment of link cost coefficients to reflect respective priorities is a necessary condition for any effective applications of not only NFP but LP-based water allocation models. Most common applications directly assign the cost coefficients related to the links of carryover storage or water consumption to represent the priorities of associated stakeholders. However, there are situations when internal links other than demand or storage have to be assigned with nonzero costs in order to achieve specific allocation requirements, such as water conveyance preference or surplus water diversion. This type of assignment is not straightforward for practitioners with little theoretical background, especially when forced to deal with a regional system of multiple reservoirs, water conveyance routes, instream flow requirements and transbasin water diversions.
The concept of developing a method for establishing cost coefficients of NFP models to adequately represent water allocation priorities was originally proposed by Israel and Lund (1999). Ferreira (2007) further broadened the scope for more general LP problems by demonstrating how different types of side constraints and variables in the LP formulation may affect the priorities defined by the cost coefficients of links in the NFP subset. These previous works represented the priority requirement as a set of rules. The rules were compiled into an LP problem that is solved as a means of initializing the actual allocation model (Ferreira, 2007). The present study follows and expands upon this principle with the proposal of additional allocation rules and a path-enumeration algorithm to facilitate automation of the cost-determination procedure. The presented rules allow one to simulate such water allocation priorities as reservoir rule curves, storage allocation among multiple reservoirs, preferred water mains, and transbasin diversion of surplus water. Path enumeration analysis is adopted to convert user-specified water supply allocation rules into a set of constraints; solving these constraints yields the cost coefficients that adhere to all specified rules. Further, an approach to prune the original system into a reduced network is proposed to establish the precise constraints of nonzero cost coefficients, which can then be efficiently solved. This pruned procedure thus functions successfully to efficiently initialize an effective application of water allocation models.

Alternative approaches: linear programming vs. network flow programming
The following presentation of methodology uses a NFP framework to demonstrate the procedure of determining cost coefficients. This concept is helpful to interpret the establishment of an objective function for more generalized LP-based models. One of the major differences between these alternative optimization approaches in modeling water resources allocations is how the non-NFP constraints, which cannot be represented by Eqs. (2) and (3), are incorporated. These constraints usually originate from the need to simulate physical water movement processes, such as return flows, flow losses, reservoir evaporation, and channel routing effects. In pure NFP-based models, these features have been handled through the use of successive iterations (Ilich, 2008). These iterative processes are external to the algorithmic solving procedure. Usually the lower or upper limits of links are iteratively adjusted to meet non-NFP constraints; thus, the priorities specified by link costs are unchanged during iterations. By contrast, an LP solver can directly incorporate non-NFP features into the formulation and the algorithmic solving procedure. However, this flexibility may impair the characteristic of priority-based water allocation of NFP. One simple example is that water may be allocated to a junior-priority demand with less flow loss, rather than a senior demand with greater flow loss, if the objective function is not appropriately set up in the LP formulation. Another example is the effect of channel flow routing, which may be easily modeled by the Muskingum method and incorporated into an LP formulation. Suppose that there are two demands located at the upstream and downstream ends of a river channel, respectively, with junior and senior priorities. The travel time required for water to flow through the channel from the location of upstream (junior) demand to downstream (senior) demand exceeds the unit time period of an LP-based simulation model. The portion of water that does not reach the point of downstream demand cannot explicitly contribute to the objective function in the current unit time period. The solution to this issue, similar to that for the flow loss case, consists of allocating water to the junior demand first instead of maximizing satisfaction of the senior downstream demand, if the discrepancy between their assigned cost coefficients is not large enough to compensate for the retained and ineffective portion of water. While NFP-based models are still widely utilized, several general software packages have updated their optimization engines with LP solvers to manage the rising demand for simulating non-NFP constraints and variables. Some examples include CALSIM (Draper et al. 2004), OASIS (Hydrologics, Inc., 2009) and WEAP (Stockholm Environment Institute, 2011. Nonetheless, the impacts of non-NFP features on water allocation have not been adequately discussed; only Israel and Lund (1999), Labadie andBaldo (2001), andFerreira (2007) have addressed this topic. Since non-NFP constraints must be strictly satisfied, they could be regarded as a higher level of priorities that would supersede and may disturb the priorities originally defined in the NFP subset as stated by Ferreira (2007). A desired resolution may be to achieve a simultaneous satisfaction of these two levels of priorities, if such a condition is feasible. In order to achieve this goal, the impact of non-NFP features on the allocation mechanism must be explicitly incorporated into the costdetermining procedure. Such as the two non-NFP constraints mentioned above, water transmission loss and flow routing can be modeled as the portion of water that is lost or delayed while allocating water to senior demands. This por-tion of water is ineffective to the objective function; the assigned link costs should be able to withstand these impacts to preserve the priorities of water allocation. For practical purposes, however, the present study focuses solely on determination of link costs for NFP-based modeling. Future research may extend to derive a comprehensive approach for more generalized LP-based models, thus accounting for all types of non-NFP constraints that may be encountered in real-world applications.

Framework of a network flow programming-based allocation model
NFP-based water allocation models can be used to allocate water over single or multiple time steps. For models that allocate water across multiple time steps, links connect reservoir nodes in different time periods to represent carryover storage. These models have been applied in reservoir sizing (Kuczera, 1989;Khaliquzzaman and Chander, 1997), capacity expansion (Martin, 1987;Gondolfi et al., 1997), the derivation of reservoir operating rules (Lund and Ferreira, 1996;Bessler et al., 2003), water transfer during droughts (Cheng et al., 2009), and the optimal real-time flood control operation of reservoirs (Braga and Barbosa, 2001). Single time step models allocate water only within an operational unit period, but the allocation is sequentially solved in every step during the simulation time horizon. Routing results produced in this manner are useful for quantifying the expected water supply situation and the risks of water shortage under the simulated conditions. This study discusses the assignment of cost coefficients for the single time step model. Figure 1 illustrates a water resources system as a network during a unit operational period. Virtual links illustrated by dotted lines satisfy Eq. (2), which specifies continuity equations of nodes, by conveying water into and out of the system. These virtual links signify the inflow of system, initial and carryover storage of reservoir, water consumed by the stakeholders, and the water body that receives surplus flow.

Principle in assigning cost coefficients and the necessity of preprocessing analysis
The cost coefficients of links, generalized in Fig. 1, quantify the relative priority of each respective water user. These cost coefficients must reflect the flow priorities associated with demand or storage under predefined operating conditions. One straightforward way to achieve this is to assign decreasing unit costs for demand/storage links of higher priority to ensure that the highest priority stakeholder is satisfied first in the cost minimization problem (Israel and Lund, 1999). The costs of internal links other than demand or storage can be kept as zero, thus the allocation will be solely driven by the relative value of costs on the virtual links as shown in Fig. 1. Nevertheless, there are situations when only assigning cost coefficients on demand or storage links is not enough to achieve the allocation requirements. One simple example is that minor costs such as −1 or +1 are commonly assigned on links where flow is to be encouraged, such as hydropower plant, or discouraged, such as routes with high transmission loss.
Another example is the transbasin diversion of surplus water, which requires diverting the required surplus water of a system into the adjacent system to enhance the efficiency of water utilization. An intuitive way to achieve this requirement is to use the iterative approach suggested by Labadie and Baldo (2001). This approach recommends a conceptual "flow-through" demand to be placed in the transbasin tunnel. This demand is given a lower priority than all demands or storage in the system to be diverted, which guarantees that transbasin diversions only occur once all demands in the original system are satisfied. According to the water supplied to the flow-through demand, iterations are then performed to artificially inject this diverted water into the adjacent system. Thus, transbasin diversion will work as long as the original system has surplus water, regardless of the hydrological condition of the other system. However, there is no need to perform diversion when both systems are in abundance of water, for the diverted flow will become surplus to the other system. Although the "flow-through" approach is capable of simulating physical water movement process such as nonconsumptive water usage, it may not properly model the operational features, such as adequate timing of diversion in this situation. This is especially critical when the transbasin diversion is charged with money; thus, unnecessary diversions should be avoided. Inevitably, satisfying the condition of surplus wa-14 comprehensive analyzing framework as shown in Fig. 2. Water allocation rules and cos 255 determining procedure is described in detail in the following section.

Rule 1: Trans-basin diversion of surplus water 262
Generally, the development of a new trans-basin water diversion project must n 263 impact existing users of the system. Fig. 3 depicts a simple example, in which on 264 surplus water in the system associated with reservoir B can be diverted for storage 265 reservoir A. Thus, the first rule allows users to specify a link in the network representin 266 ter diversion requires assigning a positive cost on the link of the transbasin tunnel, without using the flow through demand approach.
The determination of cost becomes more complicated if a combination of various allocation rules is involved, such as different operating rule curves for individual reservoirs, preferences of water conveyances in multiple locations, the allocation of multireservoir storage, and transbasin water diversions. When multiple links in the system have to be assigned with nonzero cost coefficients, the accumulation of costs along a flow path to a demand/reservoir might impair its priority, which is originally dictated by the cost of the virtual link. The connectivity between links of nonzero costs has to be identified to ensure that the sum of cost coefficients in paths to a water usage of higher priority is always less than the total costs of any path to a lower priority stakeholder. If the user cannot ensure assigning nonzero costs on the links to achieve the allocation requirements, a general preprocessing analysis will have to assume that the cost coefficient of every link in the system is unknown.
This study develops a procedure to establish the objective function of NFP-based water allocation models, in which all representative allocation rules encountered are considered. The allocation associated with reservoir operating rule curves and multireservoir storage balancing was preliminarily addressed in Chou and Wu (2011). These two rules are more elaborated in this paper, with two additional rules, transbasin surplus diversion and water conveyance preference, being proposed to constitute the comprehensive analyzing framework as shown in Fig. 2. The water allocation rules and costdetermining procedure are described in detail in the following section.

Rule 1: transbasin diversion of surplus water
Generally, the development of a new transbasin water diversion project must not impact existing users of the system. Figure 3 depicts a simple example in which only surplus water in the system associated with reservoir B can be diverted for storage in reservoir A. Thus, the first rule allows users to specify a link in the network representing a way of distributing water with last priority. The priorities of all paths through this specific link are junior to any other paths to demands and storage in the system.
Let L be the set of all links, L D be the set of virtual demand links, L S be the set of virtual storage links in the network, and (L D + L S ) be the union of L D and L S . Define a path as a sequence of links without the repetition of head nodes, i.e., with no cycle in the path. Use R LP to represent the set of paths containing the specific link for the diversion of surplus water, and R L D +L S to represent the set of paths with the final links belonging to (L D + L S ). The mathematical formulation of priority requirement for surplus water diversion can be expressed as where (R L D +L S −R LP ) is the same as R L D +L S but excluding R LP , cost is a function used to calculate the sum of the cost coefficients of the links in a path, and cost (R LP ) represents the set of total costs for all paths in R LP . Equation (4) states that the largest cost conducted by paths that do not pass from the transbasin link is less than the lowest cost by passing from the transbasin link. Because the lowest priority should correspond to the largest cost under the framework of NFP, a set of cost coefficients that satisfies this condition should guarantee that the transbasin link will work only in case of surplus. For a total of np 1a paths in R LP where the kth path is represented as P 1a k , a Kronecker delta function can be used to represent if P 1a k contains link (i, j ): Suppose that (R L D +L S − R LP ) contains np 1b links and P 1b k represents the kth path in (R L D +L S − R LP ). Another Kronecker delta δ1b (i,j ) k can be used to represent if P 1b k contains link (i, j ): Equation (4) can then be expressed by the following constraints: where c (i,j ) is the cost coefficient per unit flow of link (i, j ), C Min1 represents the lower bound of the total costs of paths in R LP , C Max1 represents the upper bound of the total costs of paths in (R L D +L S − R LP ), and ε is an arbitrary positive integer specified by the user.

Rule 2: priorities between water usages and reservoir storage
The basic framework of water allocation in the water resources system is the priorities between water usages and reservoir storage. The priorities may be defined by water rights, judicial or legislative actions to protect specific water usages, private agreements between stakeholders or the operating rule curves of reservoirs. Chou and Wu (2011) illustrated the setting of priorities between demands and storage for the operating rule curves commonly adopted in individual reservoir operating systems of Taiwan. The proposed mathematical formulation was as follows.
Assume that (L D + L S ) is the set that consists of all virtual demand and storage links. (L D + L S ) (k) is the link prioritized kth among (L D + L S ). Equation (10) prioritizes all virtual demand and storage links that comprise a water supply network as follows: In Eq. (10), the set R L D +L S (k) consists of all potential flow routes with final link as L D + L S (k), R LP is the same as defined in Eq. (4) of Sect. 3.1; and m d +m s represents the number of links in (L D + L S ). Equation (10) states that the largest cost among paths to a senior priority demand or storage is less than the lowest cost conducted by paths to a junior priority water usage. It thus guarantees finding coefficients that project the defined priorities. The following constraints can be established from the concept of Eq. (10), derived by a similar process of converting Eq. (4) into Eqs. (7)-(9) as shown in Sect. 3.1. (i,j ) k,l c (i,j ) ≤ C Max2 k l = 1, .., np 2,k , where C Max2 k and C Min2 k define the feasible range of net conveyance costs for flow paths in R L D +L S (k) − R LP ; the Kronecker delta function δ2 (i,j ) k,l indicates whether the lth flow path of R L D +L S (k) − R LP includes the link (i, j ), np 2,k is the number of paths that exist in R L D +L S (k) − R LP , and ε is the same as in Eq. (9), which is used to maintain an interval of costs between consecutive priorities.

Rule 3: preferences in water conveyance
Although there are multiple ways to meet a demand, for water the routes with less transmission loss, lower operating costs, and the potential for additional hydropower generation are generally preferred. This rule allows users to specify the priorities of water conveyance through paths between two specific nodes. For example, possible paths between the reservoir and demand nodes in Fig. 4 are listed in the sequence of their priorities as follows: (1) Suppose that there are np 3 possible paths between the specified source and target nodes. We assume that these paths are arranged in sequence according to their conveyance priorities, i.e., if P 3 k represents the kth path, then water conveyance through P 3 k should be prior to P 3 k+1 . The function δ3 (i,j ) k indicates whether P 3 k includes the link (i, j ). The following constraints can then be established: where C Max3 k and C Min3 k represent the upper and lower bounds of costs associated with the paths between the specified source and target nodes.

Rule 4: priorities in multireservoir storage allocation
The operation of a multireservoir system involves allocating water from multiple reservoirs to satisfy the joint demand. The respective priority rankings for carryover storage of each reservoir determine which reservoir should be used first to satisfy demand throughout a multireservoir system. The operation of a multi-reservoir system involves allocating 354 multiple reservoirs to satisfy the joint demand. The respective priority r 355 carryover storage of each reservoir determine which reservoir should be u 356 satisfy demand throughout a multi-reservoir system. For example, Fig.5 depi  357 with two parallel reservoirs, Reservoirs A and B, which both can provide w 358 joint demand. Operating rules of this two-reservoir system dictate that joint 359 supplied by allocating water from available sources in the following order: (1 360 system dictate that joint demand be supplied by allocating water from available sources in the following order: (1) first from Weir C until it has been emptied; (2) then from Reservoir A, provided that its water level is over its lower limit of rule curve; and (3) finally, from Reservoir B. Accordingly, the storage components can be listed in the sequence of their associated priorities as (1) the storage under the lower limit of Reservoir A, (2) the storage of Reservoir B, (3) the storage over the lower limit of Reservoir A and (4) the storage of Weir C.
Assume that L S (k) represents the kth-priority link in the set of storage links, L S . The priority constraint for allocating storage in a multireservoir system can be expressed as follows:

Fig. 5 Example of a multi-reservoir system
Assume that L S (k) represents the k th -priority link in the set of storage links, L S .
The priority constraint for allocating storage in a multi-reservoir system can be expressed as follows:  15) is explained as follows: suppose that there is one unit of water initially stored in the reservoir for each of the storage links. The water can either be released to satisfy the joint demand or retained in the reservoir to contribute to the associated carryover storage. The left-hand side of Eq. (15) represents the largest cost induced by storing water in the senior storage link (index k) and releasing water from the junior storage (index k+ 1) to supply joint demand. However, the right-hand side represents the lowest cost induced by storing and releasing water in the converse way. The inequality ensures that a junior storage will release water in a higher priority to supply joint demand.
According to similar process as shown from Eq. (4) to Eqs. (7)-(9), the following constraints can be established: (i,j ) k,l c (i,j ) ≤ C Max4a k l = 1, ..., np 4a,k ; k = 1, ..., m s , where C Max4a k and C Min4a k define the feasible range of net conveyance costs for flow paths represented by (R L S (k)→JD − R LP ); C Max4b k and C Min4b k define the feasible range of net conveyance costs for flow paths represented by (R L S (k) − R LP ); the functions δ4a (i,j ) k,l and δ4b (i,j ) k,l indicate whether the lth flow path of (R L S (k)→JD − R LP ) and (R L S (k) − R LP ) include the link (i, j ) respectively; np 4a,k and np 4b,k are the numbers of paths in (R L S (k)→JD − R LP ) and (R L S (k) − R LP ) respectively; and ε is the same as in Eqs. (9) and (11).

Rule 5 (default): minimization of surplus water
The proposed method penalizes any water into the final receiving body by the following requirements: where L T is a set that includes all terminal links originated from the node representing the water receiving body; R L T is a set that consists of all possible flow paths, each of which has a final link that belongs to L T . Equation (19) states that the lowest cost by paths that include the virtual terminal link is greater than zero, and Eq. (20) states the largest cost to a virtual demand or storage link is less than zero. In this manner, the NFP algorithm will then try to allocate unregulated flows to water users, and release spill flows from reservoir only if absolutely necessary to prevent inducing a positive cost. The following inequalities can then be established: where, δ5 (ij ) k and δ6 (i,j ) k are Kronecker delta functions to represent whether link (i, j ) is in the kth path in R L D +L S and R L T , respectively; np 5 is the number of paths in R L D +L S ; and np 6 denotes the number of paths in R L T .
Furthermore, we assume that the cost coefficients of all links other than demand, storage and terminal are greater than 0:

Linear programming for determining cost coefficients
The constraints, Eqs. (7)-(9), (11)-(12), (13)-(14), (16)- (18) and (21)- (23), define the feasible region for cost coefficients. Linear programming (LP) can be employed to solve the problem, by coupling the constraints with the following objective function: Equation (24) will keep the costs of links other than storages, demands and terminals to be zero as long as feasible. Only a few links will be assigned with nonzero costs when absolutely necessary. For example, rule 3 may require assignment of nonzero costs on particular links to discourage flow through routes with high loss rates. The assigned cost will then be minimized to be +1 based on the objective function and Eqs. (13) and (14). Under this setting, the allocation of water will be primarily dictated by the costs of virtual links, while the minor costs on particular nonvirtual links guide local flow conveyance.

Determination of values of the Kronecker delta functions
The Kronecker delta functions for each link as described in Sects. 3.1-3.5, can be established using the path enumeration algorithm of Kroft (1967). Here a path refers to a sequence of nodes such that from each node there is a link to the next node in the sequence. Furthermore, there should be no cycle, i.e., repetition of nodes, in the path. Repeated identification of possible paths between different associate nodes can help determining the values of the above Kronecker delta functions. The computing procedure of Kroft's algorithm is provided in Appendix A.

Case study
The proposed method was applied to determine cost coefficients of the NFP model for simulating the joint water allocation of the Hsintein and Tahan rivers' water resources system of northern Taiwan. This case study simulates projected conditions of the given system in 2021. The Feitsui Reservoir, with an effective storage capacity of 336 × 10 6 m 3 , is located on Peishih Creek, one of the two major upstream tributaries of the Hsintein River. It serves mainly to supply the demand for domestic water in the Taipei (TP) district. Downstream from the confluence of Peishih and Nanshih Creeks are the Cihukeng, Chihtan, and Chintan weirs, which serve to regulate upstream flow and raise the water level for the diversion of water into three treatment plants. The Cihukeng weir also serves to raise the water level to divert flow into the off-channel Cihukeng hydropower plant through a manmade canal. The tail-water from the hydropower plant is then diverted to the downstream Chintan weir. The other river in the joint operating system, the Tahan River, has its own reservoir, the Shihmen Reservoir. The capacity of the Shihmen Reservoir is 215 × 10 6 m 3 according to the survey in 2011. It was designed for irrigation, hydropower generation, public water supply, and flood moderation. Downstream from the Shihmen Reservoir are its afterbay and the Yuanshan weirs, which serve to regulate the reservoir release. The Shanshia pumping station on the Shanshia River, which is a tributary of the Tahan River, can also support public water supply in this region.
The primary demands for water in the Shihmen Reservoir system are irrigational and the public demand of southern, northern Taoyuan (TY) and Pan-Hsin (PH) districts. The Pingcheng, Longtang, and Shihmen treatment plants withdraw raw water from the Shihmen Reservoir and supply the southern TY district. The northern TY district is supplied by the Danan treatment plant, which withdraws raw water from the Yuanshan weir.
The Tahan River and Hsintien River systems jointly supply the public demand from the PH district. The Panhsin treatment plant receives raw water from both the Yuanshan 27 called Limogan Weir, and a trans-basin tunnel upstream of Nanshih Creek. It aims to 477 divert surplus water from Nanshih Creek to an upper section of Sanshia River, thereby 478 increasing the water utilization efficiency through joint operations. The network of this 479 water resources system is depicted in Fig. 6. weir and Shanshia pumping station. The Hsintien River system will provide a maximum of 1.01 million m 3 day −1 of treated water to the PH district after 2016 through the underconstruction transbasin pipeline of the "Pan-Hsin Water Supply Improvement Plan, Phase II" (PH-Phase II).
There is also a transbasin raw water diversion project being planned for Nanshih Creek in the upstream of the Hsintein River, which will focus on building a diversion weir, called Limogan weir, and a transbasin tunnel upstream of Nanshih Creek. It aims to divert surplus water from Nanshih Creek to an upper section of the Sanshia River, thereby increasing the water utilization efficiency through joint operations. The network of this water resources system is depicted in Fig. 6.

Priority requirement for transbasin water diversion
The diversion link of the Limogan weir is specified as the last priority link of rule 1, because it should only divert surplus water from Nanshih Creek. This setting ensures that the transbasin tunnel will not withdraw water originally intended to meet the demands of the Hsintein River system.

Priority requirement for reservoir operating rule curves
The rule curves of the Feitsui Reservoir include the severe limit (SL), lower limit (LL), middle limit (ML) and upper limit (UL). The Feitsui Reservoir administration specifies the following conditions for operation in 2021 (Chou and Wu, 2011): 1. While reservoir water level is below the SL, it only has to provide 80 % of TP demand.
2. While reservoir level is above the SL but below the LL, it only has to provide 80 % of TP and PH demands.  Shihmen Reservoir operating rule curves must comply with the following criteria: 513 1. While reservoir level is below the SL, it only has to provide 50% of irrigational and 514 80% of TY and PH demands. 515 2. While reservoir level is above the SL but below the LL, it only has to provide 75% of 516 irrigational and 90% of TY demands. 517 3. While the reservoir level is above the LL, 100 % of the TP and PH demands should be satisfied.

4.
While the reservoir level is raised to range between the ML and UL, extra water may be released for peak-hours hydropower generation.
5. Sufficient water should be released to support fullcapacity hydropower generation while reservoir level exceeds the UL.  -190 (= -270+80), which is larger than the cost of simply storing that water in 560 the storage facilities in the Tahan River system. 561 Assume that both the Feitsui and Shihmen Reservoirs each have one unit of water 565 and that the Feitsui water level is higher than its SL. If the water from Shihmen 566 Reservoir is allocated to supply 80% of the joint demand, the other one unit of water can 567 be stored in Feistui Reservoir to achieve the minimum unit cost of -280. On the other 568 hand, the unit cost of supplying joint demand with Feitsui Reservoir water (and thus 569 retaining Shihmen Reservoir storage) is only -290. Hence, minimum-cost NFP-based 570 water allocation ensures that the joint demand will be satisfied by the Feitsui storage in a 571 higher priority, provided that its water level exceeds the SL. 572 Figure 8. Assigned coefficients based on conditions specified in  unable to attain the SL. Under these conditions, the alternate supply source, the Shihmen Reservoir, will supply 80 % of the PH district demand. The cost of supplying the remaining the PH demand would be −190 (= −270 + 80), which is larger than the cost of simply storing that water in the storage facilities in the Tahan River system.
Assume that both the Feitsui and Shihmen reservoirs each have one unit of water and that the Feitsui water level is higher than its SL. If the water from the Shihmen Reservoir is allocated to supply 80 % of the joint demand, the other one unit of water can be stored in the Feistui Reservoir to achieve the minimum unit cost of −280. However, the unit cost of supplying joint demand with the Feitsui Reservoir water (and thus retaining the Shihmen Reservoir storage) is only −290. Hence, minimum-cost NFP-based water allocation ensures that the joint demand will be satisfied by the Feitsui storage in a higher priority, provided that its water level exceeds the SL.
The transbasin diversion link in Fig. 8 has a positive cost coefficient of +180. The minimum total cost of paths through this link is −180, which is the sum of the costs of the diversion link and the highest priority demand in the Tahan River system. The lowest priority in the Hsintein River system is storage in weirs, each of which has a cost of −210. Thus the model will not allocate water from Nanshih Creek unless all of the weirs of the Hsihtein River are full. In other words, the transbasin tunnel will only divert surplus water from Nanshih Creek.
In the joint operation of Fig. 8, the Feitsui Reservoir is the primary regular source and the Shihmen Reservoir provides the backup source for the PH district. Another operating strategy is to maintain the storage of these two reservoirs at the same intervals as their individual rule curves. For instance, the storage zones between the LL and SL of both reservoirs would share the same priority. Based on this concept of storage balancing joint operation, the virtual storage links are listed in the order of their associated priorities as 34 591 Fig. 9 Cost coefficients for storage balancing of two reservoirs 592 593 Based on Fig. 9, possible joint operating scenarios include the following: 594 1. Any water over the UL in the Shihmen Reservoir will be allocated to the PH district to 595 meet 80% of its full demand, provided that Feitsui level does not exceed its UL. 596 2. When the level of Shihmen Reservoir is between its UL and LL, the Feitsui Reservoir 597 will satisfy the joint demand as long as the its level exceeds the ML. However, if 598 Feitsui storage is unable to attain the LL, then water from the Shihmen Reservoir will 599 be allocated to meet 80% of PH district demand in a higher priority. 600 3. Provided that the Shihmen Reservoir water level ranges between the SL and LL and 601 the water level in the Feitsui Reservoir exceeds the LL, water from Feitsui Reservoir 602 will be allocated to PH district demand. Shihmen Reservoir water will be released to 603 Under this setting, the reservoir with the higher storage is charged with supplying the joint demand to maintain the storage of the two reservoirs in the same interval. The analyzed cost coefficients based on the storage balancing joint operation are illustrated in Fig. 9.
Based on Fig. 9, possible joint operating scenarios include the following: 1. Any water over the UL in the Shihmen Reservoir will be allocated to the PH district to meet 80 % of its full demand, provided that Feitsui level does not exceed its UL.
2. When the level of the Shihmen Reservoir is between its UL and LL, the Feitsui Reservoir will satisfy the joint demand as long as the its level exceeds the ML. However, if Feitsui storage is unable to attain the LL, then water from the Shihmen Reservoir will be allocated to meet 80 % of the PH district demand in a higher priority.
3. Provided that the Shihmen Reservoir water level ranges between the SL and LL and the water level in the Feitsui Reservoir exceeds the LL, water from the Feitsui Reservoir will be allocated to the PH district demand. Shihmen Reservoir water will be released to independently satisfy 80 % of joint demand only when the Feitsui water level drops below its SL.
4. When the Shihmen Reservoir water level drops below the SL, the Feitsui Reservoir will independently fulfill the PH district demand provided that its own water level exceeds the SL. If the Feitsui Reservoir water level is below the SL, then the Shihmen Reservoir water will be allocated to ensure that 80 % of the PH demand is satisfied.
In addition to the allocation priorities defined by operating rule curves and joint operating rules, preference for flow through a hydropower plant can be simulated by directly assigning a negative unit cost to the links connecting to the runof-river or reservoir hydropower plants to encourage associated flows. Because the interval of costs between consecutive priorities of demands or storage is set as 10, this unit cost will not impair the priority requirements by the above rules, as long as the accumulations of minor costs to demands or reservoirs are within the range between −10 and 10.

A pruned analysis procedure
In the aforementioned analysis procedure, the bulk of the computational load is expended on network path enumeration analysis. For a complete network, in which every pair of distinct nodes is connected by a unique link (as an extreme example), if there are n nodes in the network, then the number of links will be 2 × C n 2 , resulting in n−2 i=1 C n−2 i paths between any two distinct nodes. This means that the number of paths would grow exponentially with an increase in the number of nodes for such a dense network. The enormous number of resulting paths would not only require considerable time for enumeration, but would also expand the size of the subsequent LP problem. Path enumeration is required because the cost coefficient of every link is assumed to be unknown in the default condition. If additional conditions could be included, such as the assignment of only a few links with nonzero costs and the costs of other links set at 0, then a simpler analysis procedure could be employed to reduce the required computational load.
Using G(N, L) to present the network under analysis, which is defined by a set N of n nodes and a set L of m links, suppose that there are m P nonvirtual links within L that are assigned with nonzero costs and m P < m. Define L P as the set containing these specified links, N PT and N PH as the sets of tail and head nodes of links in L P , respectively, and (N D +N S +N T ) as the set that contains all nodes which represent demands, reservoirs or final water receiving bodies in N , and (L D + L S + L T ) as the set of demand, storage or terminal links. Then the cost-determining procedure can be simplified as follows: 1. From each of the nodes that convey inflow into the system, use the depth first search (DFS) algorithm to identify the downstream reachable nodes in G(N, L − L P ).
The detail of the DFS algorithm can be found in Ahuja et al. (1993).

2.
A fictitious node, denoted as node f , is created. If node i ∈ (N D + N S + N T ) is identified to be reachable from inflow nodes in the previous step, then a fictitious link (f , i) is created. This fictitious link serves to replace all paths to node i which consist of only links with zero cost in G(N, L); define L F as the set that contains these fictitious links.
3. Use DFS to identify the downstream reachable nodes in G(N, L − L P ) from the head node of each link in L P .
4. Suppose that link (i, j ) belongs to L P and node k belongs to either N PT or (N D + N S + N T ). If k can be reached from j in G(N, L − L P ), then a fictitious link (j ,k) is created and added into L F . These fictitious links represent the connectivity between links with nonzero costs.
5. Establish a reduced network G (N , L ) in which N is the union of N PT , N PH , (N D + N S + N T ) and node f , and L is the union of L P , L F and (L D + L S + L T ).
6. The same procedure described in Sect. 3 can be followed to determine the cost coefficients of links in L P and (L D + L S + L T ), except that G(N, L) is replaced by G (N , L ).
The above procedure takes advantage of the fact that total costs of a path are determined only by the links with nonzero cost coefficients in the path. Thus the enumeration of paths containing all links in L can be reduced to only enumerating feasible combinations of links in L P and (L D +L S +L T ). Because DFS is a basic algorithm with worst-case complexity as only O(m), the reduced network G can be efficiently established from the original network G. The scale of G should be much less than G because typically m P m. Thus, enumerating paths in G will require much less computational time and the size of the consequent LP problem can be greatly reduced.
This pruned procedure was employed to finally evaluate the two illustrative problems in Sect. 4. In these final evaluations, only the transbasin diversion link and the links connecting to 20 % joint demand are specified with nonzero costs. The original system was pruned into a reduced network similar to the schematic shown in Fig. 8. For each problem, the number of constraints in the LP formulation was reduced from the original 3227 to only 486. The analysis' results using the pruned procedure were identical to those as illustrated in Figs. 8 and 9.

Conclusions
This paper presents a methodology for determining the cost coefficients of the objective function of a NFP-based model for simulating river-reservoir system operations and associated water allocation. This issue is of great importance because adequate simulation of water allocation rules is the key to successful implementations of any water allocation model. Among the many studies on water allocation within reservoir-river systems in the literature, this paper is one of the very few that explicitly studies how to appropriately set up the objective function for a NFP-based simulation model. The assignment of cost coefficients was usually performed intuitively, as practices of art, by researchers. This issue is F. N.-F. Chou and C.-W. Wu: Priority-based water allocation linear programming model treated in a scientific manner in this paper, with systematic presentations of representative allocation rules encountered in real-world applications. A general procedure is proposed to solve the problem. Although additional analysis efforts are required, the obtained coefficients guarantee that the allocation requirements are satisfied. Thus the possibly timeconsuming trial and error process to check the validity of assigned costs can be avoided.
For an experienced analyst, the adequate assignment of cost coefficients may be done without any preprocessing procedure. But this is not necessarily true for practitioners with less theoretical background, especially when they are dealing with systems of complex networks and allocation rules.
To a system consists of multiple reservoirs and a transbasin diverting tunnel or pipe as shown in the case study, achieving surplus water diversion and storage allocation inevitably requires assigning nonzero costs on internal links other than demands or storage. This practice is not as straightforward as for systems with simple allocation priorities on demands or reservoir storage. Even for an experienced practitioner, there is always a risk of the wrong assignment of costs due to the variety and complexity of water resources systems. The proposed procedure can also serve to validate the effectiveness of the intuitively assigned costs.
Furthermore, if the links to be assigned with nonzero costs can be specified in advance, a simpler procedure can be employed to reduce the computing effort of preprocessing analysis. This procedure prunes the original system into a reduced network. Thus, the time required to establish and solve the constraints of cost coefficients can be greatly shortened, which further increases the merit of the proposed method.