1 Introduction

The optimal design of water distribution networks (WDN) is a classical problem which has traditionally consisted in selecting pipe diameters, for a given network topology, with the objective of minimizing design cost. Recently, renewed interest for the problem of optimal WDN design has emerged, as a result of the need to expand, rehabilitate or sectorize networks, to cope with the deterioration of assets and increasing customer demands (Giudicianni et al. 2020). In this context, the problem formulation no longer assumes a fixed network topology, but aims to simultaneously optimize network design and control, taking into account a wider range of objectives. In the UK, the requirements of the water services regulation authority to minimize water losses and interruptions of customer supply (Ofwat 2019) have led to the implementation of pressure control and network rehabilitation measures. This work investigates the optimal design-for-control (DfC) of existing WDNs with the objectives of minimizing pressure induced background leakage, which makes up a large part of total water losses, and maximizing network resilience, which reflects the ability of a WDN to maintain continuous customer supply. To the best of the authors’ knowledge, this problem has not been considered in previous literature.

The presented design-for-control problem consists in simultaneously selecting elements to be added to existing WDNs from predefined sets of candidate valves (CNV) and pipes (CNP), and optimizing the controls of new and existing pressure control valves. WDNs are modeled as directed graphs, subject to the conservation of mass and energy. Flows through network links as well as pressure heads at network nodes, and head differences across pipes and valves are represented by continuous variables, while the selection status of CNVs and CNPs is represented by binary variables. Friction losses across network links are modeled using the empirical Hazen-Williams or Darcy-Weisbach equations. Finally, we consider the objectives of minimizing pressure induced background leakage and maximizing network resilience, represented respectively by the average zone pressure (AZP) (Wright et al. 2015) and the resilience index (\(I_r\)) (Todini 2000). This results in a non-convex bi-objective mixed-integer non-linear program (BOMINLP).

As AZP and \(I_r\) are both increasing functions of pressure, minimizing AZP and maximizing \(I_r\) are conflicting objectives. As a result, solving the DfC problem requires identifying the best compromises between the objectives, which are characterized using the notions of efficiency and non-dominance: a solution in the decision space is said to be efficient or Pareto optimal if its image in the objective space is non-dominated, i.e. no other feasible solution can improve on one of the objectives without making the others worse (Miettinen 1998). In the case of bi-objective mixed-integer linear programs, the non-dominated set (or Pareto front) consists of a subset of the convex hull of a finite set of non-dominated points (Yu and Zeleny 1975). As a result, it is possible to obtain a complete description of the Pareto front using multi-objective branch-and-bound (Mavrotas and Diakoulaki 1998). More recently, De Santis et al. (2020) proposed a branch-and-cut method to find the finite set of non-dominated points of purely integer bi-objective programs. However, Pareto fronts of multi-objective mixed-integer non-linear programs (MOMINLP) may be composed of disconnected segments of continuous curves and isolated points (Das 2000). Hence, obtaining a complete description of the efficient and non-dominated sets (resp. Pareto set and Pareto front) can be impractical in the case of MOMINLPs. The heuristic algorithm proposed by Cacchiani and D’Ambrosio (2017), which returns a subset of non-dominated points, is the first general purpose method for the solution of convex MOMINLPs. In this paper, we are interested in methods to approximate the Pareto front of the considered non-convex bi-objective DfC problem.

Heuristic approaches, based on genetic algorithms (GAs) or simulated annealing, and mathematical optimization methods have been applied to the design and/or control of WDNs, including optimal pipe sizing (Xu and Goulter 1999; Prasad and Tanyimboh 2009; Cunha and Marques 2020), network rehabilitation (Walters et al. 1999; Alvisi and Franchini 2009) and valve placement and control (Nicolini and Zovatto 2009; Eck and Mevissen 2012; Dai and Li 2014; Pecci et al. 2017a, c) problems. GAs and simulated annealing require a significant number of hydraulic simulations and objective function evaluations. Moreover, they do not provide global optimality (in the single-objective case) or efficiency/non-dominance (in the multi-objective case) bounds to guarantee the quality of the computed solutions. Global optimization methods, on the other hand, return, for single-objective problems, a feasible solution with certified global optimality bounds (Eiger et al. 1994; Gleixner et al. 2012; D’Ambrosio et al. 2015; Pecci et al. 2018). For specific cases (e.g. convex mixed-integer problems or non-convex continuous problems), global optimization methods have been extended to generate guarantees of global efficiency/non-dominance for multi-objective problems, using scalarization approaches (Fernández and Tóth 2007; Ehrgott and Gandibleux 2007) or multi-objective branch-and-bound (De Santis et al. 2019; Niebling and Eichfelder 2019). Instead of a strict subset of the Pareto front, these methods return a set of feasible solutions and a guarantee of their non-dominance, within a given tolerance, in the form of a tight superset of the efficient and/or non-dominated sets.

In this paper, we present a method to approximate the non-dominated set of the design-for-control of existing water distribution networks problem using a superset of the Pareto front. The main contributions of this work include the formulation of the bi-objective design-for-control problem for the joint minimization of AZP and maximization of \(I_r\) in WDNs, the development of a method to approximate the Pareto front of the resulting BOMINLP and compute potentially non-dominated solutions with global non-dominance bounds, and, finally, an equivalent linear reformulation to facilitate the computation of global optimality bounds for the maximization of the resilience index \(I_r\).

In Sect. 2, we formulate the DfC problem as a non-convex BOMINLP, and Sect. 3 presents a novel tailored method to approximate its Pareto front. First, the BOMINLP is scalarized into a sequence of single-objective mixed-integer non-linear problems (MINLP) using an \(\epsilon \)-constraint method; a tailored spatial branch-and-bound (sBB) method is then applied to compute feasible solutions of the single-objective MINLPs with certified global optimality bounds. In particular, the formulation of the resilience index includes linear-fractional terms, resulting in a non-linear function. As a consequence, previous literature has solved optimization problems involving the resilience index using meta-heuristic approaches (Todini 2000; Atkinson et al. 2014; Alvisi 2015; Creaco et al. 2016). We adapt a previous sBB algorithm from the literature (Pecci et al. 2018) and propose an equivalent linear reformulation of \(I_r\), which enables the linearization of the parametrized MINLPs and the computation of lower bounds throughout the sBB search. We show that the proposed approach returns a set of potentially non-dominated solutions as well as global non-dominance bounds, in the form of a superset of the Pareto front. In Sect. 4, we apply this method to two case study networks. We show that the proposed sBB method outperforms state-of-the-art global optimization solvers SCIP and BARON, and converges to good quality feasible solutions. We comment on the computational results and hydraulic characteristics of the network configurations corresponding to the obtained set of potentially non-dominated solutions, before presenting our conclusions in Sect. 5.

2 Problem formulation

In this paper, we consider a water distribution network with \(n_n\) demand nodes or junctions, \(n_0\) water sources and \(n_p\) links, among which \(n_{\text {V}}\) existing pressure control valves, \(n_{\text {CNV}}\) candidate links for valve installation and \(n_{\text {CNP}}\) candidate pipes for installation. Pumps are not considered and we assume known the location of existing valves, CNVs and CNPs, as well as the length, diameter and material of CNPs. The problem consists in simultaneously selecting, from the predefined sets of CNVs and CNPs, the candidate valves and pipes to add to the network, and optimizing the control of existing and new valves.

We model the network as a directed graph with \(n_n + n_0\) nodes (vertices) and \(n_p\) links (edges). A positive flow through a link is defined as directed from start to end node. We refer to V as the set of network nodes (vertices) and E as the set of network links (edges). Index subsets \(E_{\text {CNV}}\) and \(E_{\text {CNP}}\) correspond to CNV and CNP locations, respectively. We consider \(n_t\) different time steps of network operation and, for \(k=1,\dots ,n_t\), the vector of known nodal demands is given by \(d^k \in \mathbb {R}^{n_n}\) and known hydraulic heads at water sources by \(h_0^k \in \mathbb {R}^{n_0}\). Let \(\zeta \) and \(h^{\min } \in \mathbb {R}^{n_n}\) be the vectors of known elevations and minimum allowed hydraulic heads at network nodes, respectively. For the case study networks considered here, we set, for \(i=1,\dots ,n_n\),

