Optimization and validation of pumping system design and operation for water supply in high-rise buildings

Abstract

The application of mathematical optimization methods for water supply system design and operation provides the capacity to increase the energy efficiency and to lower the investment costs considerably. We present a system approach for the optimal design and operation of pumping systems in real-world high-rise buildings that is based on the usage of mixed-integer nonlinear and mixed-integer linear modeling approaches. In addition, we consider different booster station topologies, i.e. parallel and series-parallel central booster stations as well as decentral booster stations. To confirm the validity of the underlying optimization models with real-world system behavior, we additionally present validation results based on experiments conducted on a modularly constructed pumping test rig. Within the models we consider layout and control decisions for different load scenarios, leading to a Deterministic Equivalent of a two-stage stochastic optimization program. We use a piecewise linearization as well as a piecewise relaxation of the pumps’ characteristics to derive mixed-integer linear models. Besides the solution with off-the-shelf solvers, we present a problem specific exact solving algorithm to improve the computation time. Focusing on the efficient exploration of the solution space, we divide the problem into smaller subproblems, which partly can be cut off in the solution process. Furthermore, we discuss the performance and applicability of the solution approaches for real buildings and analyze the technical aspects of the solutions from an engineer’s point of view, keeping in mind the economically important trade-off between investment and operation costs.

1 Introduction

More than half of the world’s population (55 %) lived in urban areas in 2018, as shown in UN (2018). Based on this data published by the United Nations, it is estimated that 2.5 billion people will be added to the urban population by 2050 leading to an increase by more than half of the number of people living in urban areas today. This continuous trend towards urbanization will lead to an increasing number of mega-cities, cities with more than 10 million inhabitants.

As the space for living in cities is limited, one effective way to cope with the steady urbanization and to create living space are high-rise buildings. In those buildings distributed pressure boosting pumping stations, also called booster stations, provide water for every floor, especially in higher pressure zones. Consequently, their need will increase proportionally with the number of high-rise buildings. Already today there is a rapid growth of completed buildings of 200 meters or greater height, CTBUH (2018).

Water supply systems are a key component of any urban infrastructure, cf. Coelho and Andrade-Campos (2014), and urbanization will pose increasing challenges for these infrastructures in the future: Whilst a secure supply of water for the growing number of citizens is to be ensured, energy consumption for water supply must be reduced on the way to a more sustainable future. Thus, efficiency optimization of water supply systems has been in the focus of ongoing research, for a review cf. e.g. Coelho and Andrade-Campos (2014).

While existing research has so far tended to focus on the energy consumption of larger centralized water infrastructures, in this contribution, we want to investigate the energy consumption of the mentioned booster stations for buildings. Especially in tall buildings, booster stations have a considerable energy consumption, in particular if they are wrongly dimensioned and/or operated. Therefore, the potential for energy savings through the application of several pressure zones or a decentralized pump placement and an optimal selection of pump types is high, cf. e.g. Altherr et al. (2019).

In order to unlock this optimization potential, a large number of pipe network topologies and pump configurations must be considered in the planning process. However, with the increasing number of degrees of freedom, especially in tall buildings, it is a difficult challenge to investigate all possible system variants, as the number of possible design choices can be overwhelming. Moreover, one not only has to check whether a promising design is able to fulfill the different load scenarios occurring in the building, but also for an economic solution by balancing investment and operating costs against each other. To achieve these goals, we present an optimization-aided system design approach based on a problem-specific solving algorithm. In order to ensure that the underlying optimization models represent the system under consideration with sufficient accuracy, we also present an experimental validation. To model the pressure losses due to friction we use own measurements as well as literature data and compare the optimization results of both approaches.

In the following Sect. 2, we introduce relevant literature for the optimization. Afterwards, we present the technical application in detail and show different formulations of the optimization model as well as a problem specific solving algorithm. We then discuss the performance of different solving approaches as well as technical aspects of the solutions and close with an experimental validation. Finally we draw a conclusion and give an outlook for future work.

2 Related work

In this article we present a holistic approach to the optimal synthesis of technical systems, which we illustrate by means of a pumping system for water supply. In the context of this research focus, the state-of-the-art is presented in two different areas:

• (i) Methods used in engineering for the design and efficiency optimization of pumping systems, in particular pressure boosting systems in buildings.

• (ii) Application of mathematical optimization methods to general technical systems, in particular methods for two-stage optimization problems.

In the next two subsections we will give a short overview of relevant literature in these two areas.

2.1 Design and efficiency optimization of pumping systems

When optimizing the efficiency of a pumping system, a first step is to assemble the system from efficient components. For the transport of drinking water, booster stations consisting of centrifugal pumps are commonly used. Therefore, we focus on this pump type within this literature overview.

The development and hydraulic optimization of centrifugal pumps is already at a very high level and multiple standard works are focusing on centrifugal pump design as well as performance, cf. e.g. Japikse and Marscher (1997), Brennen (2011), or Gülich (2010). In recent decades, also through the use of computer-aided methods such as Computational Fluid Dynamics (CFD), cf. e.g. Garg (1998), the efficiency of centrifugal pumps, in particular at their design point, has been continuously improved.

In real life applications, however, it can be seen that centrifugal pumps are predominantly not operated at their design point, but in the partial load range for a large part of the time, cf. Hirschberg (2014). In this case, the excess energy is either dissipated via throttling valves, or—if variable speed pumps, so called RPM-regulated pumps, are employed—they are operated at a lower rotational speed than their nominal rotational speed. Both leads to an operation at lower efficiency, often far away from their best efficiency point. Therefore, when classifying and optimizing pump performance, a suitable approach for assessing efficiency in practical usage is crucial.

The product approach, which is part of an european standard for energy-efficiency evaluation of water pumps, c.f. EN 16480:2016 (2016), is a first step in this direction. This approach no longer considers merely the pump as a component, but the pump and electric motor as part of a system. The efficiency of this system is assessed according to the Minimum Energy Efficiency Index (MEI) which rates the pump’s efficiency not only at its design point, but also at partial load and overload. Therefore, the so-called “house of efficiency” scheme is applied, which sets efficiency criteria for 75%, 100% and 110% of the volume flow at the design point.

While this is a first step towards the right direction, when pumps are used in a real building, not only their volume flow varies, but so does the required pressure increase in relation to the demanded flow. In addition, when assessing overall energy efficiency, it is important to take the time slice of the respective volume flow and its corresponding pressure requirements into account. For this reason, the Extended Product Approach, cf. Stoffel (2015), was developed within the scope of energy efficiency guidelines of the European Union. It considers a load profile which links volume flow and pressure demands, as well as corresponding time profiles. The energy consumption of the Extended Product (pump and electric motor) considering these load profiles then yields the Energy Efficiency Index (EEI), cf. Lang et al. (2013). In this work, we present a method which considers typical load profiles within buildings and is in accordance with the EEI.

When building pumping systems, several technical aspects have to be considered. The standards DIN 1988-500 (2011) and DIN 1988-300 (2012) provide codes of practice for drinking water installations. While DIN 1988-500 (2011) focuses on the dimensioning of pressure boosting systems with RPM-regulated pumps, DIN 1988-300 (2012) describes the dimensioning of the pipe network. Particularly, the selection of suitable pipe diameters in regard to restrictions such as a peak load volume flow and a maximal water velocity for reasons of noise is outlined. In addition to a central pump configuration, the standard DIN 1988-500 (2011) also considers some decentralized pump configurations due to the specification of different pressure zones within the building. E.g., lower pressure zones of the building can be supplied directly by the public water mains and solely higher pressure zones are supplied via a pressure boosting system. Furthermore, the installation of several booster stations is described so that each pressure zone is supplied by its own pumping system. However, more sophisticated decentralized configurations, which could be automatically assessed and evaluated by our approach, are not considered.

Whereas efficient components are certainly a basic prerequisite for an efficient overall system, they are not yet a guarantee for this. With a system-wide assessment in mind, the scope of the international standard “Pump system energy assessment” is to set “the requirements for conducting and reporting the results of a pumping system energy assessment [...] that considers the entire pumping system, from energy inputs to the work performed as the result of these inputs”, ISO/ASME 14414 (2019). However, the standard neither specifies how to design a pumping system nor how to validate the energy savings that would result from implementing recommendations based on this assessment. That are two aspects that we cover in this work.

In general, there are several system-wide measures to increase the energy efficiency of pumping systems:

• The optimal selection of pumps, i.e. usage of variable-speed pumps, cf. Ferreira et al. (2010), and avoiding over-dimensioning, cf. e.g. Fisher (2001), Tindall and Pendle (2015), Weber and Lorenz (2017), Weber and Lorenz (2018).

• The allocation of several pressure zones within the building and the usage of decentralized pump layouts, cf. e.g. Norgaard and Nielsen (2010); Leise and Altherr (2018).

• The optimization of the pumps’ operation, cf. e.g. Pedersen and Yang (2008), Groß et al. (2017), and Nowak et al. (2018).

The approach presented in this contribution allows to cover all three of these aspects, since not only the optimal pumping system design, but also its optimal operation can be considered.

While in this work we particularly focus on the optimization of the efficiency and investment costs, other possible objectives are shown in the literature. One example is the optimization of resilient water supply networks as shown by Herrera et al. (2016), Hartisch et al. (2018), and Altherr et al. (2019).

2.2 Mathematical aspects of two-stage optimization problems in engineering

Besides literature from the engineering domain covering the technical aspects of optimizing pumping systems, also previous works regarding the mathematical aspects are relevant.

When designing an economically optimal pumping system, the overall goal often is to minimize the sum of investment costs and energy costs for the expected life cycle. Since the load within residential buildings varies, cf. Sect. 2.1, its operation at different load points has to be considered to assess the energy costs of the pumping system.

The corresponding mathematical optimization problem consists of two stages: First, finding a low-priced investment decision, i.e. a pump and pipe configuration. Second, operating a subset of the chosen pumps for each quasi-stationary demand scenario of the given load profile, such that the system satisfies this demand and at the same time the efficiency is maximized.

