Optimal Booster Station Design and Operation under Uncertain Load

Given industrial applications, the costs for the operation and maintenance of a pump system typically far exceed its purchase price. For finding an optimal pump configuration which minimizes not only investment, but life-cycle costs, methods like Technical Operations Research which is based on Mixed-Integer Programming can be applied. However, during the planning phase, the designer is often faced with uncertain input data, e.g. future load demands can only be estimated. In this work, we deal with this uncertainty by developing a chance-constrained two-stage (CCTS) stochastic program. The design and operation of a booster station working under uncertain load demand are optimized to minimize total cost including purchase price, operation cost incurred by energy consumption and penalty cost resulting from water shortage. We find optimized system layouts using a sample average approximation (SAA) algorithm, and analyze the results for different risk levels of water shortage. By adjusting the risk level, the costs and performance range of the system can be balanced, and thus the system's resilience can be engineered.


Introduction and State of the Art
While for pumps used in automotive engineering or in private households, the production price and thus the investment costs are prominent, in industrial applications, life-cycle costs are in focus. Over the life span of a pump system, the costs for its operation, maintenance and for forced downtimes typically far exceed the purchase price. To find an optimal pump configuration with minimal life-cycle costs, a method optimizing the design and control of technical systems simultaneously is needed.
In this work, we adopt the method Technical Operations Research (TOR) [1,2], which is based on Mixed-Integer Programming (MIP). While MIP has been applied for the optimization of water distribution networks, c.f. [3] for an overview, uncertainty has rarely been investigated in this context.
In [4], a Quantified Linear Program formulation, where decision variables may not only be existentially but also universally quantified, c.f. [5], has been suggested to optimize pump systems under uncertain load. It allows to find the best-priced system with respect to a handful of frequent loads, out of all system configurations that are guaranteed to fulfill a huge number of infrequent loads, whatever the operational costs. Uncertainty due to pump failures is adressed in [6], where a booster station is optimized given uncertain load and multiple failure scenarios. The probability of the failure scenarios is derived based on the pumps' Mean Time Between Failures, and the system configuration is optimized for minimal investment, operation, and downtime cost.
With this contribution, we want to address the problem of designing resilient pump systems while preventing overdimensioning due to uncertainty. Therefore, we propose Stochastic Mixed-Integer Programming in combination with chance constraints. With this approach, we are able to design an energy-efficient and cost-optimal booster station under uncertain load demand. Moreover, it allows us to specify a risk level for water shortage. By adjusting this level, one may balance investment cost and resilience, and find the optimal system for one's specific use case.