$$\begin{aligned} \begin{aligned} h^{\min }_i= {\left\{ \begin{array}{ll} \zeta _i + 10 \text {m}, \quad \text {if} \underset{k=1,\dots ,n_t}{\max }(d^k_i)>0, \\ \zeta _i \text {m}, \qquad \quad \text {otherwise}. \end{array}\right. } \end{aligned} \end{aligned}$$

At time step k, we consider the continuous decision variables \(h^k \in \mathbb {R}^{n_n}\) and \(q^k \in \mathbb {R}^{n_p}\), representing the hydraulic heads at the network demand nodes and the vector of flow rates, respectively. We use vector \(\eta ^k \in \mathbb {R}^{n_p}\) to model additional head losses introduced by the action of the valves or head differences across closed pipes. Discrete decision variables \(z \in \{0;1\}^{(n_{\text {CNV}}+n_{\text {CNP}})}\) represent the selected (1) and unselected (0) status of CNVs and CNPs.

2.1 Constraints

At time step k, hydraulic heads \(h^k\), flows \(q^k\), and additional head loss across links \(\eta ^k\) are subject to energy and mass conservation constraints:

$$\begin{aligned} \begin{aligned} A_{12}h^k + A_{10}h^k_0 + \phi (q^k) + \eta ^k = 0, \quad k=1,\dots ,n_t \end{aligned} \end{aligned}$$
(1a)
$$\begin{aligned} \begin{aligned} A_{12}^T q^k - d^k = 0, \quad k=1,\dots ,n_t \end{aligned} \end{aligned}$$
(1b)

where \(A_{12} \in \mathbb {R}^{n_p \times n_n}\) and \(A_{10} \in \mathbb {R}^{n_p \times n_0}\) are the link node incidence matrices for demand nodes and water sources respectively. Friction head losses associated with pipe flows \(q^k= [q_1^k \dots q_{n_p}^k ]^T\), represented by the vector \(\phi (q^k)=[\phi _1(q_1^k) \dots \phi _{n_p}(q_{n_p}^k) ]^T\), are commonly described either by the Darcy–Weisbach (D–W) or Hazen–Williams (H–W) equations. As neither equation is smooth, it is convenient to use a quadratic approximation of the formulae over the operational range of flows in each link. This produces a pair of positive coefficients \((a_j,b_j) \in (\mathbb {R}_+)^2\) for each pipe j in the network and the vector of friction head loss at time step k is then given by \(\phi (q^k)=[\phi _1(q^k_1) \dots \phi _{n_p}(q^k_{n_p})]^T\), where \(\phi _j(q^k_j) = q^k_j(a_j\left| q^k_j\right| +b_j)\). This is a common approach for the formulation of head loss constraints in WDN optimization literature (Dai and Li 2014; Eck and Mevissen 2015; Pecci et al. 2017b; Zamzam et al. 2019). Moreover, the results presented in “Appendix F” suggest that the use of a quadratic approximation does not affect the quality of the solutions computed for the case study networks considered in this work. The error in the objective values, in particular, does not exceed \(2\%\).

To isolate the non-convex friction head loss term, we introduce the vector of auxiliary variables \(\theta ^k \in \mathbb {R}^{n_p}\) and constraint (1a) becomes:

$$\begin{aligned} \begin{aligned} A_{12}h^k + A_{10}h_0^k + \theta ^k + \eta ^k = 0, \quad k=1,\dots ,n_t \end{aligned} \end{aligned}$$
(2a)
$$\begin{aligned} \begin{aligned} \theta ^k = \phi (q^k), \quad k=1,\dots ,n_t \end{aligned} \end{aligned}$$
(2b)

where (2b) is the non-convex friction loss equality constraint, while constraints (1b) and (2a) are linear. Constraint (3)

$$\begin{aligned} \begin{aligned} A_{10}^T q^k \le 0, \quad k=1,\dots ,n_t\ \end{aligned} \end{aligned}$$
(3)

ensures water is flowing out from the network sources. We also define vectors \(q_U^k\), \(q_L^k\), \(h_U^k\), \(h_L^k\), \(\theta _U^k\), \(\theta _L^k\), \(\eta _U^k\) and \(\eta _L^k \in \mathbb {R}^{n_p}\) of upper and lower bounds on variables \(q^k\), \(h^k\), \(\theta ^k\) and \(\eta ^k\), respectively, such that

$$\begin{aligned} \begin{aligned} q_L^k \le q^k \le q_U^k, \quad k=1,\dots ,n_t\ \\ h_L^k \le h^k \le h_U^k, \quad k=1,\dots ,n_t\ \\ \theta _L^k \le \theta ^k \le \theta _U^k, \quad k=1,\dots ,n_t\ \\ \eta _L^k \le \eta ^k \le \eta _U^k, \quad k=1,\dots ,n_t\ \\ \end{aligned} \end{aligned}$$
(4)

where \(h_U^k=\mathop {\max }\nolimits _{i=1,\ldots ,n_0}((h_0^k)_i) \times \mathbb {1}_{n_n}\), \(h_L^k=h^{\min }\) and \(q_U^k\), \(q_L^k\), \(\theta _U^k\), \(\theta _L^k\), \(\eta _U^k\) and \(\eta _L^k\) are defined for pipes, valves and CNVs, and CNPs by the expressions (A.1), (A.2) and (A.3) in “Appendix A”.

To model the status of the candidate valves and links, we introduce the binary variable \(z \in \mathbb {R}^{n_{\text {CNV}}+n_{\text {CNP}}}\): for \(j \in E_{\text {CNV}} \cup E_{\text {CNP}}\), the status of j is represented by \(z_{l_j}\), \(l_j \in \{1 \dots n_{\text {CNV}}+n_{\text {CNP}}\}\). For \(j \in E_{\text {CNV}}\), either \(z_{l_j}=1\) and the candidate link j is selected for valve installation and modeled as a pressure control valve, or \(z_{l_j}=0\) and link j is unselected for valve installation, i.e. modeled as an open pipe. For \(j \in E_{\text {CNP}}\), either \(z_{l_j}=1\) and the candidate location j is selected for pipe installation and modeled as an open pipe, or \(z_{l_j}=0\) and j is modeled as a closed pipe. Additional head losses across unselected CNVs and flows across closed pipes are forced to zero, while hydraulic heads at the start and end nodes of closed pipes are decoupled.

This is enforced by the following big-M constraints on variables \(q^k\), \(\theta ^k\) and \(\eta ^k\):

$$\begin{aligned} \begin{aligned} q^k + Q_U^k (\mathbb {1}_{(n_{\text {CNV}}+n_{\text {CNP}})}-z) \le q_U^k, \quad k=1,\dots ,n_t\ \\ -q^k - Q_L^k (\mathbb {1}_{(n_{\text {CNV}}+n_{\text {CNP}})}-z) \le -q_L^k, \quad k=1,\dots ,n_t\ \\ \end{aligned} \end{aligned}$$
(5a)
$$\begin{aligned} \begin{aligned} \theta ^k + \Theta _U^k (\mathbb {1}_{(n_{\text {CNV}}+n_{\text {CNP}})}-z) \le \theta _U^k, \quad k=1,\dots ,n_t\ \\ -\theta ^k - \Theta _L^k (\mathbb {1}_{(n_{\text {CNV}}+n_{\text {CNP}})}-z) \le -\theta _L^k, \quad k=1,\dots ,n_t\ \\ \end{aligned} \end{aligned}$$
(5b)
$$\begin{aligned} \begin{aligned} \eta ^k + N_U^k (\mathbb {1}_{(n_{\text {CNV}}+n_{\text {CNP}})}-z) \le \hat{\eta }_{U}^k, \quad k=1,\dots ,n_t\ \\ -\eta ^k - N_L^k (\mathbb {1}_{(n_{\text {CNV}}+n_{\text {CNP}})}-z) \le -\hat{\eta }_{L}^k, \quad k=1,\dots ,n_t\ \\ \end{aligned} \end{aligned}$$
(5c)

where diagonal matrices \(Q_U^k\), \(Q_L^k\), \(\Theta _U^k\), \(\Theta _L^k\), \(N_U^k\) and \(N_L^k \in \mathbb {R}^{n_p \times (n_{\text {CNV}}+n_{\text {CNP}})}\) are defined by (A.4a)–(A.4c) and vectors \(\hat{\eta }_{L}^k\) and \(\hat{\eta }_{U}^k \in \mathbb {R}^{n_p}\) by (A.5) in “Appendix A”.

Further constraints on binary variables z, to prevent specific combinations of CNVs and CNPs from being selected for instance, can be enforced with

$$\begin{aligned} \begin{aligned}&Az \le \mathbf {e} \\ \end{aligned} \end{aligned}$$
(6)

where \(A \in \mathbb {R}^{n_A \times (n_{\text {CNV}}+n_{\text {CNP}})}\) is a constant matrix of binary inequalities representing combinatorial constraints on the feasible sets of candidate valves and pipes and \(\mathbf {e} \in \mathbb {R}^{n_A}\) is a \({n_A}\) by 1 vector with fixed entries. Finally, we define the vector of decision variables \(x=(z,q,h,\eta ,\theta )\), where \(q=(q^k)_{k=1,\dots ,n_t}\), \(h=(h^k)_{k=1,\dots ,n_t}\), \(\eta =(\eta ^k)_{k=1,\dots ,n_t}\), \(\theta =(\theta ^k)_{k=1,\dots ,n_t}\), and the set \(\text {F}(Q) \subset \{0,1\}^{n_{\text {CNV}}+n_{\text {CNP}}} \times \mathbb {R}^{n_p} \times \mathbb {R}^{n_n} \times \mathbb {R}^{n_p} \times \mathbb {R}^{n_p}\) given by (1b), (2a), (3), (4), (5) and (6), which depends on the bounds \(Q {:}{=}\{ q_L^k,q_U^k \}_{k=1,\dots ,n_t}\).

2.2 Objective functions

In this work, we consider the bi-objective problem of the optimal design-for-control of water distribution networks, with the objectives of minimizing pressure induced background leakage and maximizing resilience.

Pressure induced background leakage Water utilities in the UK are incentivized, by means of financial penalties, to reduce leakage throughout their systems. Background leakage, in particular, which represents a large part of total water losses, depends on pressure levels in the distribution network (Van Zyl and Clayton 2007). Here, we use average zone pressure (Wright et al. 2015) as a surrogate measure for pressure induced leak flow:

$$\begin{aligned} \begin{aligned} \text {AZP}(x)&=\frac{1}{n_{t}W} \sum _{k=1}^{n_t} w^T (h^k-\zeta ) \\ \end{aligned} \end{aligned}$$
(7)

where \(w \in \mathbb {R}^{n_n}\) and \(W \in \mathbb {R}\) are the vector of AZP weights and the AZP normalization factor, respectively. For \(i \in V\), we define \(J_i\), the set of all pipes incident to node i. Then \(w_i=\sum _{j \in J_i} L_j/2\) and \(W=\sum _{i=1}^{n_n}w_i\), where \(L \in \mathbb {R}^{n_p}\) is the vector of pipe lengths.

Resilience Another key performance commitment concerns the reduction of interruptions of supply. The resilience of a WDN, which refers to its ability to maintain continuous customer supply under failure conditions, depends on the availability of alternative supply paths (topological redundancy) and pressure surplus (energy redundancy). The resilience of a WDN can be measured directly as its performance under failure conditions. This approach however suffers from combinatorial limitations, as it requires the simulation of a wide range of failure events. In comparison, surrogate measures of resilience depend only on the structure of the network, or its performance under normal operating conditions. A commonly adopted surrogate measures of resilience is the resilience index \(I_r\) (Todini 2000), which is defined as the ratio of surplus energy to total energy input measured under normal operating conditions:

$$\begin{aligned} \begin{aligned} I_r(x)&= \frac{\sum _{k=1}^{n_t} \left( (d^k)^T h^k - (d^k)^T h^{\min } \right) }{\sum _{k=1}^{n_t} \left( (-A_{10} h_0^k)^T q^k - (d^k)^T h^{\min } \right) }. \end{aligned} \end{aligned}$$
(8)

Both the numerator and denominator in (8) are affine functions of the vectors of hydraulic heads h and flows q respectively: in the general case, \(I_r\) is a linear fractional function of the problem variables. Moreover, for all \(x=(z,q,h,\eta ,\theta )\) satisfying (1b), (2), (3), (4), (5), (6) and \(z \in \{0;1\}^{n_{\text {CNV}}+n_{\text {CNP}}}\), the denominator of \(I_r(x)\) is strictly positive and \(0 \le I_r(x) \le 1\)—see “Appendix B”.

We observe, from (7) and (8), that both AZP and \(I_r\) are increasing functions of the hydraulic heads at the network nodes. As a result, minimizing pressure induced background leakage and maximizing WDN resilience are conflicting objectives. The optimal DfC problem for the joint minimization of AZP and maximization of \(I_r\) (i.e. minimization of \(-I_r\)) can be formulated as a bi-objective (non-convex) mixed-integer non-linear program:

figure a

where \(f_1\) and \(f_2\) represent either objective functions AZP or \(-I_r\), (2b) is the non-convex potential-flow coupling constraint and F(Q) is the set defined by the linear constraints (1b), (2a), (3), (4), (5) and (6). Solving \(\text {BOMINLP}(Q)\) requires identifying the best compromises between the objective functions \(f_1\) and \(f_2\), or non-dominated set:

Definition 1

  1. 1.

    A feasible solution \(x^*\) of \(\text {BOMINLP}(Q)\) is efficient if there is no x feasible for \(\text {BOMINLP}(Q)\) such that \((f_1(x),f_2(x)) \ne (f_1(x^*),f_2(x^*))\), \(f_1(x) \le f_1(x^*)\) and \(f_2(x) \le f_2(x^*)\).

  2. 2.

    The point \((f_1(x^*),f_2(x^*))\) in the objective space is non-dominated for \(\text {BOMINLP}(Q)\) if \(x^*\) is an efficient solution of \(\text {BOMINLP}(Q)\). The set of non-dominated solutions of \(\text {BOMINLP}(Q)\) is called the non-dominated set, or Pareto front of \(\text {BOMINLP}(Q)\).

In the following section, we present a method to approximate the Pareto front of \(\text {BOMINLP}(Q)\).

3 Solution method: \(\epsilon \)-constraints and sBB

The bi-objective design-for-control problem aims to identify trade-off solutions representing efficient compromises between conflicting objectives, as no single solution jointly optimizes all objectives. In particular, we are interested in deterministic methods for computing potentially non-dominated solutions of \(\text {BOMINLP}(Q)\) with guarantees of global non-dominance.

In Fernández and Tóth (2007); Ehrgott and Gandibleux (2007), the authors present scalarization-based approaches to compute outer approximations of the Pareto set of bi-objective non-linear problems (NLP) and the Pareto front of bi-objective mixed-integer linear problems (MILP). Previous studies (Cacchiani and D’Ambrosio 2017; De Santis et al. 2019; Niebling and Eichfelder 2019) obtain similar guarantees of global efficiency and/or non-dominance by extending the early work of Mavrotas and Diakoulaki (2005); Fernández and Tóth (2007) on multi-objective branch-and-bound to compute supersets of the Pareto sets and/or Pareto fronts of a wider range of multi-objective optimization problems, including convex MINLPs and continuous NLPs with non-convex objective functions. None of these approaches has, however, considered BOMINLPs with non-convex constraints.

This work investigates the application of scalarization, together with global optimization methods, to approximate the Pareto front of \(\text {BOMINLP}(Q)\). We implement the method of \(\epsilon \)-constraints (Ehrgott 2006) to scalarize the bi-objective problem (see Sect. 3.1) and use spatial branch-and-bound to compute upper and lower bounds on the optimal values of the resulting sequence of parametrized single-objective problems (see Sect. 3.2). We show that, using this approach, we can compute a superset of the Pareto front of \(\text {BOMINLP}(Q)\), i.e. potentially non-dominated solutions as well as guarantees of their global non-dominance (see Sect. 3.3).

3.1 Scalarization using the method of \(\epsilon \)-constraints

Scalarization approaches convert a multi-objective problem into a sequence of parametrized single-objective problems. By varying parameter values, different solutions from the efficient and/or non-dominated sets are produced. For a review of scalarization methods, see Ehrgott (2006). In this work, we focus in particular on the method of \(\epsilon \)-constraints. In \(\epsilon \)-constraints, one objective function is selected to be minimized while upper bounds are enforced on other objective function. We scalarize \(\text {BOMINLP}(Q)\) by introducing an \(\epsilon \)-constraint on the second objective function \(f_2\) and reformulate it into a sequence of \(n_{\epsilon }\) (non-convex) parametrized MINLPs:

figure b

where \(\epsilon _j \in [0,1 ]\), for all \(j=1,\dots ,n_{\epsilon }\), and the anchor point \(x_{f_i}^*\), \(i=1,2\), is an optimal solution of the single-objective non-convex problem \(\text {MINLP}_{f_i}(Q)\):

figure c

As a result, the quality of the approximation of the Pareto front depends on the order of the objective functions in \(\text {BOMINLP}(Q)\) and choice of scalarization parameters \(\epsilon _j,\ j=1,\dots ,n_{\epsilon }\). We discuss this aspect, as well as the scalability challenges associated with the extension of the proposed method to multi-objective problems with more than two objectives, in Sect. 4.3.

Feasible solutions with certified optimality bounds can be computed, for the single-objective problems \(\text {MINLP}_{f_i}(Q)\), \(i=1,2\), and \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\), \(j=1,\dots ,n_{\epsilon }\), using global optimization methods such as the tailored spatial branch-and-bound (sBB) procedure presented in Pecci et al. (2018); Ulusoy et al. (2019)—see Sect. 3.2. For a review of global optimization methods, including spatial branch-and-bound, see Belotti et al. (2013).