In this case, the considered load profile can be seen as a discrete distribution of uncertain load parameters, with the time portion of each load scenario corresponding to its probability of occurrence. The optimization program, which minimizes investment costs as well as expected value of operation costs, can then be seen as the Deterministic Equivalent of a general two-stage stochastic program with recourse, cf. e.g. Wets (1974).

As many optimization problems in the area of engineering design are subject to uncertain parameters, e.g. uncertain material characteristics or load scenarios, two-stage stochastic programs have been frequently studied in this context, and a good overview can be found in Popela (2010) and Popela et al. (2014). The life cycle of products and systems in engineering includes the two phases of planning (first stage) and operation (second stage). In two-stage stochastic optimization problems the operation is already anticipated in the planning phase, which leads to better investment decisions. The operation itself is then typically ensured by closed- or open-loop controllers.

An important aspect for the difficulty of solving these two-stage optimization problems to global optimality is whether integer variables are (1) present in the first stage solely, or (2) also present in the second stage. Regarding design and operation optimization of technical systems, integer variables in the first stage, and more precisely binary variables, are often used to model purchase decisions. Moreover, in this work, binary variables in the second stage are used to model the activation or deactivation of active components, i.e. an on/off-switch, leading to problems of case (2).

Theoretical aspects of these problems, also called complete mixed-integer recourse problems, have e.g. been investigated by Schultz (1992), and sophisticated solution approaches have been proposed, cf. e.g. Carøe and Tind (1998). Given the formulation as Deterministic Equivalent Program, standard solution approaches based on Branch & Bound with LP relaxation may be applied. Yet these might suffer from high computation times, especially if multiple (load) scenarios are considered in the second stage.

3 Technical application

As mentioned in the previous sections, we generate the Deterministic Equivalent Programs of two-stage stochastic optimization models to compute energy-efficient centralized and decentralized booster stations in high-rise buildings and validate these models on a test rig. The physical properties of these systems and topology decisions are modeled within the constraints of the optimization programs. Within this contribution, the fresh cold water supply is considered solely, while waste water, hot water, and water for firefighting are not considered. Thereby, the objective is to minimize the simplified life-time cost, which consists of investment and energy costs within the lifetime. These two parts represent the economically most important goals, cf. Tolvanen (2007), when designing pumping systems. When modeling these systems, two aspects are of importance. First, reasonable computation times of model instances are required. Second, the model’s level of detail has to represent the real system with an acceptable accuracy to allow for the transfer of the optimization results. The computed optimum within the model instance, should represent the true optimum of the underlying technical system. Usually, a trade-off between these two goals is necessary, as the former tends to lead to models that are as simple as possible, while the latter demands for models that are highly accurate, but computationally not tractable. To achieve a model that is as simple as possible and at the same time as accurate as needed, we conduct an experimental validation step based on a test rig that represents a scaled building with five pressure zones.

We present all relevant basics of the pump and booster station modeling, the test rig, the system modeling and the preselection of pumps in the following subsections.

3.1 Pump modeling

The most important components to fulfill the water supply in high-rise buildings are the pumps. All pumps are modeled with polynomial approximations for the pressure–volume flow characteristic, as well as the power–volume flow characteristic. This approach is based on Ulanicki et al. (2008). We model the pressure increase $\varDelta p = f(q,n)$ with $f:{\mathbb {R}}^{+}_{0}\times {\mathbb {R}}^{+}_{0} \rightarrow {\mathbb {R}}^{+}_{0}$ by using a quadratic approximation relating the pressure head $\varDelta p$, the volume flow $q \in {\mathbb {R}}^{+}_{0}$ and the normalized rotational speed $n \in [{\underline{N}},{\overline{N}}]$, as given by

$$\begin{aligned} \varDelta p= \sum _{m=0}^{2} \alpha _{m} \,{q^{\mathrm {}}_{}}^{2-m} \, {n^{\mathrm {}}_{}}^m, \end{aligned}$$

(1)

where the parameters $\alpha _{m}$ and the lower bound of the normalized rotational speed ${\underline{N}}$ depend on the pump type. The rotational speed is normalized by using its maximum value. Therefore, we get an upper bound ${\overline{N}}=1$. For better readability, we refer to the normalized rotational speed in the following as rotational speed or speed, only. The power $po = g(q,n)$ with $g:R^{+}_{0}\times {\mathbb {R}}^{+}_{0} \rightarrow {\mathbb {R}}^{+}_{0}$ is modeled by using a cubic approximation,

$$\begin{aligned} po = \sum _{m=0}^{3} \beta _{m} \,{q^{\mathrm {}}_{}}^{3-m} \, {n^{\mathrm {}}_{}}^m. \end{aligned}$$

(2)

The parameters $\alpha _{m}$ and $\beta _{m}$ are derived from manufacturer provided component characteristics for each pump type. Exemplified characteristic diagrams for a pump are shown in Fig. 1.

Fig. 1

Characteristics for an example pump based on the used pump model. These characteristics are usually limited in each domain according to specific pump requirements. The different curves describe several different rotational speeds

The shown characteristics are often limited by maximum power as well as minimum and maximum volume flow requirements for each pump. These limits are added in the underlying optimization as additional constraints to model a realistic pump behavior. For simplicity they are not shown in Fig. 1.

3.2 Booster station modeling

The function of a booster station is to provide sufficient pressure to overcome the geodetic height of the floor to be provided with water as well as the friction in the pipes feeding the respective floor. We distinguish three different booster station topologies, which differ in the alignment of the considered pumps. Today, pumps are commonly aligned in parallel at the lowest pressure zone of the building, as the product range of various manufacturers (Wilo 2020; KSB SE Co KGaA 2020; Grundfos Pumps Ltd 2016) and the discussion in (Gülich 2010, pp. 666–669) shows. This leads to a central, parallel booster station design providing the pressure and volume flow for the whole building. Varying from this approach, also the serial arrangement of pumps in a central booster station is possible, leading to higher degrees of freedom in the design process and therefore allowing for higher energy efficiency of the booster station. For the centralized topology the building’s demand can be modeled as solely one sink with one available water source. Additionally to these centralized system topologies, a decentralized booster station topology is possible. In this approach, pumps can be positioned at any pressure zone of the building and therefore a more accurate modeling of the volume flow demand, i.e. the flow demand per pressure zone, is necessary. In the following, we comprise floors of similar height into individual pressure zones. Each pressure zone is modeled as an individual sink which has a different pressure demand, depending on the height of the floors. The decentralized topology allows for an even higher degree of freedom of pump arrangement and a reduction of the throttle losses, facilitating the improvement of the overall efficiency even further.

Fig. 2

Exemplified system design of the underlying optimization program

The system design approach with an exemplified booster station is shown in Fig. 2. In this exemplary system three pumps are used. The first pump is used in series with the remaining two pumps. The second and third pump supply different pressure zones in the building individually.

3.3 Test rig for validation

As indicated before, an experimental validation of the optimization results is crucial to guarantee the technical functionality. A validation using real buildings is not possible for reasons of high effort. In addition, novel concepts, i.e. the decentralized booster station design, shall be investigated, which cannot be realized in already existing buildings.

Fig. 3

Sketch (left) and photo (right) of the modular test rig with a central booster station design

For this reason, a down-scaled, modular test rig has been set up, cf. Figure 3 and Müller et al. (2019b). The down-scaled system’s booster station has the same characteristics as those of high-rise buildings. The maximum geodetic height difference is 5 m divided into five different pressure zones, which can be supplied individually. The water is pumped from a tank with constant pressure to the five pressure zones by the booster station. As in real buildings, the water flows in at a constant supply pressure and flows out in an open drain pipe at ambient pressure (open system). The water returns to the tank to keep the water level and therefore the supply pressure constant. The test rig is equipped with sensors for measuring the volume flow in each pressure zone as well as the pressure head and power consumption of each pump. It is possible to control the speed of the pumps and the valve position on each pressure zone and thus to control the volume flow on each pressure zone. Overall, a set of thirteen pumps, consisting of six different types, is available to build and test different booster station designs.

The pumps are centrifugal circulator pumps for heating systems with an integrated frequency converter (FC) and therefore speed controllable. The sizes are chosen to fit the needs for the down-scaled booster system, even though the ratio of pressure head to volume flow is lower than in real-world booster station pumps. The pumps can be connected in any way and up to six pumps can be operated simultaneously as a booster station for the test rig. Compared to real pressure booster stations, the friction caused by installation fittings of the pumps (e.g. for pressure measurement) on the test rig is quite high. Therefore, those pressure losses are taken into account in the loss coefficient $\zeta _\mathrm {inst}$.

3.4 System modeling and load scenarios

We consider three different real-world high-rise hotel buildings in addition to the previously described test rig as examples for the presented system optimization approach. We consider each building’s geodetic height, base area and type to estimate the volume flow demand and the respective pressure increase. The estimation of volume flow demands is highly uncertain if no consumption patterns are known previously to the optimization. An overview of water demand for hotel buildings is given by Gössling et al. (2012). Besides these different measurements, Hirschberg (2014) also presented an approximation method to estimate the volume flow demand in buildings of different types. Another possibility to approximate the volume flow in buildings is given by DIN 1988-300 (2012). It shows a computational approach to estimate the volume flow demand in buildings based on a standardized estimation procedure.

Table 1 Characteristics of the chosen real-world hotel buildings and the down-scaled test rig

We estimated the volume flow making use of the estimation method shown by Hirschberg (2014) for building 1 and 2. The volume flow demand in building 3, as shown in Table 1, is based on real consumption data in this building. In the following, the different buildings are referred to as B1, B2 and B3, respectively. We assume that all buildings in Table 1 are divided into 5 pressure zones similarly to the test rig.

The demand of the test rig system is derived by scaling the demands of real buildings according to the geodetic height and the friction losses of the test rig.

The maximum load scenario—given by the maximum volume flow demand ${\overline{Q}}$ and maximum pressure head $p_\mathrm {s}$—allows to lay out the booster station for the peak load. This load scenario occurs only for a fraction of the operation time. To design the booster station more realistically, we also consider the events of partial load. The volume flow demand in partial load is derived according to Hirschberg (2014) and all scenarios are shown in Fig. 4.

Fig. 4

Load scenarios based on Hirschberg (2014) and their according time shares