Theoretical Background
The investigated optimization problem consists of a two-stage stochastic program combined with a chance constraint. We give a short overview of the underlying mathematical concepts but do not go into details. The interested reader is referred to [7]. Two-stage stochastic programming (TSSP) is a widespread optimization method that allows to consider uncertain parameters. The general mixed-integer nonlinear two-stage stochastic program is given by [8]: where ξ(ω) is a random vector with probability distribution P , and E P [Q(x, ξ(ω))] is the optimal expected value of the second-stage (recourse) problem: Here, all functions f 1 , are general linear or nonlinear functions and f 2 (·, ω), g 2 i (·, ω), h 2 i (·, ω) are measurable in ω for any fixed first argument. The variables x denote the first stage decisions which have to be determined prior to the realization of the uncertain parameters ω, whereas the variables y represent the recourse decisions that can be made after uncertainty is disclosed. Problem (1) seeks a first-stage decision x that minimizes first stage costs f 1 (x) and the expected cost E P [Q(x, ξ(ω))] of second-stage (recourse) decisions of Problem (2).
Chance-constrained programming (CCP). While stochastic programs allow to model uncertain parameters, their constraints are strict and thus required to be met in all scenarios. Since this may lead to solutions that are too conservative, chance constraint programming [9] has been investigated in engineering applications like water management [10] or unit commitment problems [11,12]. In CCP, solutions may violate so-called chance constraints to some extent. The general form is given by: where x ∈ X denotes the deterministic feasible decisions, and Eq. 4 is a chance constraint, which is fulfilled with confidence level 1 − ε, where ε ∈ [0, 1]. The feasible region of a chance-constrained problem depends on the user's choice of the risk level. For safety-critical applications, the risk level can be lowered which makes the feasible region smaller. This typically results in more costly solutions, which are on the other hand more resilient.
Sample Average Approximation. Sample average approximations have been used to obtain good candidate solutions for both stochastic programs (e.g. [13]) as well as chance-constrained programs (e.g. [14]). In both cases, the actual probability distribution of random variables (for stochastic programs the distribution of the uncertain data / for chance-constrained programs the distribution in the chance constraint) is replaced by an empirical distribution using computer simulation, e.g. Monte Carlo simulation. In case of the two-stage stochastic program described above, a Monte Carlo sample { ξ 1 , . . . , ξ N } is generated, and the expectation objective E P [Q(x, ξ)] of the original problem is replaced by the sample average N −1 ∑ N k=1 Q(x, ξ k ). In case of the chance-constrained problem, the chance constraint is replaced by an indicator function Once the sample has been generated, the corresponding SAA problems become deterministic optimization problems. This deterministic approximation is then solved with a suitable algorithm and serves as a candidate solution to the true problem. By repeating the sampling-optimization procedure several times, one can derive statistical confidence intervals on the obtained candidate solutions. Moreover, one can show that the number of samples needed to get a sufficiently accurate solution with high confidence is not too large. While traditionally SAA has been employed to solve either stochastic or chance-constrained problems, in [15], a first algorithm for solving models that contain both chance-constrained, as well as two-stage stochastic program features has been presented and applied to solve a unit commitment problem addressing uncertain wind power output. While this algorithm focuses on problems with linear objective and constraints, due to the technical and physical characteristics of the investigated system, we deal with nonlinear constraints and use the proposed algorithm as a primal heuristic.

Optimization Problem
In this paper, a nonlinear chance-constrained two-stage stochastic program (CCTS) is developed to optimize a booster station under uncertain load. The first stage of this program corresponds to finding optimal topology decisions with minimal investment cost. In the second stage of the program, goal is to find an optimal operation strategy for a given first-stage decision, while saving operation and penalty cost incurred by power consumption and water shortage, respectively. The chance constraint in the second stage ensures the feasibility of the first-stage decisions at a certain risk level.
As topology options, we consider all possible parallel connections that can be built, given a construction kit of different pump types, cf. Fig. 1. Specifically, the pump type selection and the number of pumps installed in parallel per pump type are considered as first stage decisions. The second stage decisions include operation decisions like the number of activated pumps, the rotating speed and the flow rate of each pump. The basis of our mathematical model follows [16]. In the pump system, water flows from the inlet to the outlet via pumps connected in parallel. Check valves are installed behind each pump to avoid reverse flow. The pressure at the outlet is considered higher than at the inlet. The system is modeled as a graph G = (V, E), representing the superstructure of all possible topology options. Its edges E represent technical components like pipes and pumps. Its vertices V represent connection ports. For each vertex, the flow conservation holds.
The water demand q at the sink is uncertain and modeled as a normally distributed random vector, i.e. q ∼ N (µ, σ 2 ), cf. Fig. 2.