The use of the \(\epsilon \)-constraints scalarization method has two main benefits. First, we show (see Sect. 3.2.2) that relaxing the potential-flow coupling constraints in \(\text {MINLP}_{f_i}(Q)\) or \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\) produces a problem which can be equivalently reformulated into an mixed-integer linear program (MILP). This provides a convenient mean to compute tight lower bounds throughout the sBB procedure. In addition, we show that the sets of lower and upper bounds on the optimal value of the single-objective problems \(\text {MINLP}_{f_i}(Q)\), \(i=1,2\), and \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\), \(j=1,\dots ,n_{\epsilon }\), computed by sBB define the boundary of a superset of the Pareto front of \(\text {BOMINLP}(Q)\)—see Sect. 3.3.

3.2 Spatial branch-and-bound

In this section, we present a tailored spatial branch-and-bound (sBB) algorithm to compute good quality feasible solutions, with certified global optimality bounds, for a non-convex mixed-integer non-linear program \(\text {MINLP}(Q)\), of the form \(\text {MINLP}_{f_i}(Q)\) or \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\), defined on the initial bound set \(Q=Q^{\text {init}}\). The sBB method described in Algorithm 1 recursively divides the initial bound set \(Q^{\text {init}}\) into increasingly smaller sets Q stored in the working list \(\mathcal {Q}\), creating a tree of descendant subproblems \(\text {MINLP}(Q)\). Lower bounds on the optimal value of the subproblems are obtained by solving mixed-integer linear relaxations. In particular, we propose a tailored lower bounding procedure, which relies on an equivalent linear reformulation of the resilience index \(-I_r\). In comparison, upper bounds correspond to locally optimal solutions of the NLPs obtained by fixing the values of the integer variables in the non-convex MINLP subproblems. Global optimality is then certifiably reached when the relative optimality gap between upper bound and lower bound is smaller than a specified tolerance, \(\delta \) (here we set \(\delta =1e^{-2}\)—see Algorithm 1).

For \(i=1,2\) (resp. \(j=1,\dots ,n_{\epsilon }\)), the sBB algorithm returns a feasible solution \(x_{f_i}\) (resp. \(x_{f_1,f_2,\epsilon _j}\)), an upper bound \(\text {UB}_{f_i}\) (resp. \(\text {UB}_{f_1,f_2,\epsilon _j}\)) and a lower bound \(\text {LB}_{f_i}\) (resp. \(\text {LB}_{f_1,f_2,\epsilon _j}\)) on the optimal value of \(\text {MINLP}_{f_i}(Q^{\text {init}})\) (resp. \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q^{\text {init}})\)), and the set of nodes \(\mathcal {Q}_{f_i}\) (resp. \(\mathcal {Q}_{f_1,f_2,\epsilon _j}\)) explored during the search. Finally, to compute the right-hand side of the \(\epsilon \)-constraint in \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\), we define \(x_{f_1}^*=x_{f_1}\) and \(x_{f_2}^*=x_{f_2}\).