The necessary pressure at the sink of a central booster station is calculated using the geodetic height of the building H, the minimum flow pressure $p_\mathrm {flow} = 1\,\hbox {bar}$, the supply pressure $p_\mathrm {sup} = 3\,\hbox {bar}$ (i.e. provided by the connection to the regional water supply system) and the friction, which depends quadratically on the flow rate. According to Hirschberg (2014), the pressure loss due to friction in the maximum load case is approximately $\sigma = 52 \, \%$ of both the geodetic height and minimum flow pressure. This leads to:

$$\begin{aligned} \begin{aligned} p_s&= \varrho g H + p_\mathrm {flow}+ \frac{1}{2} \varrho \zeta \left( \frac{Q_s}{A} \right) ^2 - p_\mathrm {sup} \\&= (\varrho g H + p_\mathrm {flow}) \left( 1+\sigma \left( \frac{Q_s}{\bar{Q}} \right) ^2\right) - p_\mathrm {sup} \quad \forall s \in {\mathcal {S}}. \end{aligned} \end{aligned}$$

(3)

For the decentralized system design, we derive individual load scenarios for each pressure zone, which sum up to the load scenarios in Fig. 4. The pressure at the different sinks is derived by scaling the pressure for the central booster station according to the geodetic height and the volume flow of the individual pressure zones.

3.5 Preselection of pumps

In this work, we apply a model series of 200 different market-available pumps which were explicitly designed for the water supply of buildings. The combinatorial possibilities resulting from this wide range of variants is challenging within reasonable computation time. In addition, most pumps are not suitable for each specific building. Therefore, we aim to reduce the number of considered pumps based on domain-specific knowledge. To the best of the authors knowledge there exists no general purpose method to select multiple appropriate pumps for consideration within an optimization program. Sets of pumps for the respective buildings are currently only selected based on expert knowledge as for example shown by Larralde and Ocampo (2010a) and Larralde (2010). Therefore, a crucial step in the optimal design of booster stations for high-rise buildings is the selection of an appropriate set of available pumps, as shown by Leise and Altherr (2018). For this, we implemented a heuristic, as shown in Algorithm 1, for the pump selection, which is based on a two-step process.

Fig. 5

Characteristic diagrams of all used pumps at their maximum rotational speed for a building B1 of 100 m height, b building B2 of 80 m height, c building B3 of 40 m height, and d the test rig . The considered scenarios are shown as black squares. The scaled points to select pumps in Algorithm 1 are shown as gray circles for the three buildings

In the first step of the heuristic selection of pumps we define additional partial load scenarios that should be met by individual pumps selected in the following step. These additional scenarios are derived from the original scenarios (represented by black square markers in Fig. 5) by taking advantage of the generally valid relationship of pressure and volume flow for in series and parallel connected pumps:

On the one hand, for in series connected pump systems with N pumps the volume flow in each pump is equivalent to the total system volume flow, while the pressure in each pump is given by $\varDelta p/N$ (if we assume that the pressure increase is equal in each pump).

On the other hand, a parallel usage of N pumps of the same type results in total in the same pressure increase for each pump but the volume flow for each pump reduces to q/N if all parallel pumps are operated with the same speed.

These two effects for in series and parallel system design result in a reduction of the original load scenarios to fractions of the original load scenarios in q and/or $\varDelta p$ direction.

In the selection heuristic we start by generating all fractional loads that can be derived for parallel combinations with up to 6 equivalent pumps based on the original load scenarios. Afterwards we then derive all loads that are based on in series connections of up to 6 equivalent pumps for the originally given loads and newly computed ones.

For each building we want to consider a predefined number of pumps (in our case 15 in total). As the number of fractional loads in the first step exceeds this predefined number and also results in a clustering of fractional load scenarios with small power, we use a set distance $d^{\text {max}}$ in the $\varDelta p$–q domain to select only loads which have a predefined distance from each other. We then iteratively delete load scenarios with a lower hydraulic power demand within the range of the set distance. With this reduction step, we assure that every deleted load scenario point remains feasible in the solution, as we could use a pump with a higher hydraulic power and a frequency converter to reduce the speed according to the desired fractional demand. Additionally, it would be possible in a real building to throttle the flow to achieve the desired partial load scenario. This simple procedure ensures to derive the predefined number of pumps as well as a great coverage of the overall domain.

In the second step of the heuristic selection of pumps we select an appropriate pump for each of the remaining load scenarios in the $\varDelta p$–q domain which are given by gray circles and black squares in Fig. 5. A pump is selected if it can fulfill the desired demand at one fractional load and has the lowest electric power of all considered pumps at full speed. In this manner we derive a final pump catalog that consists of a predefined number of pumps (in our case 15), which are suitable to supply each original load scenario by a single pump or a combination of in series/parallel connected pumps. Additionally, each selected pump has the lowest power demand to supply the partial load it was selected for. By implementing and applying this heuristic, we remove all pumps that are not suitable and all suitable pumps with a higher power demand for the supply of a respective building. In the following, we consider solely a set number of preselected pumps which specifically supply a certain building efficiently.

4 Optimization model

This section deals with the mathematical formulation of the previously introduced technical design problem. The technical specifications are modeled in a basic stochastic optimization model (Sect. 4.1).

We reduce the two-stage stochastic problem to its Deterministic Equivalent Program by using multiple load scenarios. In our case, the first stage variables represent investment and design decisions while second stage variables represent control decisions. Within this scenario based optimization approach binary variables are present in both stages. Even in the second stage as the model takes into account that not all pumps are always switched on to meet the demand of the different load scenarios. In contrast, some pumps may also be switched off and bypassed for efficiency reasons, especially in load cases with lower volume flow requirements.

Since pumps have nonlinear characteristics, the mixed-integer nonlinear program (MINLP) pump model (Sect.  4.2) as well as two linear pump models, allowing the formulation of a mixed-integer linear program (MILP), are presented. In a first approach, a piecewise linearized pump model (Sect. 4.3) is derived and in a second approach, a piecewise linear relaxed pump model (Sect.  4.4) is established.

4.1 Basic optimization model

All derived optimization models in this contribution have a common mixed-integer linear foundation, which we present in this section. The complete system design with all possibilities is modeled as a graph ${\mathcal {G}} = ({\mathcal {V}}, {\mathcal {E}})$, in which the edges ${\mathcal {E}}$ model the pumps as well as the pipes. The junctions of the booster station are modeled as vertices ${\mathcal {V}}$. The basic model’s variables, which are indicated by lower case letters, as well as sets and parameters, indicated by upper case and Greek letters, are shown in Table 2.

Table 2 Decision variables, sets and parameters

The shared objective of all presented optimization models is the maximization of the net present value, cf. Meck et al. (2019). This results in a minimization of lifetime cost, which are approximated considering investment costs for pumps and frequency converters as well as discounted energy costs. The objective is shown in Obj. (5a). We use the annuity-present value factor $APF$ to discount all costs to the present day.

We model an upper bound of pumps, Constr. (5b), to solely select a suitable subset of all possible pumps in ${\mathcal {P}}$ in the final solution, i.e. maximum $N^\mathrm {buy}$ pumps may be bought. The necessary precondition to use a pump is that it was bought, Constr. (5c). Purchased pumps should be used only in their physical domain. Therefore, we limit the pump speed, power and volume flow with Constr. (5d) and (5e). Constr. (5f) models the volume flow conservation in the network, while Constr. (5g) models upper bounds for the volume flow in each pipe. We demand the volume flow in the sinks to be exactly the predefined volume flow for each scenario, as shown in Fig. 4 with Constr. (5h). It is imperative to choose at least one pump to connect the source with the sinks, Constr. (5i) and (5j). Constr. (5k)–(5m) prevent undesired topologies, as for example loops. The pressure at the inlet of a pump is set to the pressure at the source $P^{\mathrm {source}}_s$ in Constr. (5n) if a pump is connected to the source. The pressure increase from the inlet to the outlet of a pump is modeled by Constr. (5o). The pressure at the inlet and outlet of an unused pump is set to zero, Constr. (5p). The pressure of a pump’s outlet should fulfill the predefined pressure requirements at the sink if it is connected to the sink, Constr. (5q). The pressure propagation between connected pumps is modeled by Constr. (5r). If a frequency converter for a certain pump is not purchased ($y^{FC}_i=0$), the speed of the pump has to be set to the nominal speed ($n=1$) in case of activity, Constr. (5s). If the pump is turned off, or a frequency converter is purchased, this constraint has to be deactivated, Constr. (5s). If a frequency converter is bought, the pump is speed controllable, but a part of the transmitted electrical power is dissipated, Constr. (5t). In addition, the electrical power for the pump motor is added to the total power, Constr. (5u). Moreover, we can say that the electrical power of the system has to be higher than the hydraulic power divided by the best efficiency of the best pump, Constr. (5v). This provides a lower bound for the power consumption, which significantly speeds up the optimization process, as shown by Müller et al. (2019a).

In the subsequently introduced models, we allow for the option to restrict the usage of in series connected pumps, as pumps connected in parallel correspond to the current standard installation procedure, cf. Sect. 3.2. A solely parallel connection can be achieved by adding the following constraint to the basic model (5):

$$\begin{aligned} t_{i,j} = 0&\quad \forall \, i,j \in {\mathcal {P}}. \end{aligned}$$

(4)

The objective as well as all so far presented constraints, Model (5) and Constr. (4), depend linearly on the decision variables. The according mixed-integer linear program describes the generic selection and usage of pumps to supply multiple sinks with relevant connected pumps to derive a suitable booster station. Nevertheless, this model does not contain the explicit modeling of the used pumps, as shown in Eq. (1) and Eq. (2). Therefore, we present different methods to model these nonlinear relations of type $\varDelta p =f(q,n)$ and $po = g(q,n)$ in the following subsections.

4.2 Mixed-integer nonlinear pump model

The most accurate modeling approach is given by a nonlinear approach in which each considered pump is modeled with nonlinear constraints according to Eq. (1) and Eq. (2). We add all pumps by adding