104
Uncertainty in Mechanical Engineering III Mathematical Model. We formulate a nonlinear CCTS which combines a two-stage stochastic nonlinear mixed-integer problem with a chance constraint. Some assumptions made when building this mathematical model should be clarified: (i) Pumps are the only technical components considered. The number of available types, their operation characteristics and prices are known; (ii) The unit cost incurred by energy consumption and penalty for water shortage are fixed; (iii) We only consider parallel connections, and limit the number of pumps of the same type. (iv) For pumps of the same type, flow rate and rotating speed should be equal to ensure the highest efficiency [17]. (v) The water demand at the outlet of the booster station is the only parameter considered uncertain. We approximate the continuous probability distribution by a N -sample of discrete scenarios l = 1, . . . , N , but leave out scenarios l with a water demand q(l) < 0.
The list of decision variables, parameters and sets used in the CCTS is given by Tab. 1, the complete model is given by Eqs. 5 -25. We denote decision variables by capital letters, parameters by lower case letters. The only exception is the number of discrete scenarios N , which is also a parameter but represented by a capital letter.
The first-stage problem in Eqs. 5 -8 comprises the topology decision, i.e. which pump types should be chosen from the construction kit and how many pumps of each chosen type should be installed. Specifically, Eq. 6 makes sure that at least one pump type t is installed in the system. Eq. 7 constrains the maximum number of pumps of the same type. If pump type t is selected for the topology, Eq. 8 ensures that at least one pump of this type is bought.
The second-stage problem in Eqs. 9 -25 comprises the operation decisions, i.e. which pumps should operate in which demand scenario and if a demand should be fully fulfilled or not, given the penalty costs. The expected costs of the second stage, E [F (X , q)], are specified in Eq. 9, as the sum of operation cost incurred by energy consumption and the penalty cost incurred by water shortage over all discrete water demand scenarios l, weighted with their respective probability p(l).
Eqs. 10 -25 are the constraints of the second-stage problem. Here, Eq. 10 couples first and second stage by restricting the number of activated pumps in each load scenario by the number of pumps purchased in the first stage. Eqs. 11 -13 model the operation ranges of the different pump types available in the construction kit. For each pump type t, Eqs. 15 and 17 sum volume flow and power consumption of the number N t of pumps of this type installed. Eqs. 14 and 16 sum the total volume flow and power consumption of all pumps in the system.
The characteristic curves of each variable speed pump of type t are given by Eqs. 18 and 19, in which the fitting parameters a t,H/P , b t,H/P , c t,H/P , d t,P are determined in pre-processing based on the pumps' characteristic curves for their nominal rotating speed ω t,max . The flow rate of each pump type is limited by a linear restriction, cf. Eq. 20. The head and power consumption of each pump type are limited by Eq. 21-22. This limit is dependent on the flow rate and is fitted based on the pump type's working area, cf. Fig. 4. The minimum head of activated pumps is restricted by the load requirement in Eq. 23.
The chance constraint in Eq. 25 ensures that the volume flow demand q(l) of each scenario l is satisfied with a confidence level of at least 1 − ε. Small deviations of 1% of q(l) are allowed. Eqs. 24 and 25 are reformulated as Eqs. 26 -28 to integrate them into a Mixed-Integer Nonlinear Program (MINLP).
Eq. 26 models that the actual water supply Q total (l) may deviate from the demand q(l) by ∆(l). If this deviation is bigger than 1% of q(l), we assume that there is severe water shortage in scenario l, cf. Eq. 27. We introduce a binary variable Z(l) which indicates this water shortage. It is coupled to the deviation ∆(l) in Eqs. 27-28. The term q max := max l∈L q(l) is calculated in pre-processing once the N -sample of scenarios is drawn, and thus a parameter within the MINLP. Eq. 29 counts the times of severe water shortage and restricts their number by the risk level ε multiplied by the number of sampled scenarios N .
s.t. : where X = {X 1 , . . . , X T , Y 1 , . . . , Y T } is the set of the first-stage variables, and E [F (X , q)] is the optimal expected value of the second-stage problem: where Eqs. 24, 25 are reformulated by 106 Uncertainty in Mechanical Engineering III  [q min,t , q max,t ] Q t (l) sum of flow rates of pumps of type t in load scenario l R + Q total (l) sum of flow rates of all pumps in load scenario l R + Ω t (l) rotational speed of pumps of type t in load scenario l [ω min,t , ω max,t ] H t (l) head of pumps of type t in load scenario l [h min,t , h max,t ] P pump t (l) power of single pump of type t in load scenario l R + P t (l) sum of power of pumps of type t in load scenario l R + P total (l) total power consumption of all pumps in load scenario l R + ∆(l) water shortage in scenario l in m 3 /h R + Z(l) binary indicator for water shortage in scenario l {0, 1}