3.2.1 Selection strategy, branching strategy and bound tightening

After initializing the sBB search, and as long as the relative optimality gap is greater than the set tolerance \(\delta \) and the time limit \(T_{\max }\) has not been exceeded, the set of bounds \(Q^s\) corresponding to the subproblem \(\text {MINLP}(Q^s)\) with the smallest lower bound L\((Q^s)\) is removed from the working list \(\mathcal {Q}\) and divided. The resulting descendant bound sets are then reduced using an Optimality-Based Bound Tightening (OBBT) method—see Puranik and Sahinidis (2017) for a review of domain reduction techniques—and stored in the working list \(\mathcal {Q}\). Although OBBT is associated with a significant computational overhead, the implementation of domain reduction methods throughout the sBB search was shown to improve final bounds on the optimal value of non-convex problems more efficiently than branching alone (Castro and Grossmann 2014; Ulusoy et al. 2019). In this work, we implement the domain reduction procedure described in Algorithm 1 in Pecci et al. (2018), with \(K_{\max }=1\). For more details about the selection and branching strategies, and implementation of OBBT, see Pecci et al. (2018); Ulusoy et al. (2019).

3.2.2 Lower bounding

To obtain lower bounds L(Q) on the optimal value of the subproblems \(\text {MINLP}(Q)\) generated throughout the sBB search, we relax the subproblems by computing linear outer approximations of the non-convex potential-flow coupling constraint (2b) and reformulate the resilience index \(-I_r\) to produce relaxed MILPs. The solution of the resulting linear problems provides lower bounds on the optimal value of subproblems \(\text {MINLP}(Q)\).

Relaxation of the potential-flow coupling constraints To relax \(\text {MINLP}(Q)\), we compute a polyhedral envelope of the potential-flow coupling constraint (2b) implementing a method for monomials of odd degree developed by Liberti and Pantelides (2003). The original potential-flow coupling constraint (2b) is replaced by the linear constraint:

$$\begin{aligned} \begin{aligned} R^k q^k + E^k \theta ^k \le r^k, \forall k=1,\dots ,n_t, \end{aligned} \end{aligned}$$
(9)

where matrices \(R^k\) and \(E^k \in \mathbb {R}^{m \times n_p}\) and vector \(r^k \in \mathbb {R}^{m}\) depend on the bounds Q and number of linearizations used to compute the linear approximation of the convex envelope of constraint (2b). For the sake of brevity, we omit a detailed description of these relaxations, which can be found in Pecci et al. (2018); Ulusoy et al. (2019). The relaxation of constraint (2b) in \(\text {MINLP}(Q)\) results either directly in an MILP or, for \(\text {MINLP}_{f_i}(Q)\) with \(f_i(x)=-I_r(x)\), in the mixed-integer linear fractional problem:

$$\begin{aligned} \begin{aligned}&\underset{x=(z,q,h,\eta ,\theta )}{\text {minimize}}&-I_r(x) \\&\text {subject to}&R^k q^k + E^k \theta ^k \le r^k,\ \forall k \in \{1,...,n_t\},\quad (9) \\&&(q^k,h^k,\theta ^k,\eta ^k) \in \text {F}(Q). \\ \end{aligned} \end{aligned}$$
(10)

In this case, we propose a linear reformulation of the objective function \(-I_r\) for (10).

Linear reformulation of \(-I_r\) In Charnes and Cooper (1962), the authors propose a linear reformulation for the continuous non-linear problem (NLP) of maximizing a linear fractional function subject to linear constraints. Similarly, we reformulate the linear fractional objective function \(-I_r\) in (10), and introduce new variables \(\tilde{t} \in \mathbb {R}_+^*\) and \(\tilde{w} \in \mathbb {R}^{n_{\text {CNV}}+n_{\text {CNP}}}\) to produce the following MILP:

$$\begin{aligned} \begin{aligned}&\underset{\tilde{x}=(\tilde{z},\tilde{q},\tilde{h},\tilde{\eta },\tilde{\theta },\tilde{t},\tilde{w})}{\text {minimize}}&-\sum _{k=1}^{n_t} \left( (d^k)^T \tilde{h}^k - \tilde{t} (d^k)^T h^{\min } \right) \\&\text {subject to}&R^k \tilde{q}^k + E^k \tilde{\theta }^k \le r^k t,\ \forall k=1,\dots ,n_t, \\&&(\tilde{z},\tilde{q},\tilde{h},\tilde{\theta },\tilde{\eta },\tilde{t},\tilde{w}) \in \tilde{\text {F}}(Q),\\ \end{aligned} \end{aligned}$$
(11)

where \(\tilde{q}=(\tilde{q}^k)_{k=1,\dots ,n_t}\), \(\tilde{h}=(\tilde{h}^k)_{k=1,\dots ,n_t}\), \(\tilde{\eta }=(\tilde{\eta }^k)_{k=1,\dots ,n_t}\), \(\tilde{\theta }=(\tilde{\theta }^k)_{k=1,\dots ,n_t}\), and \(\tilde{\text {F}}(Q) \subset \{0;1\}^{n_{\text {CNV}}+n_{\text {CNP}}} \times \mathbb {R}^{n_p} \times \mathbb {R}^{n_n} \times \mathbb {R}^{n_p} \times \mathbb {R}^{n_p} \times \mathbb {R} \times \mathbb {R}^{n_{\text {CNV}}+n_{\text {CNP}}}\) is the set defined by the linear constraints (C.3), (C.5a)–(C.5d), (C.6), (C.8), (C.9), (C.10) and (C.11) in “Appendix C”.

We can then show that

Proposition 1

If \(x=(z,q,h,\eta ,\theta )\) is a feasible solution of (10), then \(\tilde{x}=(\tilde{z},\tilde{q},\tilde{h},\tilde{\eta },\tilde{\theta },\tilde{t},\tilde{w})\), given by \(\tilde{t}=\left( \sum _{k=1}^{n_t} (-A_{10} h_0^k)^T q^k - (d^k)^T h^{\min } \right) ^{-1}\) and \(\tilde{z}=z,\ \tilde{q}=\tilde{t}q,\ \tilde{h}=\tilde{t}h,\ \tilde{\eta }=\tilde{t}\eta ,\ \tilde{\theta }=\tilde{t}\theta \ {and}\ \tilde{w}=\tilde{t}z\), is a feasible solution of (11). Conversely, if \(\tilde{x}=(\tilde{z},\tilde{q},\tilde{h},\tilde{\eta },\tilde{\theta },\tilde{t},\tilde{w})\) is a feasible solution of (11), \(x=(z,q,h,\eta ,\theta )\) is a feasible solution of (10), where \(z=\tilde{z},\ q=\tilde{q}/\tilde{t},\ h=\tilde{h}/\tilde{t},\ \eta =\tilde{\eta }/\tilde{t}\ {and}\ \theta =\tilde{\theta }/\tilde{t}\).

Proof

See Proposition 1 in “Appendix C”.

In addition, it follows from the above proposition that

Corollary 1

If x is a global optimal solution of (10), then \(\tilde{x}\), as given by Proposition 1, is a global optimal solution of (11), and vice versa.

Proof

See Corollary 1 in “Appendix C”.

For the sake of brevity, we refer the reader to the “Appendix” for the full description of the MILP reformulation and technical proofs of Proposition 1 and Corollary 1. Hence, to solve the linear fractional problem (10), it is necessary and sufficient to solve the corresponding linear reformulation (11). Note that a similar reformulation can be applied to problems \(\text {MINLP}_{-I_r,AZP,\epsilon _j}\), \(j=1,\dots ,n_{\epsilon }\).

As a result, we are able to compute lower bounds L(Q) on the optimal value of subproblems \(\text {MINLP}(Q)\) throughout the sBB search by relaxing potential-flow flow constraints (2b), reformulating the resilience index \(-I_r\) if necessary, and applying an MILP solver to the resulting linear problem.

3.2.3 Upper bounding

Finally, upper bounds \(\text {U}(Q)\) on the optimal value of subproblems \(\text {MINLP}(Q)\) are computed throughout the sBB search by fixing the values of the binary variables z in \(\text {MINLP}(Q)\) to the solution of the corresponding lower bounding MILP relaxations and computing a locally optimal solution to the resulting non-convex NLPs. This is done directly using an NLP solver.

figure d

3.2.4 Note on the initialization of the sBB search for \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\)

For \(j=1,\dots ,n_{\epsilon }-1\), let \(x_{f_1,f_2,\epsilon _j}\) be the best feasible solution of \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\) returned by the sBB procedure. Then, if \(\epsilon _{j+1}\) is greater than or equal to \(\epsilon _j\), since \(x_{f_2}\) is the best feasible solution found for \(\text {MINLP}_{f_2}(Q)\), we have that \(f_2(x_{f_1})-f_2(x_{f_2}) \ge 0\), and

$$\begin{aligned} \begin{aligned} f_2(x_{f_1,f_2,\epsilon _j})&\le f_2(x_{f_2}) + \epsilon _j (f_2(x_{f_1})-f_2(x_{f_2})), \\&\le f_2(x_{f_2}) + \epsilon _{j+1} (f_2(x_{f_1})-f_2(x_{f_2})). \end{aligned} \end{aligned}$$