$$\begin{aligned} \varDelta p^{\mathrm {}}_{i,s} = (\alpha _{i,0} - \zeta ^\mathrm {inst.}) \,{q^{\mathrm {}}_{i,s}}^{2} + \sum _{m=1}^{2} \alpha _{i,m} \,{q^{\mathrm {}}_{i,s}}^{2-m} \, {n^{\mathrm {}}_{i,s}}^m&\quad \forall i\in {\mathcal {P}}, s \in {\mathcal {S}}, \end{aligned}$$

(6a)

$$\begin{aligned} po^{\mathrm {pump}}_{i,s} = \sum _{m=0}^{3} \beta _{i,m} \,{q^{\mathrm {}}_{i,s}}^{3-m} \, {n^{\mathrm {}}_{i,s}}^m&\quad \forall i\in {\mathcal {P}}, s \in {\mathcal {S}}, \end{aligned}$$

(6b)

to the basic model in Eq. (5), where $\zeta ^\mathrm {inst.}$ is the pressure loss factor of the fittings used for installing the pump into the system. The benefits of this approach are a fast implementation and a highly accurate model, while the drawbacks are potentially high solution times.

4.3 Piecewise linearized pump model

A common approach to solve nonlinear non-convex optimization models with integer variables, is to make use of approximations based on piecewise linearization (PWL) techniques, cf. Geißler et al. (2012). Linearization of all nonlinear constraints results in a mixed-integer linear program, which can then be solved by state-of-the-art solvers as for example GUROBI (GUROBI 2020) or CPLEX (IBM 2020). Piecewise linearization methods are therefore applied in a high variety of applications, both within classical operations research, as well as technical applications, e.g. Misener et al. (2009), Morsi et al. (2012), Rausch et al. (2016) and Mikolajková et al. (2018). The underlying modeling to reformulate a nonlinear functional relationship as linear constraints can be implemented by multiple methods, as shown by Misener and Floudas (2010), Vielma et al. (2010) and Geißler et al. (2012).

Table 3 Additional variables, sets and parameters for PWL

In this contribution, we use the convex combination method, as described by Vielma et al. (2010), to model the nonlinear relations of the pressure head $\varDelta p = f(q,n)$ and of the power $po = g(q,n)$ as piecewise linearized function approximations. All additional variables, sets and parameters for the PWL are shown in Table 3. Instead of using Eq. (6) to model the nonlinear pump characteristics, we model a piecewise linearized approximation by using the set of constraints (7).

The index sets ${\mathcal {Q}}$ and ${\mathcal {N}}$ represent grid points in the (q, n) domain in which the pressure $\varDelta P_{u,v,i,s}$ and power $Po^{\text {pump}}_{u,v,i,s}$ are evaluated for each pump $i\in {\mathcal {P}}$ before optimization according to Eqs. (1) and (2). The pressure losses and additional power for frequency converters are added to these values in the preprocessing of the optimization to derive an equivalent description, as shown in Constr. (6). Besides these grid points, we have to use additional binary variables $a_{i,s,k} \in \{0,1\}$ and continuous variables $\lambda _{u,v,i,s} \in [0,1]$ to ensure a piecewise linearized approximation. The set ${\mathcal {K}}$ in Constr. (7f) contains all segments of the piecewise linearized approximation. In turn, the subset ${\mathcal {K}}(u,v, i,s)$ in Constr. (7g) covers all segments that contain the grid point (u, v) for pump i and scenario s. These additional constraints are added to the optimization program in model (5), which results in a mixed-integer linear program that can be solved by common MILP solvers.

4.4 Piecewise linear relaxed pump model

As a third possibility we use piecewise linear relaxed pump models within the final optimization program. Relaxations are often used in the MINLP solution process, as shown by Puranik and Sahinidis (2017), to derive valid lower bounds. A general introduction to relaxation methods for nonlinear network design problems is given by Humpola and Fügenschuh (2015). Relaxation methods with piecewise linear underestimators are also used in MINLP solvers to improve the solution speed, cf. Misener and Floudas (2012), Misener and Floudas (2013), Misener and Floudas (2014). Besides these multi-purpose implementations, these methods can also be used to improve the solution process for domain specific optimization models, like for example shown by Pecci et al. (2019) in the domain of water supply networks. The authors derived global bounds on the placement of valves by using a Branch-and-Bound algorithm with relaxed nonlinear constraints. As shown in all previous references, the usage of piecewise relaxed underestimators can improve the solution process. For this reason, we developed a relaxed pump model to improve the solution speed of the solving algorithm, cf. Sect. 5.

In order to underestimate the total cost it is necessary to overestimate the pressure increase and to underestimate the power consumption of the pumps:

$$\begin{aligned} { \varDelta p_{i,s}} \le \varDelta p_{i,s}^\mathrm {relaxed}&\quad \forall i \in {\mathcal {P}},s \in {\mathcal {S}}, \end{aligned}$$

(8a)

$$\begin{aligned} { po^\mathrm {pump}_{i,s}} \ge po_{i,s}^\mathrm {relaxed}&\quad \forall i \in {\mathcal {P}},s \in {\mathcal {S}}. \end{aligned}$$

(8b)

For these estimations, piecewise defined planes are determined, which serve as linear under- or overestimators. They have to be defined piecewise, since $\varDelta p =f(q,n)$ is not concave in n and $po^\mathrm {pump} = g(q,n)$ is not convex in q. The functions are shown in Figs. 6 and 7 along with the over- and underestimations.

Fig. 6

Piecewise linear overestimation of the quadratic function of the pressure $\varDelta p =f(q,n)$. Shown are a the function $\varDelta p(q)$ for different speeds n and b the function $\varDelta p(n)$ for different volume flows q. $R=6$ elements in q direction and $K=2$ elements in n direction are used. The dashed lines in a show the values of q for which the lines were plotted in b and vice versa

Fig. 7

Piecewise linear underestimation of the cubic function of the power consumption $po =g(q,n)$. a shows the nonlinear function and b the piecewise underestimation. c shows the function po(q) for different speeds n and d po(n) for different volume flows (q). $R=4$ elements in q direction and $K=4$ elements in n direction are used. The dashed lines in c show the values of q for which the lines were plotted in d and vice versa

We define R equidistant segments in q- and K equidistant segments in n direction. The boundaries of these segments are described by the sets ${\mathcal {Q}}^\mathrm {seg} = \lbrace {Q}^\mathrm {seg}_1,...,{Q}^\mathrm {seg}_{1+R} \rbrace$ and ${\mathcal {N}}^\mathrm {seg}= \lbrace {N}^\mathrm {seg}_1,...,{N}^\mathrm {seg}_{1+K} \rbrace$. For each segment (u, v) and each function ($\varDelta p =f(q,n)$, $po^\mathrm {pump} =g(q,n)$) a plane is defined, which should approximate the function as good as possible. Each plane is defined by the parameters $\gamma _{i,k,u,v}$ ($\varDelta p =f(q,n)$) and $\delta _{i,k,u,v}$ ($po^\mathrm {pump} =g(q,n)$) with $i\in {\mathcal {P}}$, $k\in \lbrace 1,2,3 \rbrace$, $u \in \lbrace 1,...., R \rbrace$ , $v \in \lbrace 1,...,K \rbrace$.

To solve the overall MINLP, cf. Sect. 4, the constraints of the pump characteristics (6) are replaced by the following constraints:

$$\begin{aligned}&\varDelta p_{i,s} - \left( \gamma _{i,1,u,v} \, {q_{i,s}} + \gamma _{i,2,u,v} \, {n_{i,s}} + \gamma _{i,3,u,v} \right) \le (1-x^{\gamma }_{i,s,v}){\overline{M}}_1 \nonumber \\&\forall u \in \lbrace 1,..., R \rbrace , v \in \lbrace 1,..., K \rbrace , i \in {\mathcal {P}}, s \in {\mathcal {S}} \end{aligned}$$

(9a)

$$\begin{aligned}&po^\mathrm {pump}_{i,s} - \left( \delta _{i,1,u,v} \, {q_{i,s}} + \delta _{i,2,u,v} \, {n_{i,s}} + \delta _{i,3,u,v} \right) \ge - (1-x^{\delta }_{i,s,u}){\overline{M}}_2 \nonumber \\&\forall u \in \lbrace 1,..., R \rbrace , v \in \lbrace 1 ,..., K \rbrace , i \in {\mathcal {P}}, s \in {\mathcal {S}} \end{aligned}$$

(9b)

$$\begin{aligned}&x^{\gamma 1}_{i,s,v} + x^{\gamma 2}_{i,s,v} \le x^{\gamma }_{i,s,v} \quad \forall v \in \lbrace 1,..., K \rbrace , i \in {\mathcal {P}}, s \in {\mathcal {S}} \end{aligned}$$

(9c)

$$\begin{aligned}&n_{i,s} \ge ({\mathcal {N}}^\mathrm {seg}_v + \varepsilon ) (1 - x^{\gamma 1}_{i,s,v}) \quad \quad \forall v \in \lbrace 1,..., K \rbrace , i \in {\mathcal {P}}, s \in {\mathcal {S}} \end{aligned}$$

(9d)

$$\begin{aligned}&n_{i,s} - ({\mathcal {N}}^\mathrm {seg}_{v+1} - \varepsilon ) \le (1 - x^{\gamma 2}_{i,s,v}) \quad \forall v \in \lbrace 1,..., K \rbrace , i \in {\mathcal {P}}, s \in {\mathcal {S}} \end{aligned}$$

(9e)

$$\begin{aligned}&x^{\delta 1}_{i,s,u} + x^{\delta 2}_{i,s,u} \le x^{\delta }_{i,s,u} \quad \forall u \in \lbrace 1,..., R \rbrace , i \in {\mathcal {P}}, s \in {\mathcal {S}} \end{aligned}$$

(9f)

$$\begin{aligned}&q_{i,s} \ge ({\mathcal {Q}}^\mathrm {seg}_u + \varepsilon ) (1 - x^{\delta 1}_{i,s,u}) \quad \forall u \in \lbrace 1,..., R \rbrace , i \in {\mathcal {P}}, s \in {\mathcal {S}} \end{aligned}$$

(9g)