Optimization Algorithm
The CCTS is solved by an algorithm adapted from [15]. We use this algorithm as a primal heuristic to find good candidate solutions for our nonlinear model. The single steps are depicted in the flow chart in Fig. 3 and are the following:

Initialization of parameters
Initialize sample sizes N 1 , N 2 , the number of iteration steps I, the risk level ε.
Repeat for i = 1, 2, ..., I (1) Generate a sample {q(1), q(2), ..., q(N 1 )} of the normally distributed water demand q and solve the respective SAA problem. The corresponding candidate solution, its firststage solution and objective value are denoted by S opt,i a bigger sample {q(1), q(2), ..., q(N 2 )} to validate the feasibility of the candidate solution. In a first step, check for each demand q(1), . . . , q(N 2 ) if it lies within the performance range of the candidate topology (i.e. could be fulfilled, regardless of the respective power consumption). The percentage of the N 2 demands that cannot be fulfilled by a certain topology X i N 1 is denoted by p N 2 (X i N 1 ). (3) To assess the feasibility of the candidate topology X i ) /N 2 , where 1.96 is the approximate value of the 97.5 percentile point of the normal distribution, i.e., the corresponding coefficient for a confidence level of 95%.
If the feasibility check is not passed, go to the next iteration. If it is passed, set the first-stage decisions according to X i N 1 and solve the remaining problem of minimizing the expected operation and penalty cost for N 2 scenarios.
(6) Denote the solution of this approximation by U i N 2 (X i N 1 ). Note that its objective is a lower bound for the objective value of the original CCTS problem.
(7) Pick the minimal U i N 2 (X i N 1 ) of all I iterations and denote it and its respective solution by U * and X * , respectively.

Computational Results and Discussion
We compute optimal topologies and operation parameters for different risk levels. The construction kit of pumps is given by Tab. 2, Fig. 4 shows their respective working areas and power consumption.  We restrict the maximum number of parallel pumps of the same type t by five, i.e. m min,t = 0 and m max,t = 5. We generate the volume flow demand scenarios based on a normal distribution, N (30 m 3 /h, 15 2 m 3 /h). The minimum head of each activated pump is required to be ∆h min = 3 m. The pump system is assumed to be working for s = 5 years with an operation cost coefficient of 0.3 e/kWh and penalty cost coefficient of 250 e/10 3 m 3 h −1 if the supplied volume flow deviates from the demand by more than 1 %. The algorithm is implemented in Matlab using the toolbox YALMIP [18]. SCIP version 5.0.1 [19] is used for solving the MINLPs.
Considering [14], [15], and the computational effort, the parameters of the algorithm are chosen as follows: N 1 = 20, N 2 = 100, I = 10, J = 10, and β = 0.01. Results for different risk levels ε are shown in Tab. 3: Here, "Weighted Water Shortage" is the probability-weighted proportion of water demand not covered by the optimal solution. It is calculated by 1 − E P [ min(q, Q total )/q ] . "Realizable Scenario Coverage" is the proportion of the 100 sampled scenarios, the optimal topology could fulfill, regardless of its power consumption and thus operation costs. The higher this value, the larger the system's performance range, cf. [20], and thus its resilience.
It can be found that the objective value decreases for increasing ε, since the chance constraint in Eq. 29 is less restrictive. For ε = 100 %, the chance constraint does not have to hold, and in order to save energy and investment cost, the weighted water shortage could rise up to 100 %. Yet, due to penalty cost, the optimal solution is to still keep water shortage around 26.2 %. Generally, for increasing risk level ε, the number of scenarios covered decreases. These overall lower costs for higher ε come with a price: The system's performance range, and thus its resilience decreases. Fig. 5 shows the optimized topologies for four different risk levels. For higher risk levels, it is beneficial to pay more penalty cost rather than to buy another or a bigger pump. Specifically, Fig. 6 shows the trade-off between investment, operation and penalty cost under different risk levels. The subfigure above depicts the first-stage investment cost, the two subfigures underneath show the second-stage objectives operation and penalty cost. For higher risk levels, it is advantageous to cut the cost of investment and operation, but to increase the expenditures on penalty cost. By varying ε, a system designer may find a cost-optimized system suitable for the risk level of his specific application.