Since all other constraints are the same, \(x_{f_1,f_2,\epsilon _j}\) is a feasible solution of \(\text {MINLP}_{f_1,f_2,\epsilon _{j+1}}(Q)\). As a result, we define the sequence of parameters \((\epsilon _j)_{j=1,\dots ,n_{\epsilon }}\) such that \(\epsilon _j \le \epsilon _{j+1}\), for all \(j=1,\dots ,n_{\epsilon }-1\), and we apply sBB to the sequence of parametrized problems \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\) for \(j=1,\dots ,n_{\epsilon }\), in increasing order of \(\epsilon _j\), initializing the upper bound UB of every problem \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\), \(j=1,\dots ,n_{\epsilon }\), at the best feasible solution of the previous problem in the sequence, i.e. \(\text {UB} \leftarrow \min (\text {U}(Q),f_1(x_{f_1,f_2,\epsilon _{j-1}}))\) on line 4 in Algorithm 1. In particular, \(x_{f_1,f_2,\epsilon _{0}}=x_{f_2}\), the best feasible solution of \(\text {MINLP}_{f_2}(Q)\).

3.3 Approximation of the Pareto front

In this section, we show that, combining the methods described in Sects. 3.1 and 3.2, we obtain a superset of the Pareto front of \(\text {BOMINLP}(Q)\). We define \(\mathcal {L}^{\text {FS}}_{f_1,f_2}\) as the image in the objective space of the set of feasible solutions of \(\text {MINLP}_{f_i}(Q)\), \(i=1,2\), and \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\), \(j=1,\dots ,n_{\epsilon }\), returned by Algorithm 1:

$$\begin{aligned} \begin{aligned} \mathcal {L}^{\text {FS}}_{f_1,f_2} =\ \left( \bigcup _{i=1,2} \{ (f_1(x_{f_i}),f_2(x_{f_i})) \} \right) \cup \left( \bigcup _{j=1\dots ,n_{\epsilon }} \{ (f_1(x_{f_1,f_2,\epsilon _j}),f_2(x_{f_1,f_2,\epsilon _j})) \} \right) . \\ \end{aligned} \end{aligned}$$

Next, we perform a pairwise comparison of the elements of \(\mathcal {L}^{\text {FS}}_{f_1,f_2}\), called a Pareto filter, to eliminate dominated solutions. The resulting subset of potentially non-dominated solutions, \(\mathcal {L}^{\text {NDS}}_{f_1,f_2}\), is then given by:

$$\begin{aligned} \begin{aligned} \mathcal {L}^{\text {NDS}}_{f_1,f_2} =\ \bigcup _{p=(p_{f_1},p_{f_2}) \in \mathcal {L}^{\text {FS}}_{f_1,f_2}} \{ p : \forall \tilde{p}=(\tilde{p}_{f_1},\tilde{p}_{f_2}) \in \mathcal {L}^{\text {FS}}_{f_1,f_2},\ p_{f_1}<\tilde{p}_{f_1}\ \text {or}\ p_{f_2}<\tilde{p}_{f_2} \}. \\ \end{aligned} \end{aligned}$$

For a detailed description of the Pareto filter algorithm, see Section 4.1 in Messac et al. (2003). We also define the sets \(\mathcal {L}^{\text {LBS}}_{f_1,f_2}\), \(T_{f_1,f_2}\) and \(U_{f_1,f_2}\):

$$\begin{aligned} \begin{aligned} \mathcal {L}^{\text {LBS}}_{f_1,f_2}&= \bigcup _{j=1,\dots ,n_{\epsilon }} \{ (\text {LB}_{f_1,f_2,\epsilon _j},f_2(x_{f_2}) + \epsilon _j (f_2(x_{f_1})-f_2(x_{f_2}))) \}, \\ T_{f_1,f_2}&= \bigcup _{(p_{f_1},p_{f_2}) \in \mathcal {L}^{\text {NDS}}_{f_1,f_2}} \{ (p_{f_1}+u,p_{f_2}+v) : (u,v) \in (\mathbb {R}_{+}^*)^2 \}, \\ \text {and} \quad U_{f_1,f_2}~&=\ U_1 \cup U_2, \end{aligned} \end{aligned}$$

where

$$\begin{aligned} \begin{aligned} U_1&= \{ (u,v) \in \mathbb {R}^2 : u< \text {LB}_{f_1}\ \text {or}\ v < \text {LB}_{f_2} \}, \\ U_2&= \bigcup _{(p_{f_1},p_{f_2}) \in \mathcal {L}^{\text {LBS}}_{f_1,f_2}} \{(p_{f_1}+u,p_{f_2}+v) : u \in \mathbb {R}_{-}^*, v \in \mathbb {R}_{-} \}. \end{aligned} \end{aligned}$$

We refer the reader to Figs. 1 and 3 for an illustration of the sets \(\mathcal {L}^{\text {LBS}}_{f_1,f_2}\), \(\mathcal {L}^{\text {NDS}}_{f_1,f_2}\), \(T_{f_1,f_2}\) and \(U_{f_1,f_2}\). Then, the following Proposition holds.

Proposition 2

The set \(T_{f_1,f_2} \cup U_{f_1,f_2}\) does not contain non-dominated solutions of \(\mathrm {BOMINLP}(Q)\).

Proof

Let \(\tilde{p}=(\tilde{p}_{f_1},\tilde{p}_{f_2})\) be a non-dominated solution of \(\text {BOMINLP}(Q)\), corresponding to the efficient decision vector \(\tilde{x}\), i.e. \(\tilde{p}_{f_1}=f_1(\tilde{x})\) and \(\tilde{p}_{f_2}=f_2(\tilde{x})\). We first show, by contradiction, that \(\tilde{p} \notin T_{f_1,f_2}\). Assume \(\tilde{p} \in T_{f_1,f_2}\). Then \(\exists \ p=(p_{f_1},p_{f_2}) \in \mathcal {L}^{\text {NDS}}_{f_1,f_2}\ \text {and}\ (u,v) \in (\mathbb {R}^*_{+})^2\ :\ \tilde{p}_{f_1}=p_{f_1}+u,\ \tilde{p}_{f_2}=p_{f_2}+v\). Hence, \(\tilde{p}_{f_1}<p_{f_1}\) and \(\tilde{p}_{f_2}<p_{f_2}\). Since p is a feasible solution of \(\text {BOMINLP}(Q)\) (by definition of \(\mathcal {L}^{\text {NDS}}_{f_1,f_2}\)), this contradicts the assumption that \(\tilde{p}\) is a non-dominated solution of \(\text {BOMINLP}(Q)\). We can conclude that \(\tilde{p} \in \mathbb {R}^2 {\setminus } T_{f_1,f_2}\).

Next, we show that \(\tilde{p} \notin U_{f_1,f_2}\). Assume \(\tilde{p} \in U_{f_1,f_2}\). We distinguish three cases.

  • There exists \(p=(p_{f_1},p_{f_2}) \in \mathcal {L}_{f_1,f_2}^{\text {LBS}},\ u \in \mathbb {R}_{-}^*\ \text {and}\ v \in \mathbb {R}_{-} :\ \tilde{p}_{f_1}=p_{f_1}+u,\ \tilde{p}_{f_2}=p_{f_2}+v\), i.e. \(\tilde{p}_{f_1}<p_{f_1}\) and \(\tilde{p}_{f_2} \le p_{f_2}\). By definition of \(\mathcal {L}_{f_1,f_2}^{\text {LBS}}\), \(\exists \ j=1,\dots ,n_{\epsilon }\ :\ p_{f_2}=f_2(x_{f_2}) + \epsilon _j (f_2(x_{f_1})-f_2(x_{f_2}))\). Moreover, \(\forall \ x\) feasible solutions of \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\) with parameter value \(\epsilon _j,\ p_{f_1}\le f_1(x)\). As \(\tilde{p}_{f_2} \le p_{f_2}=f_2(x_{f_2}) + \epsilon _j (f_2(x_{f_1})-f_2(x_{f_2}))\), \(\tilde{x}\) is a feasible solution of the parametrized problem \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\) with parameter value \(\epsilon _j\); This however contradicts \(p_{f_1}\) lower bound of \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\), as \(\tilde{p}_{f_1}<p_{f_1}\).

  • We have \(\tilde{p}_{f_1} < \text {LB}_{f_1}\). Since \(\tilde{x}\) is a feasible solution of \(\text {BOMINLP}(Q)\), this contradicts \(\text {LB}_{f_1}\) lower bound on the optimal value of \(f_1\): by definition of \(\text {LB}_{f_1}\), \(\forall \ x\) feasible solution of \(\text {BOMINLP}(Q)\), \(\text {LB}_{f_1} \le f_1(x)\);

  • We have \(\tilde{p}_{f_2} < \text {LB}_{f_2}\). As \(\tilde{p}\) is a feasible solution of \(\text {BOMINLP}(Q)\), this contradicts \(\text {LB}_{f_2}\) lower bound on the optimal value of \(f_2\).

As a result, \(\tilde{p} \in \mathbb {R}^2 {\setminus } \left( T_{f_1,f_2} \cup U_{f_1,f_2} \right) \).

Proposition 2 implies that the set \(\mathbb {R}^2 {\setminus } \left( T_{f_1,f_2} \cup U_{f_1,f_2} \right) \), defined by the best feasible solutions and lower bounds returned for \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\), \(j=1,\dots ,n_{\epsilon }\), by Algorithm 1, is a superset of the non-dominated set of \(\text {BOMINLP}(Q)\). In particular, as Proposition 2 holds for both \(f_1=\text {AZP}\) and \(f_2=-I_r\), and \(f_1=-I_r\) and \(f_2=\text {AZP}\), the sets T, U, \(\mathcal {L}^{\text {NDS}}\) and \(\mathcal {L}^{\text {LBS}}\), given by