$$\begin{aligned}&q_{i,s} - ({\mathcal {Q}}^\mathrm {seg}_{u+1} - \varepsilon ) \le {\overline{Q}}(1 - x^{\delta 2}_{i,s,u}) \quad \forall u \in \lbrace 1,..., R \rbrace , i \in {\mathcal {P}}, s \in {\mathcal {S}} \end{aligned}$$

(9h)

Constraints (9a) and (9b) guarantee that the pressure in each segment is underestimated and the power overestimated. Here, $x^{\gamma }_{i,s,v}$ and $x^{\delta }_{i,s,u}$ are binary activation variables which activate a plane if and only if the speed (in case of pressure overestimation) or volume flow (in case of power underestimation) is within the segment boundaries. A small value of $\varepsilon =10^{-3}$ is necessary to ensure that at least one plane is active at the boundaries of the piecewise defined planes. ${\overline{M}}_1$ and ${\overline{M}}_2$ are upper bounds which are necessary for the bigM-constraints. Constraints (9c) to (9h) are used to check the mentioned boundaries with binary auxiliary variables $x^{\gamma 1}$, $x^{\gamma 2}$, $x^{\delta 1}$ and $x^{\delta 2}$.

The parameters of the planes are determined using a MILP. The objective is to find planes within a segment which describe the function as good as possible and over- or underestimate it. To model this, a two-dimensional fine grid ${\mathcal {Q}}^\mathrm {grid}$ and ${\mathcal {N}}^\mathrm {grid}$ is defined: at each grid point it is constrained that the overestimation or underestimation is maintained and the distance to the nonlinear function is calculated. The set ${\mathcal {Q}}^{u}$ with $u \in \lbrace 1,..., R \rbrace$ and ${\mathcal {N}}^{v}$ with $v \in \lbrace 1,..., K \rbrace$ describes the grid points within a segment. The following applies:

$$\begin{aligned} {\mathcal {Q}}^{u} = \lbrace z \in {\mathcal {Q}}^\mathrm {grid} | Q^\mathrm {seg}_u \le z \le Q^\mathrm {seg}_{u+1}\rbrace&\quad \forall u \in \lbrace 1,..., R \rbrace , \end{aligned}$$

(10a)

$$\begin{aligned} {\mathcal {N}}^{v} = \lbrace z \in {\mathcal {N}}^\mathrm {grid} | N^\mathrm {seg}_v \le z \le N^\mathrm {seg}_{v+1}\rbrace&\quad \forall v \in \lbrace 1,..., K \rbrace . \end{aligned}$$

(10b)

With these definitions, both, MILP (11) for the pressure characteristic and MILP (12) for the power characteristics are formulated for each pump $i\in {\mathcal {P}}$ and for each segment (u, v) with $u \in \lbrace 1,..., R \rbrace , v \in \lbrace 1,..., K \rbrace$.

The objective (11a)/(12a) is to minimize the distance between the nonlinear function and its over-/underestimation, respectively. Constraint (11b) models the overestimation of the pressure characteristic and Constr. (12b) the underestimation of the power characteristic at every grid point. Since the number of grid points is limited, distance factors $\varepsilon ^{\varDelta p} = 0.005\,bar$ and $\varepsilon ^{Po} = 0.5\,W$ are introduced. These define the minimum distance between the linear planes and the nonlinear function at the grid points. Thus, it is ensured that the planes are over/underestimating also between grid points, what has been checked in the postprocessing. The result for the pressure characteristic is shown in Fig. 6 and for the power characteristic in Fig. 7.

5 Solving algorithm

The MINLP or piecewise linearized MILP presented in Sect. 4 can be solved with conventional solvers. This, however, leads to high computation times or large duality gaps for practically relevant problem sizes. For this reason, we developed a problem specific, exact solving algorithm, which allows significantly shorter computation times and leads to an acceptable result accuracy, even for the considered nonconvex MINLP.

The basic idea is to divide the problem into a number of smaller problems and to explore this solution space systematically. For this purpose, the set of all pumps ${\mathcal {P}}$ is reduced to a significantly smaller subset ${\mathcal {P}}^\mathrm {reduced} \subseteq {\mathcal {P}}$ and the purchase decisions for these pumps are fixed to one. This results in a fixation of the most important first stage decision variables. The number of pumps in the reduced set $|{\mathcal {P}}^\mathrm {reduced} |$ is less than or equal to the maximum number of pumps that can be purchased $N^\mathrm {buy}$, which drastically reduces the number of binary variables: $|{\mathcal {P}}^\mathrm {reduced}| \le N^\mathrm {buy} \ll |{\mathcal {P}}|$. If M is the number of different pump types in ${\mathcal {P}}$, then we can calculate the total number of tuples $I = |{\mathcal {T}}|$ with:

$$\begin{aligned} I = \left( {\begin{array}{c}M+N^\mathrm {buy}\\ M+1\end{array}}\right) =\frac{(M+N^\mathrm {buy})!}{(M+1)! \cdot (N^\mathrm {buy}-1)!}\,. \end{aligned}$$

(13)

For the investigated instances ($M=15$ and $N^\mathrm {buy}=4$) this results in $I=3876$ reduced MINLP with $|{\mathcal {P}}^\mathrm {reduced}|\le 4$. When considering a decentral topology (e.g. test instance B1.D, cf. 6.1) this yields a MINLP with 425 variables (64 binary) and 1057 constraints (41 nonlinear). In contrast to this is a single MINLP with $|{\mathcal {P}}| = M \times N^\mathrm {buy} = 60$ pumps. This yields 26581 variables (4380 binary) and 82562 constraints (601 nonlinear). The resulting reduced MINLP is therefore less challenging in terms of combinatorics and can generally be solved faster. The goal is to efficiently explore the set of possible purchase decisions ${\mathcal {T}}$, given the reduced set of available pumps, while also cutting off non-optimal solutions.

The solving algorithm is shown in Algorithm 2. For better clarity, we define the index set  ${\mathcal {I}}~=~\{1,2,\dots , I\}$ based on Eq. (13). First, the set of pump-tuples ${\mathcal {T}}$ is created and physically inadmissible tuples are removed during preprocessing. In addition, a set of lower bounds  ${\mathcal {L}} = \{L_1, L_2,\dots , L_{I}\}$ for the energy costs of each tuple $T_j \in {\mathcal {T}}$ is analytically calculated based on domain-specific knowledge. Besides, a set of all primal upper bounds ${\mathcal {U}} = \{U_1, U_2,\dots , U_{I}\}$ is defined and each part is set to $\infty$ before optimization. For the considered overall minimization problem, the primal bound is the minimal primal solution of a tuple $U^\mathrm {global} = \min ( {\mathcal {U}})$ (best found solution), which is infinite at the start. The dual bound of the overall minimization problem $L^\mathrm {global}$ should be raised. Therefore, the weakest dual bound in one of the subproblems/tuples is determined. Thus, the dual bound of the overall problem is bounded by the lowest dual bound of a tuple $L^\mathrm {global} \ge \min ({\mathcal {L}})$. It is initialized using the preprocessed bounds for the energy and the known investment costs of the pumps. Now the search for better primal and dual solutions is repeated until the overall problem has been solved to proven global optimality. The global primal and dual bounds are updated in each step and proven weak solutions or infeasible solutions are cut off. In the following, the single steps of the solving algorithm are explained in detail.

5.1 Preprocessing

First, the set of purchase decisions $T_j \in {\mathcal {T}}\quad \forall j \in {\mathcal {I}}$ is generated, each consisting of a reduced set of pumps ${\mathcal {P}}^\mathrm {reduced}_j \subseteq {\mathcal {P}}$ with $j \in {\mathcal {I}}$. For each tuple $T_j \in {\mathcal {T}}$, it is determined whether the pumps of the tuple ${\mathcal {P}}^\mathrm {reduced}_j$ can fulfill the required load demand. There are three elimination criteria:

• (i) The sum of the maximum pressure head of all pumps is smaller than required. Even a series connection of all pumps could not provide the required pressure.

• (ii) The sum of the maximum volume flow of all pumps is smaller than required. Even a parallel connection of all pumps could not provide the required volume flow.

• (iii) The sum of the maximum hydraulic power of all pumps is smaller than required. Even if all pumps would be operated at the point of maximum hydraulic power, it would not be sufficient.

Tuples that comply with any of these criteria are deleted immediately. This can efficiently reduce the number of tuples. Nevertheless, this excludes only a subset of all infeasible solutions.

The minimum cost are underestimated for each of the remaining tuples. Since the purchase decision for a certain pump set is fixed, the investment costs of this pump set can be calculated directly. The purchase decision for frequency converters is set to zero. We underestimate the power consumption and thus the energy cost. For each $T_j \in {\mathcal {T}}$ the underestimated power consumption $Po^\mathrm {lb, phys}_{s, j}$ is calculated by the hydraulic power $Po^{\mathrm {hydr}}_{s}$ and best efficiency of any pump in the tuple $\eta ^\mathrm {pump}_{\mathrm{best},j}$ in every scenario $s\in {\mathcal {S}}$. For each tuple $j \in {\mathcal {I}}$ the following underestimation is valid:

$$\begin{aligned} \sum _{s\in S}{ W_s \sum _{i \in {\mathcal {P}}}{ po_{i,s}}} \ge \sum _{s\in S}{W_s { Po^\mathrm {lb, phys}_{s, j}}} = \sum _{s\in S}{W_s { Po^{\mathrm {hydr}}_{s}} / { \eta ^\mathrm {pump}_\mathrm {best,j}} }. \end{aligned}$$

(14)

This corresponds to the hypothetical case that the pump with the highest efficiency provides the required hydraulic power in its best-efficiency point all by itself. The resulting cost according to Obj. (5a) underestimates the true cost and is thus a lower bound, which is based on physical principles. This consideration allows to cut off many tuples early. The power dissipated in the frequency converter is usually one order of magnitude smaller than the electrical power in the motor and is therefore neglected in this estimation.

Finally, the global primal and dual bounds are initialized in preprocessing. The global primal bound is the value of the best valid solution and is initially set to infinity. The global dual bound equals the minimum value of the lower bound of any tuple. It is initialized to the minimum of the values estimated by Eq. (14).

5.2 Bound improvement