110
Uncertainty in Mechanical Engineering III

C C
Optimal topologies for risk level of 50% and 100%. OUTLET

C C C
Optimal topologies for risk level of 10% and 20%.

Uncertainty in Mechanical Engineering III
A detailed analysis of the solution for risk level ε = 20 % is depicted in Fig. 7. Fig. 7 (a) shows the sampled water demands and their coverage by the optimized topology, i.e. its performance range. The scenarios which were not totally fulfilled are marked with lines. Fig. 7 (b) shows the respective deviation ∆ from the required demand. It is beneficial to not fulfill the demand in the lower region, since operation costs can be exchanged for relatively low penalty costs (note that the penalty cost increase linearly with ∆). For medium demands however, it is advantageous to fully provide the required volume flow, i.e., ∆ = 0. The operation costs for totally fulfilling a demand rise with the fluid power demanded. For two higher load scenarios, the demand is only partially satisfied. Overall, an optimal decision corresponds to a complex trade-off: Both the system's performance according to the required risk level has to be guaranteed, as well as balancing operation and penalty cost.
We investigate operating cost in Fig. 7 (c) which shows the power consumption of the system. We compare the power consumption for the optimized topology (where deviation from the load is allowed) to the power consumption which results, if the load is covered as good as possible. To calculate the power consumption for the second case, we set the first-stage variables according to the solution found, and minimize the deviation ∆. In both cases, we see a mostly linear relationship between power consumption and water demand. In general, the power consumption can be calculated by dividing the hydraulic power, which increases linearly with the volume flow, by the hydraulic efficiency of the pump system. The linear increase shows, that for the optimized topology and a minimal required head of 3 m, the system is able to operate with approximately constant efficiency for most demands. In the second case when covering the load as good as possible, a high relative power is needed to minimize the deviation for the lowest and highest demands. For these demands, it is beneficial to exchange operation costs for relatively low penalty costs. Fig. 7 (d) shows an analysis of the margin of the system in the different load scenarios. The margin of a technical system describes how close to its performance boundary the system is currently operating, cf. [20]. In our context, the system's margin is given by the distance of the actual operation points in the different scenarios to the maximum operating point possible. To investigate this distance, in a first step, we set the first-stage topology decision according to the found candidate solution (four pumps of type C for ε = 20 %) and maximize the volume flow that can be delivered by the system. Since we demanded a minimal head of 3 m, this maximal volume flow is 60 m 3 /h, cf. Fig. 4. In a second step, we plot the margin for each sampled load scenario. We can see that, except for the highest demand, all generated samples are within the performance bounds of the system.

Conclusions and Outlook
In this work, a nonlinear chance-constrained two-stage stochastic program (CCTS) was developed to optimize booster stations operating under uncertain load. The chance constraint prevents overdimensioning. It models that not all possible water demand scenarios have to be strictly fulfilled, but may be violated with a predefined probability. This allows to balance investment, water shortage and energy costs and to engineer the system's performance range and thus its resilience. To find optimized system layouts, a tailored algorithm for solving CCTS programs in the context of unit commitment problems was adapted for pump system design. For different risk levels, booster station topologies were optimized, and their cost analyzed. While the approach was used to optimize parallel system structures, it can be extended to non-parallel ones by adjusting the underlying MINLP. Since costs for water shortage and energy may exceed the initial investment by far, a pump system's optimal operation is of key importance. Thus, investing in technologies for pump failure detection [21] may be beneficial. We plan to integrate a cost/benefit analysis of such approaches in future optimization models.