$$\begin{aligned} \begin{aligned}&T =\ T_{\text {AZP},-I_r} \cup T_{-I_r,\text {AZP}}, \\&U =\ U_{\text {AZP},-I_r} \cup U_{-I_r,\text {AZP}}, \\&\mathcal {L}^{\text {NDS}} =\ \left( \mathcal {L}^{\text {NDS}}_{\text {AZP},-I_r} \cup \mathcal {L}^{\text {NDS}}_{-I_r,\text {AZP}} \right) {\setminus } T, \\ \text {and}&\mathcal {L}^{\text {LBS}} =\ \left( \mathcal {L}^{\text {LBS}}_{\text {AZP},-I_r} \cup \mathcal {L}^{\text {LBS}}_{-I_r,\text {AZP}} \right) {\setminus } U, \\ \end{aligned} \end{aligned}$$

define a tighter superset of the Pareto front of \(\text {BOMINLP}(Q)\), corresponding to the intersection of the supersets \(\mathbb {R}^2 {\setminus } \left( T_{\text {AZP},-I_r} \cup U_{\text {AZP},-I_r} \right) \) and \(\mathbb {R}^2 {\setminus } \left( T_{-I_r,\text {AZP}} \cup U_{-I_r,\text {AZP}} \right) \).

Similarly, the sets of explored nodes returned by Algorithm 1 for problems \(\text {MINLP}_{f_i}(Q)\), \(i=1,2\), and \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\), \(j=1,\dots ,n_{\epsilon }\), define a superset of the efficient set of \(\text {BOMINLP}(Q)\)—see “Appendix D”. The quality of the global efficiency bounds on solutions \(x_i\), \(i=1,2\), and \(x_{\epsilon _j}\), \(j=1,\dots ,n_{\epsilon }\) then depend on the size of the sets F(Q) in the tree of explored nodes returned by Algorithm 1.

4 Application and results

We investigate the optimal design-for-control of two WDNs, namely Net25 and pescara. For each network, we consider the bi-objective optimization problem \(\text {BOMINLP}(Q)\), with the objectives of minimizing average zone pressure and maximizing network resilience, and apply the method described in Sect. 3. As the scalarization method of \(\epsilon \)-constraints is sensitive to the order of objectives functions, we propose to solve both sequences of parametrized problems \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}\) and \(\text {MINLP}_{-I_r,\text {AZP},\epsilon _j}\), for \(j=1,\dots ,n_{\epsilon }\). We consider different values of \(n_{\epsilon }\), and for \(j=1,\dots ,n_{\epsilon }\), we define \(\epsilon _j=j/n_{\epsilon }\).

For both Net25 and pescara, we first solve the anchor point problems \(\text {MINLP}_{\text {AZP}}(Q)\) and \(\text {MINLP}_{-I_r}(Q)\) and compare the performance of Algorithm 1 on these instances to the off-the-shelf solvers SCIP and BARON. For each solution method, the time limit is fixed to 5 min (\(T_{\max }=5\) min), including pre-processing steps. In particular, all bound tightening steps in Algorithm 1 are performed within the fixed time limit. The results of the comparison, presented in Tables 3, 4, 5 and 6 in “Appendix E”, show that the spatial branch-and-bound algorithm outperforms state-of-the-art global optimization solvers on all instances of \(\text {MINLP}_{\text {AZP}}(Q)\) and \(\text {MINLP}_{-I_r}(Q)\). Moreover, a posteriori analysis of the convergence of Algorithm 1 for the solution of the anchor point problems \(\text {MINLP}_{\text {AZP}}(Q)\) and \(\text {MINLP}_{-I_r}(Q)\) shows that both upper and lower bounds improve rapidly at first. Past 5 min, however, the improvement in the optimality gap slows down significantly—see “Appendix E”. As a result, in the rest of this section, we implement Algorithm 1, with \(T_{\max }=5\) min, for the solution of the sequence of scalarized problems \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\), \(j=1,\dots ,n_{\epsilon }\), and define the sets \(T_{f_1,f_2}\) and \(U_{f_1,f_2}\) based on the best feasible solutions and lower bounds returned. Within Algorithm 1, GUROBI (v8.1) (Gurobi Optimization 2017) is used to solve the lower bounding and bound tightening MILPs, while upper bounds are computed with IPOPT (v3.12.9) (Wachter and Biegler 2006) using the Matlab interface (Currie et al. 2012). Taking advantage of the sparsity of the NLPs, we directly provide pre-computed gradients and Jacobians to the solver. The tolerance and time limit in Algorithm 1 are set to \(10^{-2}\) and 5 min. All computations are performed using MATLAB 2018b for Windows, on a 2.26GHz Intel® Xeon® CNPU E5520 with 8 cores.

4.1 Net25

Sectorization is a common pressure and leakage management practice which consists in closing off links in operational networks using boundary valves (BV) to create isolated supply areas, called district metered areas (DMA). Although sectorization reduces the topological redundancy of WDNs, improvements in resilience can be achieved at minimal cost in sectorized networks by re-opening closed BVs and aggregating DMAs. We illustrate this on the modified network Net25 from Ulusoy et al. (2019), where nodes with positive demand represent DMAs. We formulate the optimal aggregation of isolated DMAs (represented in Net25 by nodes with positive demand) as a DfC problem, where the set of CNPs correspond to closed BVs and DMA aggregation is limited to pairing, to preserve the structure of the network. The selection of CNPs for installation can then lead to the aggregation of a DMA with at most one other DMA, and the pairing of two DMAs excludes pairing opportunities with other neighboring DMAs.

The network model, provided in the supplementary material, consists of 3 reservoirs, 17 nodes and 26 links (15 pipes, 3 existing valves and 8 closed BVs, representing CNPs for installation). Friction head losses across network links are modelled using the H–W equation, and we consider 24 time steps. The maximum allowed velocity is fixed to 2 m/s for all network pipes (based on preliminary hydraulic simulations of the network) and a quadratic approximation \(\phi (\cdot )\) of the H–W formula is adopted to model head loss across the pipes (see Sect. 2). We set the pressure requirement at nodes with positive demand to 10 m and model the pairing limit on DMA aggregation with constraint (6). For the details of the computation of matrix A corresponding to the constraint of DMA pairing, see Appendix A in Ulusoy et al. (2019). Note that the combinatorial constraint limits the number of CNPs which can be simultaneously selected for installation.

Anchor points The computation of the anchor points of \(\text {BOMINLP}(Q)\) for Net25 returns \(-I_r(x_{-I_r})=-0.85\), \(\text {AZP}(x_{-I_r})=38.5\), \(-I_r(x_{\text {AZP}})=-0.16\), \(\text {AZP}(x_{\text {AZP}})=18.2\), \(\text {LB}_{-I_r}=-0.86\) and \(\text {LB}_{\text {AZP}}=14.2\).

Solution of \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q)\) and \(\text {MINLP}_{-I_r,\text {AZP},\epsilon _j}(Q)\), \(j=1,\dots ,n_{\epsilon }\) We apply Algorithm 1 to the sequences of parametrized problems \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q)\) and \(\text {MINLP}_{-I_r,\text {AZP},\epsilon _j}(Q)\), for \(j=1,\dots ,n_{\epsilon }\) and \(n_{\epsilon }=10,20,100\). After discarding the dominated solutions, we obtain the sets \(\mathcal {L}_{\text {AZP},-I_r}^{\text {NDS}}\) and \(\mathcal {L}_{-I_r,\text {AZP}}^{\text {NDS}}\) of potentially non-dominated solutions, represented in Fig. 1a and b for \(n_{\epsilon }=20\). On Fig. 1a and b, we also represent the sets \(\mathcal {L}_{\text {AZP},-I_r}^{\text {LBS}}\), \(T_{\text {AZP},-I_r}\) and \(U_{\text {AZP},-I_r}\), and \(\mathcal {L}_{-I_r,\text {AZP}}^{\text {LBS}}\), \(T_{-I_r,\text {AZP}}\) and \(U_{-I_r,\text {AZP}}\), respectively (see Sect. 3.3). The staircase-shaped areas \(\mathbb {R}^2 {\setminus } \left( T_{\text {AZP},-I_r} \cup U_{\text {AZP},-I_r} \right) \) and \(\mathbb {R}^2 {\setminus } \left( T_{-I_r,\text {AZP}} \cup U_{-I_r,\text {AZP}} \right) \) then represent the supersets of the Pareto front of \(\text {BOMINLP}(Q)\) obtained with each scalarization approach. The differences between the two scalarization approaches are discussed in greater detail in Sect. 4.3. As observed in Sect. 3.3, a tighter superset of the Pareto front is obtained by taking the intersection of \(\mathbb {R}^2 {\setminus } \left( T_{\text {AZP},-I_r} \cup U_{\text {AZP},-I_r} \right) \) and \(\mathbb {R}^2 {\setminus } \left( T_{-I_r,\text {AZP}} \cup U_{-I_r,\text {AZP}} \right) \), represented on Fig. 1c. The width of the set \(\mathbb {R}^2 {\setminus } \left( T \cup U \right) \), where \(T=T_{\text {AZP},-I_r} \cup T_{-I_r,\text {AZP}}\) and \(U=U_{\text {AZP},-I_r} \cup U_{-I_r,\text {AZP}}\), suggests that the potentially non-dominated solutions obtained for \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q)\) and \(\text {MINLP}_{-I_r,\text {AZP},\epsilon _j}(Q)\), \(j=1,\dots ,20\), provide a good approximation of the true non-dominated set of \(\text {BOMINLP}(Q)\).

Fig. 1
figure 1