After preprocessing, the gap between the primal and dual bound is attempted to be reduced over time. On the one hand, better primal solutions are searched for. On the other hand, the minimum values of the lower bound of the tuples are improved to increase the dual bound of the overall minimization problem. If the resulting value of the lower bound of a tuple is above the global primal bound, it can be cut off.

5.2.1 Upper bound search

In order to find a new, better primal solution, the objective value of each tuple $T_j \in {\mathcal {T}}$ is estimated in each step. The exact MINLP is then solved for the tuple with the best estimated objective value, i.e. the most promising set of pumps.

Estimation of the Objective Value The value of the primal solution of a tuple $T_j \in {\mathcal {T}}$ is estimated using a heuristic. It is assumed that the difference between the lower bound of the power consumption estimated during preprocessing $Po^\mathrm {lb, phys}_{s, j}$, Eq. (14), and the power consumptions’ exact value of the primal solution of a tuple $Po^\mathrm {primal}_{s,j}$ is comparable for all tuples $j \in {\mathcal {I}}$. We calculate this difference for all tuples in ${\mathcal {T}}$ for which we already know the primal solution by:

$$\begin{aligned} v_{s,j} = (Po^\mathrm {primal}_{s,j} -Po^\mathrm {lb, phys}_{s, j}) / Po^\mathrm {lb, phys}_{s, j} \quad \forall s \in {\mathcal {S}}, j \in {\mathcal {I}}. \end{aligned}$$

(15)

The expected power consumption $Po^\mathrm {exp}_\mathrm {s, j}$ of tuples for which we have not yet computed a primal solution is estimated using the mean of the deviation $v_{s,j}$ of all tuples with known primal solution:

$$\begin{aligned} \begin{aligned} \sum _{s\in S}{ W_s \sum _{i \in {\mathcal {P}}}{ po_{i,s}} }&\approx \sum _{s\in S}{W_s { Po^\mathrm {exp}_{s, j}}} \\ = \sum _{s\in S}&{ W_s \left( 1+ \frac{\sum _{m \in {\mathcal {I}}^\mathrm {primal} } v_{s,m} }{|{\mathcal {I}}^\mathrm {primal}|} \right) Po^\mathrm {lb, phys}_{s, j}}\quad \forall j \in {\mathcal {I}}\, , \end{aligned} \end{aligned}$$

(16)

where ${\mathcal {I}}^\mathrm {primal} = \lbrace j \in {\mathcal {I}}\,|\,\exists \, U_j < \infty \rbrace$ is the index set of tuples with finite primal solution. Based on this, the total cost is estimated according to Obj. (5a). Here, the exact known investment costs of the pumps are considered while the costs for the frequency converters are neglected.

Search for Better Primal Solution When searching for a new primal solution, the tuple with the lowest expected cost is selected. The set of pumps is reduced to the purchased pumps ${\mathcal {P}}^\mathrm {reduced}_j$ in the selected tuple $T_j$. The MINLP from Sect. 4 is formulated with the reduced set of pumps and the purchase decision variables are set to $y_i = 1~\forall \, i \in \lbrace 1,...,|{\mathcal {P}}^\mathrm {reduced}_j| \rbrace$.

In general, the resulting reduced problems for the tuples in the set ${\mathcal {T}}$ can be solved considerably faster than the overall problem. Nevertheless, high computation times may occur for some tuples $T_j\in {\mathcal {T}}$. Depending on the count for selection of the tuple, the computation time is limited. Starting with 15 s for first time computations, the limit is increased gradually for each reselection of a tuple. By this method, we try to avoid wasting computation time on suboptimal solutions. When solving the selected tuple, for feasible tuples, the primal bound $U_j$ and dual bound $L_j$ are stored for the tuple. In turn, when the subproblem of the tuple $T_j \in {\mathcal {T}}$ is infeasible, it is removed from the tuple list ${\mathcal {T}}$. We expect that the converged solution can already be cut off at an early stage without the need to solve the reduced problem to global optimality. In the best case, only a single reduced MINLP needs to be solved to global optimality. If a subproblem for a tuple $T_j \in {\mathcal {T}}$ is solved to proven global optimality, it will be deactivated for a search for further primal solutions. Unfortunately, it is not possible to continue or warm-start the computation with the selected programming framework¹ later on. Instead, the computation has to be re-started, this time with a longer maximum computation time ($\ge 15\,\hbox {s}$). For this reason, it is important to carefully select the maximum computation time for the individual steps.

If there is no promising tuple left, which we define as the lowest expected cost exceeding $120 \%$ of the global primal solution, the search for better primal solutions is no longer executed. However, the lower bound search (cf. Sect. 5.2.2) also still allows to find better primal solutions, if the reduced MINLP is solved in a later step.

5.2.2 Lower bound search

Alternating with the search for a better (i.e. lower) primal bound, the dual bound is attempted to be raised. The procedure is similar to the primal bound search: a suitable tuple is selected and a reduced optimization problem is solved. The search for a better (i.e. higher) dual bound is repeated three times before again continuing the search for a lower primal solution, as we found that it is generally more difficult than finding a better primal solution. In the dual bound search the tuple with the lowest individual dual bound $min ({\mathcal {L}})$ is selected (called weakest tuple). This is the weakest individual dual bound and determines the dual bound of the overall problem $L^\mathrm {global} \le min({\mathcal {L}})$. To increase the overall dual bound, the dual bound of the weakest tuple has to be increased by solving the corresponding subproblem. The dual bounds for tuples that have not yet been examined are the underestimated cost of the preprocessing (cf. Eq. (14)). In turn, for tuples that have already been examined, the dual bound of the tuple from the previous step is applied.

The solved optimization problem corresponds to the one from the primal bound search with a special feature: If the tuple is selected for the first time, the piecewise linear relaxed problem (PWR MILP, cf. Sect. 4.4) of the reduced MINLP is set up and solved (computation time limit of $60\,\hbox {s}$) instead of directly solving the reduced MINLP. The piecewise linear relaxed MILP is in general easier to solve and the so computed lower bound is often sufficiently high to cut off the tuple. If the solution could not be cut off, the tuple is re-selected at a subsequent loop iteration and then the reduced MINLP is solved, whereby the initial computation time is again limited to $15\,\hbox {s}$ and subsequently increased in later iterations.

5.2.3 Parallelization

To further reduce computation times, the presented solving algorithm is parallelized. In each step N optimization problems are solved simultaneously on different cores: Therefore, not only a single tuple with the best expected primal bound (in case of primal search) or lowest dual bound (in case of dual search) is selected, but N tuples with the best expected primal or lowest dual bounds are solved. If a tuple is solved before finishing the computation of any of the other first N tuples, the next best/weakest tuple is solved on the available core. The results of the computed tuples are synchronized after the computation of the first N tuples is finished.

6 Computational results

In this section the results of the optimization problems are discussed. First, the set of test instances, which are scenarios of high technological importance, are presented and the implementation and settings are clarified. Based on this, the performance of the different solving approaches is presented. Afterwards, a technical discussion of the results and the significance for practice are given.

6.1 Selection of test instances

The previously presented different systems (Sect. 3) and the wide range of variations are analyzed carefully to choose technologically meaningful test instances. Tested instances are numerated in Table 4, in which the different possible physical variants are listed. For each of the buildings B1–B3 and the test rig three different topologies are considered. (i) Parallel connection of the pumps for a central booster station ($|{\mathcal {F}}| =1$, Constr. (4) active), hereinafter referred to as parallel. (ii) Unrestricted topology for a central booster station ($|{\mathcal {F}}| =1$, Constr. (4) not active), hereinafter referred to as central. (iii) Unrestricted topology for a decentral booster station ($|{\mathcal {F}}| =5$, Constr. (4) not active), hereinafter referred to as decentral. The complexity increases significantly in the order of the listed topologies due to their increasing degrees of freedom. The heuristic presented in Sect. 3.5 is used to select 15 different pump types for each of the buildings B1–B3. Each pump type is available $N^\mathrm {avail} = 4$ times which leads to a total set of $|{\mathcal {P}}|=15\times 4 = 60$ pumps. In total, up to $N^\mathrm {buy} = 4$ pumps can be purchased and installed per booster station. For the test rig the set of pumps is limited to the $|{\mathcal {P}}|=13$ available pumps (6 different types), where up to $N^\mathrm {buy} = 6$ can be purchased and installed. While for the buildings B1–B3 a frequency converter has to be installed additionally to allow for speed control, all considered pumps for the test rig are already equipped with frequency converters and therefore the decision to buy a frequency converter is disregarded. The characteristic curves as well as costs of the pumps and frequency converters are obtained from manufacturer data. The loss factor $\zeta ^\mathrm {inst}$ is determined experimentally for the test rig, whereas it is set to zero for all real buildings. The load profiles match the load profiles presented in Sect. 3.

Table 4 Instances for testing the different solution approaches

The lifetime $T=15\,\text {a}$ for pumps is derived from the depreciation table (“AfA Table”) for the sector “Energy and Water Supply” of the German Federal Ministry of Finance (AfA table 1995). For the interest rate $r=2.08\,\%$ the discount rate (Bundesbank 2019) is used with a remaining lifetime of 15 years (7 year average, key date Aug 31st 2019). For the energy cost ($C_\mathrm {elect}=0.3\,\,$ € /kWh) the mean value of the household electricity price for the second half of 2018 in Germany is consulted (Destatis 2019). This results in an annuity-present value factor of $APF = \sum _{t=1a}^{T} 1 / (1+r)^t = 12.77$. For the investment costs of the pumps and frequency converters, we rely on market prices for the European Union of the year 2018. The pressure is calculated in bar and the volume flow in ${\hbox {m}^{3}}\hbox {h}^{-1}$. For the buildings B1–B3 the cost are calculated in k€and the power in kW, for the test rig in € and W. Thus, the coefficients of the problem become better scaled.

The three different optimization models, namely the nonlinear (MINLP), piecewise linearized (PWL), and piecewise relaxed model (PWR), can be solved using different off-the-shelf solvers and the developed solving algorithm. In this work the Solvers GUROBI v8.1.1², SCIP v5.0³ and BARON v19.7.13⁴ as well as the solving algorithm (Algorithm 2 , cf. Sect. 5.1) are trialed to compare computation time and solution quality. If nothing else is specified, default settings are used. For the PWL, 32 simplices are used for buildings B1-B3 and 50 simplices for the test rig, which are evenly distributed. We use a 1-4 orientation for each grid, as introduced in Misener and Floudas (2010). For the PWR a discretization of $R=6$ elements in q direction and $K=3$ elements in n direction for the pressure and $R=5$, $K=5$ for the power is used for real buildings. For the test rig a discretization of $R=6$, $K=2$ for the pressure and $R=4$, $K=4$ for the power is used. For reasons of simplicity, an equidistant distribution has been chosen. The PWR is primarily used within the solving algorithm trying to cut off solutions fast. For this reason, a comparatively small number of grid points has been chosen, which has proven to be efficient in test runs.

The optimization model was implemented in Matlab 2019a using YALMIP (Löfberg 2004). The solver SCIP is interfaced via the OPTI toolbox (Currie et al. 2012). For all other solvers, manufacturer interfaces are used. The solving algorithm presented in Sect. 5.1 has also been developed and executed in Matlab 2019a. All computations were performed on a Windows 10 machine with an Intel i7-8700 CPU (6 cores @3.20GHz) and 16GB RAM. $\varepsilon$-optimality (with $\varepsilon =0.1$) is tried to be achieved within a time limit of 3 hours.

6.2 Performance of different solution approaches