Solution of \(\text {BOMINLP}(Q)\) for Net25, where the combinatorial constraint (6) limits the aggregation of DMAs to pairing (Ulusoy et al. 2019). Network designs (sets of selected CNPs) corresponding to the solutions in \(\mathcal {L}^{\text {NDS}}\) are provided in the supplementary material

Full size image

Next, we discuss the trade-off between the objectives of AZP and \(I_r\), their evolution along the Pareto front and the hydraulic characteristics of the configurations corresponding to the potentially non-dominated solutions represented in Fig. 1c. Potentially non-dominated solutions \(p \in \mathcal {L}^{\text {NDS}}\) with lower AZP values, represented in the top left corner in Fig. 1c, correspond to more sectorized network designs, resulting in greater friction head loss, and the implementation of aggressive pressure reduction. Such configurations minimize surplus pressure at customer nodes, resulting in lower AZP and \(I_r\) values. Moving down along the Pareto front in Fig. 1c towards more resilient network configurations in the bottom right corner, pressure reduction is relaxed and more CNPs are selected for installation. This results in better flow distribution, reduced friction head loss across the network, excess pressure at customer nodes and, eventually, greater AZP and \(I_r\) values.

Finally, we focus on the geometry of the Pareto front of \(\text {BOMINLP}(Q)\), which is expected to be disconnected and/or composed of different branches. Figure 1d represents the potentially non-dominated solutions \(p \in \mathcal {L}^{\text {NDS}}\) from Fig. 1c, where solutions corresponding to different network designs are represented by different colors. Note that several points in \(\mathcal {L}^{\text {NDS}}\) correspond to the same network design: \(\mathcal {L}^{\text {NDS}}\) represents the union of 6 fronts, corresponding to 6 different network designs, i.e. sets of CNPs selected for addition. As a result, when selecting the most suitable network design, the decision-maker is not limited to a single trade-off point, but rather a set of trade-offs corresponding to different PCV control settings. For example, if the decision maker’s priority is minimizing AZP, then design 5, which corresponds to opening a single BV in the network, represents a better trade-off than design 1, which results in \(8\%\) larger minimum AZP values. On the other hand, design 1, which corresponds to the creation of 3 pairs of DMAs, results in more resilient network configurations.

4.2 pescara

Fig. 2
figure 2

Modified pescara network model from Bragalli et al. (2008). The original network model is available at http://www.or.deis.unibo.it/research_pages/ORinstances/ORinstances.htm

Full size image

The pescara model is a reduced version of the WDN of an Italian medium-size city, studied by Bragalli et al. (2008). The original network model counts 71 junctions, 3 reservoirs, 99 pipes and is defined on a single time step. Friction losses are modeled according to the H–W equation and maximum allowed velocities are fixed to 3 m/s for all network links. We investigate the problem of expanding the network by installing new pressure control valves and pipes in predefined candidate locations (see Fig. 2 and supplementary material), determined based on engineering judgement. For each candidate link location, we consider 3 CNPs with diameters 100, 150 and 200 mm. The material of the CNPs is cast iron, as for existing pipes, and we include a combinatorial constraint in the optimization problem in the form (6) to limit the selection of CNPs to 1 per location.

Anchor points The computation of the anchor points of \(\text {BOMINLP}(Q)\) returns \(-I_r(x_{-I_r})=-0.45\), \(\text {AZP}(x_{-I_r})=27.9\), \(-I_r(x_{\text {AZP}})=-0.18\), \(\text {AZP}(x_{\text {AZP}})=16.5\), \(\text {LB}_{-I_r}=-0.46\) and \(\text {LB}_{\text {AZP}}=15.2\).

Solution of \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q)\) and \(\text {MINLP}_{-I_r,\text {AZP},\epsilon _j}(Q)\), \(j=1,\dots ,n_{\epsilon }\) We apply Algorithm 1 to the sequences of parametrized problems \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q)\) and \(\text {MINLP}_{-I_r,\text {AZP},\epsilon _j}(Q)\), for \(j=1,\dots ,n_{\epsilon }\) and \(n_{\epsilon }=10,20,100\). For each approach, dominated solutions \(p \in \mathcal {L}^{\text {FS}}_{\text {AZP},-I_r}\) or \(\mathcal {L}^{\text {FS}}_{-I_r,\text {AZP}}\) are discarded. The resulting subsets of potentially non-dominated solutions \(\mathcal {L}^{\text {NDS}}_{\text {AZP},-I_r}\) and \(\mathcal {L}^{\text {NDS}}_{-I_r,\text {AZP}}\) obtained for \(n_{\epsilon }=20\) are represented in Fig. 3a and b, along with the sets \(\mathcal {L}^{\text {LBS}}_{\text {AZP},-I_r},\ T_{\text {AZP},-I_r}\) and \(U_{\text {AZP},-I_r}\), and \(\mathcal {L}^{\text {LBS}}_{-I_r,\text {AZP}},\ T_{-I_r,\text {AZP}}\), and \(U_{-I_r,\text {AZP}}\). We comment on the differences between the two approaches in Sect. 4.3. Moreover, we represent, on Fig. 3c, the sets \(T=T_{\text {AZP},-I_r} \cup T_{-I_r,\text {AZP}}\) and \(U=U_{\text {AZP},-I_r} \cup U_{-I_r,\text {AZP}}\) obtained for \(n_{\epsilon }=20\). The resulting staircase-shaped set \(\mathbb {R}^2 {\setminus } \left( T \cup U \right) \), which corresponds to the intersection of the sets \(\mathbb {R}^2 {\setminus } \left( T_{\text {AZP},-I_r} \cup U_{\text {AZP},-I_r} \right) \) and \(\mathbb {R}^2 {\setminus } \left( T_{-I_r,\text {AZP}} \cup U_{-I_r,\text {AZP}} \right) \) represented in Fig.  3a and b, provides a tighter superset of the set of non-dominated solutions of \(\text {BOMINLP}(Q)\).

Fig. 3
figure 3

Solution of \(\text {BOMINLP}(Q)\) for pescara. Network designs (sets of selected CNVs and CNPs) corresponding to the solutions in \(\mathcal {L}^{\text {NDS}}\) are provided in the supplementary material

Full size image

Potentially non-dominated solutions with lower AZP values in Fig. 3c correspond to the installation of fewer or smaller CNPs and more CNVs, and the implementation of aggressive pressure reduction, which result in greater head loss across the network, reduced surplus pressure at customer nodes and lower \(I_r\) values. In particular, solutions in the top left corner of Fig. 3c require the installation of all CNVs. Moving down along the Pareto front in Fig. 3c towards the bottom right corner, more resilient potentially non-dominated solutions \(p \in \mathcal {L}^{\text {NDS}}\) correspond to the installation of fewer CNVs and CNPs with larger diameters, which increase the number and capacity of alternative customer supply paths. For network configurations corresponding to solutions \(p \in \mathcal {L}^{\text {NDS}}\) in the intermediate segment of the Pareto front, only the CNV located at the outlet of the source of greatest fixed hydraulic head \(h_0\) is selected for installation (and the implementation of pressure reduction). The effect is to increase the relative importance of other supply sources, improving flow distribution across the network. As the energy efficiency of the network grows, so does the pressure surplus at network nodes, increasing AZP and the resilience index, \(I_r\).

Design variables in the optimal design-for-control problem of pescara include the placement of PCVs, and 3 different diameters to choose from, for each candidate pipe. As this results in a higher number of possible design combinations than for Net25, it is expected that fewer potentially non-dominated solutions \(p \in \mathcal {L}^{\text {NDS}}\) correspond to the same network design. From the perspective of the decision-maker, it might be interesting to take more objectives (e.g. cost) into account when deciding between two design solutions offering similar trade-offs between AZP and \(-I_r\).

4.3 Scalarization approach

Scalarization methods, including the method of \(\epsilon \)-constraints, are known to be sensitive to the order of the objective functions (Logist and Van Impe 2012) and the choice of scalarization parameters (Kim and De Weck 2005; Burachik et al. 2017). In the method we propose in Sect. 3, the first objective function \(f_1\) is chosen to be minimized while parametrized upper bound constraints are enforced on the value of \(f_2\). As a result, the quality of the Pareto front approximation depends on the order of the objective functions—see Figs. 1 and 3. In this section, we discuss the differences between the results presented in Tables 1 and 2, obtained for Net25 and pescara with both sequences of parametrized problems \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q)\) and \(\text {MINLP}_{-I_r,\text {AZP},\epsilon _j}(Q)\), for \(j=1,\dots ,n_{\epsilon }\) and \(n_{\epsilon }=10,20,100\). We also comment on the choice of the uniform sequence of scalarization parameters \(\epsilon _j=j/n_{\epsilon }\), for \(j=1,\dots ,n_{\epsilon }\).

4.3.1 Order of objective functions

To evaluate the scalarization approaches, we compare the sets of potentially non-dominated solutions \(\mathcal {L}^{\text {NDS}}_{\text {AZP},-I_r}\) and \(\mathcal {L}^{\text {NDS}}_{-I_r,\text {AZP}}\), and the sets of lower bounds \(\mathcal {L}^{\text {LBS}}_{\text {AZP},-I_r}\) and \(\mathcal {L}^{\text {LBS}}_{-I_r,\text {AZP}}\) obtained by solving the sequence of parametrized problems \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q)\) and \(\text {MINLP}_{-I_r,\text {AZP},\epsilon _j}(Q)\). Let \(\hat{\mathcal {L}}^{\text {NDS}}_{\text {AZP},-I_r}\) and \(\hat{\mathcal {L}}^{\text {NDS}}_{-I_r,\text {AZP}}\) be subsets of \(\mathcal {L}^{\text {NDS}}_{\text {AZP},-I_r}\) and \(\mathcal {L}^{\text {NDS}}_{-I_r,\text {AZP}}\) defined by:

where \(\mathcal {L}^{\text {NDS}}\) is defined as in Sect. 3.3. Then, the greater the cardinalities \(n(\hat{\mathcal {L}}^{\text {NDS}}_{\text {AZP},-I_r})\) and \(n(\hat{\mathcal {L}}^{\text {NDS}}_{-I_r,\text {AZP}})\) of the sets \(\hat{\mathcal {L}}^{\text {NDS}}_{\text {AZP},-I_r}\) and \(\hat{\mathcal {L}}^{\text {NDS}}_{-I_r,\text {AZP}}\), the more potentially non-dominated solutions of \(\text {BOMINLP}(Q)\) are obtained with the corresponding scalarization approach. Similarly, we define the subsets \(\hat{\mathcal {L}}^{\text {LBS}}_{\text {AZP},-I_r}\) and \(\hat{\mathcal {L}}^{\text {LBS}}_{-I_r,\text {AZP}}\) of \(\mathcal {L}^{\text {LBS}}_{\text {AZP},-I_r}\) and \(\mathcal {L}^{\text {LBS}}_{-I_r,\text {AZP}}\) by:

where \(\mathcal {L}^{\text {LBS}}\) is defined as in Sect. 3.3. The greater the cardinalities \(n(\hat{\mathcal {L}}^{\text {LBS}}_{\text {AZP},-I_r})\) and \(n(\hat{\mathcal {L}}^{\text {LBS}}_{-I_r,\text {AZP}})\) of the sets \(\hat{\mathcal {L}}^{\text {LBS}}_{\text {AZP},-I_r}\) and \(\hat{\mathcal {L}}^{\text {LBS}}_{-I_r,\text {AZP}}\), the tighter the lower bound on the Pareto front of \(\text {BOMINLP}(Q)\) obtained with the corresponding scalarization approach.

Table 1 Results of the solution of the sequence of parametrized problems \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q)\) and \(\text {MINLP}_{-I_r,\text {AZP},\epsilon _j}(Q)\), \(j=1,\dots ,n_{\epsilon }\) for the network Net25 and \(n_{\epsilon }=10,20,100\)
Full size table
Table 2 Results of the solution of the sequence of parametrized problems \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q)\) and \(\text {MINLP}_{-I_r,\text {AZP},\epsilon _j}(Q)\), \(j=1,\dots ,n_{\epsilon }\) for the network pescara and \(n_{\epsilon }=10,20,100\)
Full size table

Tables 1 and 2 show that the quality of the Pareto front approximation obtained by solving the sequence of scalarized problems \(\text {MINLP}_{f_1,f_2,\epsilon _j}(Q)\), \(j=1,\dots ,n_{\epsilon }\), depends on the order of the objective functions, AZP and \(-I_r\). In particular, the solutions of the sequence of parametrized problems \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q)\), \(j=1,\dots ,n_{\epsilon }\), define a tighter superset of the non-dominated set of \(\text {BOMINLP}(Q)\) over the range of values taken by \(\text {AZP}\) and \(-I_r\), as \(n(\hat{\mathcal {L}}^{\text {NDS}}_{\text {AZP},-I_r}) \ge n(\hat{\mathcal {L}}^{\text {NDS}}_{-I_r,\text {AZP}})\) and \(n(\hat{\mathcal {L}}^{\text {LBS}}_{\text {AZP},-I_r}) \ge n(\hat{\mathcal {L}}^{\text {LBS}}_{-I_r,\text {AZP}})\), for both networks and for all values of \(n_{\epsilon }\). This is illustrated, for pescara, by the difference in the width of the supersets \(\mathbb {R}^2 {\setminus } \left( T_{\text {AZP},-I_r} \cup U_{\text {AZP},-I_r} \right) \) and \(\mathbb {R}^2 {\setminus } \left( T_{-I_r,\text {AZP}} \cup U_{-I_r,\text {AZP}} \right) \) represented in Fig. 3a and b.

The difference between the obtained Pareto front approximations can be explained by the number of sBB iterations completed by Algorithm 1 for the scalarized problems \(\text {MINLP}_{-I_r,\text {AZP},\epsilon _j}(Q)\) and \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q)\). For instance, for pescara and \(n_{\epsilon }=20\), a single sBB iteration takes on average 5 times longer to complete for problems \(\text {MINLP}_{-I_r,\text {AZP},\epsilon _j}(Q)\) than for \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q)\). Therefore, half as many sBB iterations, on average, are completed within the fixed time limit (\(T_{\max }=5\) min) for each parametrized problem in the former sequence, resulting in larger gaps. By fixing a limit on the maximum number of sBB iterations completed (instead of computational time) in Algorithm 1, the results obtained with the different scalarization approaches are more balanced—see “Appendix G”. However, the results presented in Tables 1 and 2 show that it is more efficient to solve the sequence of parametrized problems \((\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q))_{j=1,\dots ,n_{\epsilon }}\) to produce a tight approximation of the Pareto front of \(\text {BOMINLP}(Q)\), within a prescribed tolerance or time limit.

4.3.2 Note on the choice of scalarization parameters \(\epsilon _j\)

The results presented in this section for the optimal design-for-control of Net25 and pescara were obtained by solving the sequences of parametrized problems \(\text {MINLP}_{\text {AZP},-I_r,\epsilon _j}(Q)\) and \(\text {MINLP}_{-I_r,\text {AZP},\epsilon _j}(Q)\), for \(j=1,\dots ,n_{\epsilon }\), using a predefined, uniformly distributed set of scalarization parameters \(\epsilon _j=j/n_{\epsilon }\).

Depending on the geometry of the Pareto front, traditional scalarization approaches relying on predefined distributions of parameters might fail to explore regions of the non-dominated set, resulting in redundant points and poorly distributed solutions on the Pareto front. Traditional weighted sum or \(\epsilon \)-constraint methods for instance struggle to obtain a representative approximation of Pareto fronts with steep segments (i.e. knees), while the normalized boundary intersection (NBI) method (Das and Dennis 1998) requires the non-dominated set to be simply connected (i.e. not composed of disconnected branches) to guarantee finding a Pareto optimal solution for all parametrized problems (Burachik et al. 2017). Previous works have shown that these limitations could be successfully addressed by using adaptive choices of scalarization parameters (Kim and De Weck 2005) or grid generation techniques (Burachik et al. 2017, 2019), to focus on unexplored regions of the Pareto front and obtain representative approximations of non-dominated sets with irregular geometries.

In this case, however, we observe (see Figs. 1 and 3) that the choice of scalarization parameters \(\epsilon _j\) results in an even spread of potentially non-dominated solutions \(\mathcal {L}^{\text {NDS}}\) and lower bounds \(\mathcal {L}^{\text {LBS}}\), providing a good approximation of the whole non-dominated set of \(\text {BOMINLP}(Q)\). As a result, in this work, we do not further investigate adaptive methods for the generation of scalarization parameters \(\epsilon _j\). The implementation of grid generation techniques could however be investigated in future work, in particular for multi-objective problems with more than two objectives.

5 Conclusion

To the best of our knowledge, this study is the first to consider the bi-objective optimization for the minimization of pressure induced background leakage and maximization of resilience in water distribution networks. In particular, we investigate the bi-objective optimal design-for-control of existing WDNs, which consists in simultaneously identifying new valves and/or pipes to add from a predefined set of candidates and optimizing the control of new and existing valves in the network. The mathematical formulation of the problem results in a non-convex BOMINLP.

The exact identification or approximation of the Pareto front of non-convex BOMINLPs with global bounds are challenging tasks, which have not been previously addressed in the literature. In this work, we present a method to compute an approximation of the non-dominated set of the considered bi-objective DfC problem with global non-dominance bounds. The BOMINLP is first scalarized, using the method of \(\epsilon \)-constraints, into a sequence of parametrized single-objective non-convex MINLPs. Spatial branch-and-bound is then applied to compute feasible solutions of the individual non-convex MINLPs with certified global optimality bounds. In particular, we propose a novel linear reformulation of the resilience index, \(I_r\), which allows us to derive MILP relaxations of the parametrized MINLPs and, as a result, lower bounds on the optimal value of the parametrized problems.

We show that the implementation of sBB in combination with \(\epsilon \)-constraints returns a set of potentially non-dominated solutions as well as a superset of the Pareto front of the original BOMINLP based on the feasible solutions and guarantees of global optimality returned by sBB for the individual parametrized MINLPs. In particular, the superset of the Pareto front is a guarantee of global non-dominance of the potentially non-dominated solutions found. The proposed method and results can also be extended, in future work, to study multi-objective design-for-control problems with additional objective functions (e.g. cost, water quality). The implementation of grid generation techniques might be required to address the scalability challenges associated with the method of \(\epsilon \)-constraints for multi-objective problems with more than two objectives.

We apply this method to solve the bi-objective DfC problem for Net25 and pescara, and compute a representative set of trade-offs between the objectives of pressure induced background leakage minimization and resilience maximization. Network configurations corresponding to the installation of more candidate PCVs and the implementation of aggressive pressure reduction minimize pressure surplus at network nodes and result in lower AZP and \(I_r\) values. The installation of more (and larger) candidate pipes, on the other hand, improves network efficiency by increasing the number (and capacity) of alternative supply paths, and results in greater AZP and resilience index values. The results provide insights into the complex structure of the Pareto fronts of non-convex BOMINLPs and show that the quality of Pareto front approximations obtained by solving a sequence of \(\epsilon \)-parametrized problems depends on the order of the objective functions.