We first discuss the influence of the grid-size of the PWL-MILP on the computation time and the deviation from the MINLP solution, cf. Fig. 8. Afterwards, we compare the different solution approaches and solvers, for which we have computed the instances given in Table 4. A large part of the instances have a gap $> \varepsilon$, which is why a comparison of the gap at the end of the computation time is performed. As in Pecci et al. (2019), the share of solved instances with a gap smaller $\delta$, $\rho (\delta ) = (\#\,\text {instances with gap} \le \delta ) / (\#\, \text {test instances})$, is plotted for the remaining gap $\delta$, cf. Fig. 9.

Fig. 8

Influence of grid size on deviation to MINLP solution and computation time. The 2 % limit is highlighted in gray. a B1.P, b B2.P, c B3.P, d TRE.P

Comparison of MINLP and MILP The influence of the discretization is investigated using the comparatively simple instances with only parallel pump configurations (B1.P, B2.P, B3.P, TRE.P). Figure 8 shows the deviation from the best known MINLP solution and the normalized computation time over the number of cells. The deviation converges towards 0, whereby the computation time increases as a trend. Partly an overestimation, partly an underestimation of the objective value of the MINLP occurs. On the basis of the investigation, 32 cells for real buildings and 50 for the test rig are selected. With an average deviation of less than 2 % (corresponding with the order of magnitude of the measurement errors on the test rig), the level of detail is sufficient.

The comparison of the directly solved MINLP with the directly solved PWL-MILP and the PWR-MILP (cf. Fig.  9a) shows that more instances can be solved with a satisfactory gap by using one of the MILPs.

Fig. 9

Solved instances $\rho (\delta )$ with a gap smaller $\delta$ for different solution approaches. Comparison of a directly solved MINLP, PWL and PWR, b different options of solving algorithm, c different solvers for the solving algorithm (settings: w/ PWR, w/par.), and d directly solved MINLP, directly solved PWL and MINLP solved applying the solving algorithm (settings: w/ PWR, w/par.)

The PWL-MILP is an approximation of the MINLP. The relative deviation between the primal solution of the PWL-MILP and the best known primal MINLP solution (calculated for the MINLP with the solving algorithm) is on average 3.6 % for all 12 instances. For solutions with a gap of $\le \varepsilon$ (4/12 instances) the deviation is 1.65 %. For 2/12 instances, the same pumps and FCs are purchased as for the best known MINLP solution. If the PWL-MILP is solved using the solving algorithm, lower gaps can be achieved (see Fig 9c). In this case, the mean deviation is 1.61 % (all instances) or 1.75 % (instances with gap $\le \varepsilon$) and the investment decisions are the same for 5/12 instances. In comparison to the direct solution of the MINLP, more useful solutions can be achieved and the approximation quality is acceptable with around 2 % for the given discretization. However, high computation times and large duality gaps still result, especially when using an off-the-shelf solver.

The dual solution of the PWR-MILP is compared with the best found primal MINLP solution: For solutions with a gap $\le \varepsilon$ (6/12 instances) the deviation is 7.0 % or 19.91 % (all instances) and investment decisions are the same for 5/12 of the instances. The high deviation is due to non-converged solutions for difficult instances. In this case, the dual bound of the PWR-MILP cannot be raised sufficiently, which is why it serves only as a weak lower bound.

Figure 9 further shows, that the difference between SCIP and BARON to directly solve the MINLP is not significant. Both solvers can only solve 1/12 of the instances with a gap $\le \varepsilon$. In addition, we occasionally encountered errors with both solvers during the solution process. For SCIP this can be reduced by a less aggressive presolving. A more detailed discussion can be found in Appendix A. A direct solution of the MINLP with off-the-shelf solvers is due to the large duality gap and potential wrong results an unsuitable approach. In the solving algorithm the individual problems are smaller and errors by SCIP or BARON were rare. However, for SCIP the presolving-emphasis settings for the calculation with the solving algorithm were set to “fast” to avoid any problems.

Performance of Solving Algorithm Fig. 9b shows the performance of the solving algorithm with and without the use of PWR (w/ PWR and w/o PWR) and parallelization (w/ par. and w/o par.). Without the use of PWR, the reduced MINLP is already solved within the first selection of a tuple, thus the PWR-MILP is not solved before. Without parallelization, only one tuple is solved instead of N tuples at the same time. SCIP is used as MINLP solver and GUROBI as PWR-MILP solver. All in all, the deviations of the different settings are rather small. The results show:

• The use of PWR shows only minor improvements. This could be caused by a too fine discretization, resulting in a dual bound of the PWR that can only be computed with high effort.

• The parallelization offers a slight performance increase. The reason it only yields small improvements could be that when solving N tuple at the same time, the computation of one tuple takes much longer than for the others. As long as the first N tuples are not finished, the primal and dual bounds are not updated. This holds the risk of wasting time on the computation of tuples which are cut off later anyway. In addition, parallelization causes computational overhead.

In Fig. 9c, BARON and SCIP are compared as underlying MINLP solvers for the solving algorithm (w/ PWR, w/ par.). The performance of SCIP is slightly better. In addition, the PWL-MILP is solved with the solving algorithm (GUROBI as PWL-MILP solver, w/o PWR, w/ par.). The performance is slightly below that of the MINLP. A possible cause is the computational overhead of the solving algorithm: The optimization model has to be set up for each tuple. This is more complex for the PWL model due to the increased number of equations, which can lead to a loss of performance. Nevertheless, compared to directly solving the PWL-MILP, considerably more instances can be solved (cf. Fig.  9a).

Concluding, the results for the solving algorithm (SCIP, w/ PWR, w/ par.) are compared with directly solved MINLP and PWL-MILP in Fig. 9d: the solving algorithm allows to solve 50% of the instances with a gap $\le \varepsilon$ and all instances with a gap $\le$ 48%. This is substantially more than directly solving the MINLP and also than directly solving the approximated PWL-MILP.

The progression of the dual and primal bound for the solving algorithm (SCIP, w/ PWR, w/ par.) and the directly solved MINLP (SCIP, presolving “fast”) is shown in Fig. 10. The number of possible variants in the problem increases from (a) (B2.P) to (b) (B2.C) to (c) (B2.D). The solving algorithm finds the optimal primal solution very fast and at all time steps the primal bound as well as the dual bound of the solving algorithm is better than for the directly solved MINLP.

Fig. 10

Comparison of the progression of the primal and dual bound for the solving algorithm (SCIP, w/ PWR, w/ par.) and the directly solved MINLP (SCIP, presolving “fast”) for the instances a B2.P, b B2.C and c B2.D

6.3 Technical discussion of the optimized booster station designs

In the following, the technical aspects of the solutions are discussed. In this section, we always refer to the best found primal solution of any solver.

Table 5 summarizes the relevant aspects of the optimization results for all instances. The following conclusions can be drawn:

• For central booster stations the restriction of the topology to solely parallel connections has no disadvantage for the tested real buildings. This may be due to the selected pump kit, since the existing pumps are designed for parallel connections. For the test rig, the central solution (unrestricted topology) is slightly superior.

• 20–25 % of the total cost can be saved with decentralized systems. The system performance increases by 50 % up to 100 %. The latter corresponds to halving the required power. This requires higher investment costs and more pumps, which is far outweighed by the energy savings.

• Although the load cases could already be covered by a single pump, always 2–4 pumps are bought due to the achieved energy savings. In some cases different pump types are used, in others (B3.P, B3.C) solely one pump type is used. Only about half of the pumps are equipped with FCs.

For building B3, the actual installed pumps are known. In the installed system four parallelly connected pumps of the same type (all with FC) are used. To be able to compare the existing system to the optimized system, the investment decision of the existing system is fixed and the rotational speed optimized (see column “inst.” in Table 5). In comparison, approx. 50% of the total cost can be saved by optimizing the investment costs and speed at the same time. Even if costs for an additional redundant spare pump are taken into account (2053 €, without FC), 32 % of the total cost can still be saved.

Table 5 Results of the best found primal solution for each instance

Fig. 11

Selected pumps after optimization for each building and the test rig. a building B1 b building B2 c building B3 d test rig

All selected pumps are highlighted in Fig. 11. Solid black lines indicate pumps, that are selected for any setting, while dashed black lines indicate pumps that are only selected for decentralized settings.

For building B1, the optimal topologies are shown in Fig. 12. The central topology (a) shows that only the largest pump is equipped with a frequency converter and operates for every load scenarios. For the load cases with higher power demand, further pumps are switched on additionally. An interpretation of this chosen topology can be that the additional pumps run for a small part of the operation time only and therefore their energy consumption and resulting costs do not weigh enough to equip either of these pumps with a frequency converter. A complex topology results for the decentralized solution (b). In scenarios with low hydraulic power demand, one pump solely supplies the pressure zones 1 and 2. Two pumps connected in series supply the upper pressure zones. After the first pump a partial volume flow is branched off to supply pressure zone 3. Only for load cases 4 and 5 a further pump is switched on, which has no frequency converter.

Fig. 12

Optimal topology for the instances a B1.P/C and b B1.D. The labels XS to XXL refer to the hydraulic power of the pumps at their best efficiency point

The trade-off between investment costs and energy costs for building B1 is shown by the Pareto front in Fig. 13. The computation was done using the $\epsilon$-constraint method (Ehrgott 2005, pp. 98–101), in which the maximum power consumption was added as a constraint. The investment costs were weighted very high by setting the lifetime to 1 month. The optimization was performed with different upper limits on the power consumption. The Pareto front of the parallel and central topologies coincide exactly. As expected, the decentralized solution completely dominates either central solution. The minimum necessary power, calculated by hydraulic power and component efficiency, can almost be reached. This indicates a good choice of the pump kit and thus confirms the heuristics in Algorithm 1. The optimal solution discussed before (lifetime of 15 years), is in the range of high energy weighting. It should be noted that the computation time for high weights for the energy increases significantly, since fewer solutions can be cut off.

Fig. 13

Pareto Front of the objective’s investment costs and power consumption for the instances B1.P, B1.C and B1.D. The filled markers represent the best solution for a lifetime of 15 years (corresponds to Table 5)

7 Experimental validation

In order to validate the methodology, assumptions and models, the results obtained for the test rig (Sect. 3.3) are experimentally examined. The boundary conditions applied for the optimization are obtained in two different ways. In case (i), the loss factor $\zeta _\mathrm {loss}$ and the required pressure increase in the load scenarios are determined experimentally. For case (ii), a hydraulic model is developed, for which the friction coefficients of the components are estimated on the basis of data given in literature (GF Piping System 2019; Idelchik 2007). The latter case is more practical since experimental surveys cannot always be carried out in the field. In both cases, manufacturer specifications are used for the characteristic curves of the pumps and a strictly quadratic relationship between pressure losses and the volume flow is assumed. It should be noted that not all friction influences can be taken into account in the model. Depending on the topology, elbow and tee fittings are necessary to achieve the desired configuration, however, in the optimization model their friction influence is neglected. For a more detailed investigation of the friction influence, we refer to Müller et al. (2019a).

The optimal topologies are shown in Fig. 14. Only for the central booster station (b)/(c) there are deviations in the optimal topology between experimental and literature data. Without the use of experimental data, the pressure losses are assumed to be lower, which is why slower rotational speeds are applied.

Fig. 14

Optimal topology for the test rig: a parallel (experimental data as well as solely literature data), b central (experimental data), c central (solely literature data) and d decentral (experimental data as well as solely literature data). For topology c the setup on the test rig is shown in Fig. 3

The shown topologies are set up on the test rig, trying to keep the pressure losses due to the layout of the connections as low as possible. For every scenario the rotational speeds are applied and the valves are set so that the same volume flow is achieved on each floor. The pumps run at a constant speed, i.e. without pressure or volume flow control-loop.

In Fig. 15, the power consumption (measure for the achievement of the objective) is plotted against the volume flow (measure for the satisfaction of the load demand) for all topologies and scenarios.

Fig. 15

Power consumption over volume flow for optimization and validation results for all five load scenarios. Results based on experimental data for a parallel b central and c decentral topologies, and based on literature data for d parallel e central and f decentral topologies are shown. The error bars are covered by the markers

The correlation between power and volume flow is well matched in computation and experiment even though there is a deviation between the required and the measured volume flows. As expected, the correlation is higher when using the experimental data.

Furthermore, the mean volume flow,

$$\begin{aligned} Q_\mathrm {mean} := \sum _{s \in {\mathcal {S}}} W_s \sum _{f \in {\mathcal {F}}} q_{f,s}, \end{aligned}$$

(17)

the mean power consumption,

$$\begin{aligned} Po_\mathrm {mean} := \sum _{s \in {\mathcal {S}}} W_s \sum _{i \in {\mathcal {P}}} po_{i,s}, \end{aligned}$$

(18)

and the geodetic system efficiency are investigated cf. Table 6,

$$\begin{aligned} \eta _\mathrm {geo} := \frac{\sum _{s \in {\mathcal {S}}} W_s \varrho g \sum _{f \in {\mathcal {F}}} \varDelta h_f q_{f,s}}{Po_\mathrm {mean} }. \end{aligned}$$

(19)

Table 6 Comparison of optimization and validation results

The values are calculated using the density $\varrho = 998\,\hbox {kg}\,\hbox {m}^{-3}$, the gravitational acceleration $g=9.81\,\hbox {m}\,\hbox {s}^{-2}$, and the geodetic height difference $\varDelta h_f = \lbrace {1\,\hbox {m}},\dots ,5\,\hbox {m} \rbrace$ for the different pressure zones. The geodetic system efficiency $\eta _\mathrm {geo}$ is a measure for the deviation of the energy costs and is used to determine the expected total cost:

$$\begin{aligned} C^\mathrm {exp}_\mathrm {total} := C^\mathrm {opt}_\mathrm {invest} + C^\mathrm {opt}_\mathrm {energy} \left( 1+ \frac{\eta _\mathrm {geo}^\mathrm {opt} - \eta _\mathrm {geo}^\mathrm {exp}}{\eta _\mathrm {geo}^\mathrm {opt}} \right) . \end{aligned}$$

(20)

The results show that the mean volume flow is higher than demanded, i.e. the pressure losses have been overestimated. As a result, the power consumption is also higher than predicted. In practice, the deviations would be compensated by a closed loop control of each pump’s rotational speed. The systems are nevertheless comparable on the basis of the geodetic efficiency, which is slightly lower than assumed in the optimization. The deviation is smaller for the pump model based on the experimental data. For the instances where the pump S was installed (cf. Fig. 14b and c), the deviation is higher, which indicates inaccuracies in the used characteristics of pump S (given by the manufacturers data sheet). In spite of these uncertainties, the total costs are lower for decentralized systems. The deviations in efficiency (5-10 %) are significantly smaller than the saving potential compared to conventional parallel design ($\approx 50\,\%$, cf. Table 5). In addition, deviations are also to be expected for the conventional design. This shows that mathematical optimization is suitable for the planning of pumping system’s design and operation and moreover the level of detail realized in the models is sufficient.

8 Conclusion and outlook

We presented three different modeling approaches for the optimization of booster stations in high-rise buildings that are based on a stochastic optimization approach. The underlying optimization model is mostly linear in its decision variables, but the consideration of the pump characteristics for the booster station design leads to a nonconvex mixed-integer nonlinear program, which is hard to solve for real-world instances in its original form. To overcome this drawback, we derived two additional model variants, allowing for MILP formulations. These are a piecewise linear approximation and a piecewise linear relaxation of the nonlinear pump characteristics. Furthermore, we implemented a problem-specific solving algorithm which is based on these different methods and additional domain-specific engineering knowledge to derive global-optimal results for real-world instances efficiently. The solving algorithm reduces the original MINLP by fixing specific pump purchase decisions and uses a selection heuristic to find feasible primal solutions in a very short time. This primal solution search is undertaken alternating with the dual solution search, which is based either on the reduced MINLP directly or on the implemented piecewise linear relaxation. The presented solving algorithm improves the solution times considerably in comparison to the original MINLP formulation. We presented the computation benefits of the solving algorithm for three different high-rise buildings as well as for a test rig. The test rig, constructed for validation purposes, represents the scaled water supply in a high-rise building. Therefore, we were able to conduct experiments for results of the optimization models considering different data sources as well as varying degrees of freedom. The comparison of optimization and experimental results on the test rig shows that the degree of detail is adequate and that the introduced method is transferable to real-world infrastructure. Allowing for different degrees of freedom for the alignment of pumps in the water supply system shows that the decentralized topology of booster stations can yield approx. 25 % of cost savings when compared to an optimized parallel pump alignment. The latter design is the current common practice. Additionally, frequency converters were installed for only approx. half of the pumps in booster stations, which is due to the trade-off of energy costs in operation and investment costs making up total cost. Compared to a conventionally designed booster station, the mathematical optimization could achieve a total cost reduction of up to approx. 50 %. These results show that we were able to derive significantly improved system designs by combining domain-specific engineering knowledge and discrete optimization methods. The developed approach assists engineers as early as in the design stage to improve the system design. This is possible, since this approach allows to consider the real load profile in the implemented two-stage stochastic optimization problem. For the future, we want to address further open research questions in this area, both in the algorithmic and engineering domain. For instance, to obtain reasonable problem sizes, we had to reduce the number of different pump types in the preprocessing with the help of an engineering knowledge applying heuristic. Therefore, improving the performance of the algorithm to consider more pumps could achieve even better system designs. This objective could be approached by an improved deployment of the underlying piecewise relaxation or by warm-starting already investigated subproblems. From an engineering perspective, we want to examine more complex, branched systems and want to look into the transfer of the methods to different fluid systems. Finally, we would like to examine the operation in more detail and ensure that the best possible performance is achieved even with load variations and deviations of the assumed parameters.

A Appendix

Within the computational evaluations we encountered problems with the two MINLP solvers SCIP (v5.0) and BARON (v.19.7.13), which we present in more detail in the following. Within modeling, we put a lot of effort in the conditioning of the underlying optimization problem to reduce errors based on an ill-conditioned optimization problem. A summary can be found in Table 7. While BARON occasionally generated no output for specific instances, SCIP either indicated larger dual bounds than the value of the true primal solutions or declared instances as infeasible while we were able to derive feasible solutions either with BARON or manually by engineering knowledge. BARON terminated without solution output for 4/12 instances for the original MINLP given by (5) and (6). SCIP output an incorrect solution for 5/12 instances in the original MINLP with default presolving settings. We were able to reduce the number of false results for SCIP by adjusting the presolving emphasis settings to “fast” or “off”, as indicated in Table 7. A newer version of SCIP or BARON can potentially solve the encountered issues.

Table 7 Specification of the instances for which an error occurred in the direct solution of the MINLPs