Field

The present disclosure relates to controlling hydraulic conditions within a liquid conduit system, such as a water distribution network. In particular, but not exclusively, the disclosure can provide a method and apparatus for increasing the resilience of a system, the reliability and life expectancy of individual assets and reducing the risk of failures.

Background

The optimal operation of water distribution networks (WDNs) requires the satisfaction of multiple criteria, some of which may be conflicting, in order to meet increasing water demand in a cost efficient manner. Some objectives include the reduction of leakage, improvements in network resilience, increasing the reliability and life expectancy of critical assets, optimization of water quality, and the efficient use of energy in pumping.

The control of pressure has been employed to manage leakage. A decrease in the average pressure within WDNs has been shown to result in significant reductions in leakage and bursts. The sectorization of water distribution systems (the division into discrete supply areas) has also provided benefits in this respect. Networks are subdivided into smaller sectors as shown in Figure 1 , called District Metered Areas (DMAs). The installation of kept-shut isolating valves is done to construct sectors (areas) which have a single supply inlet (a single-feed from a water transmission main) so that flow into each DMA (and out of each DMA for cascading DMAs) is continuously monitored, improving the ability to manage pressure and consequently leakage. The sectorization of WDNs in single-feed sectors (DMAs) has led to the implementation of simplistic statistical based feedback control methods such as US 20110290331 A1. Such control methods rely upon the pressure measured at the hydraulically critical point (CP) within a DMA (sector).

The introduction of DMAs has severely reduced the redundancy in network connectivity, affecting system resilience and water quality negatively. Furthermore, the installation and operation of valves to sectorize the network and reduce the average zone pressure may generate high diurnal pressure variability within a supply network. This reflects the fact that demand is stochastic, the energy losses are non-linear and the hydraulic conditions within water supply systems are not steady state, although the design of the physical components within such systems and their respective simulation models often assume steady state flow conditions. In addition, hydraulic conditions within water distribution networks are frequently quasi-steady and unsteady due to a myriad of factors such as the operations of valves, pumps, malfunctioning surge protection devices, air valves, intermittency of supply, variations in demand and occasional bursts. The action of control valves to achieve desired average zone pressures as described above can contribute to this variability; and therefore, network operators require more advanced methods of model based control.

Current pressure control methods as described by US201 10290331 A1 minimise pressure at a critical point in a network but introduce significant pressure variations in supply critical pipes (hydraulic links). Such pressure variations cause pressure-induced stress on high- consequence components of the distribution system, potentially accelerating deterioration mechanics and leading to failures. This outcome runs counter to the desired objective of current control systems for minimizing the average zone pressure. As a result, current control methods as described in US201 10290331 A1 do not always present optimal solutions, as current control methods used to maintain minimum average zone pressures may lead to consequent intra- and inter-zone variations in pressure. The CP-based feedback control of current pressure management schemes may reduce the reliability and life expectancy of high-consequence assets.

There is therefore an ongoing need to improve the management of liquid conduit systems, such as water distribution networks, in order to increase efficiency and reliability, and thus offer improved cost effectiveness for the supply system as a whole.

Wright, R., Stoianov, I., et al. "Adaptive water distribution networks with dynamically reconfigurable topology." Journal of Hydroinformatics 16.6 (2014): 1280-1301 propose an alternative approach for managing sectorized water distribution networks. This new approach includes the adaptive aggregation of sectors by dynamically changing the network connectivity and topology. The dynamic network reconfiguration is achieved by the replacement of closed boundary valves with automatic control valves that modify the network topology and continuously monitor the dynamic hydraulic conditions (Figure 2), resulting in multi-feed pressure control zones. The network is reversed to its original DMA topology during periods of low demand (i.e. 2am-4am) to monitor the minimum night flow within small discrete areas and maximise the management of leakage. The DMAs are then aggregated (i.e. 4am-2am) into larger multi-feed pressure controlled zones to maximise redundancy, reliability and resilience of supply while at the same time achieving improvements in pressure management due to smaller energy losses occurring in networks with additional connectivity. This model-based pressure control is unlike alternative approaches, where pressure is controlled simply by modulating (reducing) the output pressure of an automatic control valve based on a feedback pressure signal from a critical point (CP).

Summary

According to a first aspect, there is provided a method for controlling conditions within a liquid conduit system, comprising:

defining a zone within the liquid conduit system, wherein pressure within the zone is influenced by one or more actuator valves (for example, automatic control valves);

controlling the one or more actuator valves in dependence on a Pareto efficient solution to the minimisation of functions of the average pressure within the zone (AZP: average zone pressure) and the pressure variability within the zone (PVZ: pressure variability within the zone).

According to a second aspect, there is provided an apparatus for controlling a liquid conduit system comprising a zone containing one or more actuator valves, the apparatus comprising:

one or more processors configured to calculate a Pareto efficient solution to the minimisation of functions of the average pressure within the zone and the pressure variability within the zone and configured to control the one or more actuator valves in dependence on the Pareto efficient solution.

By adopting Pareto efficient solutions for the control of actuator valves it can be ensured that both the average zone pressure (AZP) and pressure variability (PVZ) are taken into account, thereby avoiding unnecessary excesses in either variable and/or allowing the operator to select optimal trade-offs between these variables. The actuator valves may be controlled to comply with the Pareto efficient solution. By managing and reducing the values of both variables, reductions in both assets degradation and leakage within the liquid conduit system may be achieved.

Pareto efficient solutions are those that are not "dominated" by any other possible approach, where to be dominated implies that there is another feasible solution which improves one objective without causing another to become worse.

As a consequence, an effective pressure management is provided which can consider both the minimization of average zone pressure (AZP) and pressure variability (PVZ). Pressure management for both single and multi-feed systems (i.e. systems with sectors that have multiple inflows as shown in Figure 2) can be provided by considering both the minimization of AZP and PVZ simultaneously. The pressure control can be actuated by inlet control valves, which modulate the pressure and flow at the inlet(s) to the zone(s), and boundary control valves, which modulate the pressure, flow and valve position including fully open, closed or anything in between the two at the boundaries between the sectors. In standard approaches and on rare occasions where multi-feed control zones are commissioned, these are installed and operated based on "trial and error" empirical knowledge with a high risk for sub-optimal hydraulic operation and failures. In contrast, the present disclosure proposes departing from the standard practice and considering the problem of optimizing the location of the actuators (i.e. inlet and boundary control valves), and the optimal pressure settings (AZP and PVZ) simultaneously - what is referred to as a co-design optimization problem. Note that once the co-design multiobjective optimization problem is solved, different control schemes can be applied for the operational management of the optimally placed automatic control valves in different demand and application scenarios.

Controlling the actuator valves may comprise opening or closing the valves and modulating for pressure, flow and position. Their functions can change and may also comprise reducing and sustaining pressure. Furthermore, their phased commissioning can also be optimally planned with the proposed method.

In addition to controlling the actuator valves, their location may be chosen as dependent on the Pareto efficient solution. Choosing optimal valve placement locations may comprise selecting the actuator valves from a set of valves at various locations or designing the system itself to include valves at the desired locations, for example. As such, in some preferred embodiments, the location of the one or more actuator valves is selected based on the Pareto efficient solution to the minimisation of functions of the average pressure within the zone (AZP) and the pressure variability within the zone (PVZ). By allowing the location of the actuator valves to be selected in this way, improved solutions can be achieved.

The decision maker (network operator) can further refine and choose the preferred solution from the set of Pareto efficient solutions (Pareto frontier) according to a plurality of criteria which may not be possible to formalise in mathematical terms. These preferences include, but are not limited to, the following:

o Impact of traffic on the installation locations (ease of installation and maintenance). o Permission from local government to execute the installation.

o Risk of damage to other infrastructure already installed at a particular location. o Ease of accessibility of a particular location.

o Operational risk associated with the area of a particular location.

o Additional preferences based on operators' knowledge.

Preferably, the zone comprises a plurality of conduits, each conduit being provided with a weighting which impacts the recommended amplitude of pressure variability. In this manner, conduits (such as pipes) for which it is particularly desirable to minimise pressure variability may be given a greater weighting than others. For example, the weighting may be generated based on the criticality of the conduits (consequence of failure) and/or the vulnerability of the conduits. Criticality may reflect the number of people reliant upon a particular conduit or the purpose it serves (such as supply to a hospital or school), for instance. Vulnerability may reflect the conduit's current condition and its likelihood and consequence of failure, and factors may include the conduit's age, material type, environmental conditions, internal and external loads that affect pipe failure mechanisms.

The weighting may also take into account additional objective functions such as cost of installation for each candidate location and/or TOTEX (total expenditure) cost related indicators, etc.

In a preferred approach, the Pareto efficient solution is obtained solving a series of single- objective mixed integer non-linear programs. It has been found that this approach can offer a computationally efficient technique to obtain Pareto efficient solutions to the multi- objective problem associated with both average zone pressure (AZP) and pressure variability (PVZ).

Preferably, the function of average pressure within the zone is formulated as a weighted sum of nodal pressures within the zone, where weights represent the lengths of conduit connected to corresponding nodes. This function allows ready computation of the average zone pressure and may enable ready monitoring of such pressures using sensors placed at one or more nodes.

In some preferred embodiments, the Pareto efficient solutions are obtained using evolutionary approaches based on genetic algorithms (GA). Results from GA can, for example, be used to compute good initial conditions for the preferred mathematical methods. It can also be appreciated that the invention can be implemented using computer program code. Indeed, according to a further aspect there is therefore provided a computer program product comprising computer executable instructions for carrying out the method of the first aspect. The computer program product may be a physical storage medium such as a Read Only Memory (ROM) chip. Alternatively, it may be a disk such as a Digital Versatile Disk (DVD-ROM) or Compact Disk (CD-ROM). It could also be a signal such as an electronic signal over wires, an optical signal or a radio signal such as to a satellite or the like. The invention also extends to a processor running the software or code, e.g. a computer or a distributed low power embedded systems configured to carry out the method described above.

Other optional features are defined in the dependent claims in the appended claim set. Brief Description of the Figures

A preferred embodiment of the present disclosure will now be described with reference to the accompanying drawings, in which:

Figures 1 and 2 schematically illustrate the adaptive sectorization of water distribution networks with dynamically reconfigurable topology, in which: Figure 1 shows the segregation of a supply zone into DMAs for leakage detection during a period of low demand (2am-4am); Figure 2 shows aggregation of DMAs into large pressure management zones (4am-2am) for improved resilience and pressure management; and, Figure 3 shows a process for calculating the control profiles and locations for each automatic control valve to achieve the desired minimization of average pressure and pressure variability which constitute the cumulative pressure induced stress for a specific asset (pipe, etc) or a zone. This model based control method is done for either of these network topology settings (single and multi-feed zones).

Detailed Description

Referring to Figures 1 and 2, two alternative topologies of a water distribution network (WDN) are shown comprising a plurality of single feed zones 10 (Figure 1) and/or a multi- feed zones 15 (Figure 2). The WDN is an example of a liquid conduit system. Within each zone 10 and 15, a plurality of pipes 60 act as conduits to connect nodes 50. Also provided are inlet control valves 20, boundary control valves 25 and kept-shut boundary valves 30. The inlet actuator valves 20 may control fluid flow, pressure and valve position and may comprise valves located at the inlet/s into the zone. An inlet into a zone is generally defined as a supply feed from a transmission main. The boundary control valves 25 may also control in dual-direction the fluid flow, pressure and valve position and may comprise valves located within and at the boundaries of the zone with other zones. The actuator valves 20 and 25 may control liquid flow and pressure within the system together with optimal settings of pumps and reservoirs. A plurality of sensors 40 may also be provided to monitor pressure within the system.

The actuator valves may be connected to a communications network 70. The communications network may enable control of the actuator valves by control apparatus 80, which can be remotely located. The control apparatus 80 may comprise one or more processors arranged to control the actuator valves 20 and 25 and zonal pumps 66 and booster pumps 67 to perform the method described below in order to obtain Pareto efficient solutions. In this way, the pressure characteristics of a zone that include both the average zone pressure and pressure variability within the water distribution network may be controlled.

Accordingly, optimal pressure management can be achieved by controlling the actuator valves, which regulate the pressure and flow at the inlets to the zone/s, and boundary valves which can also modulate for flow, pressure and valve position (fully open, closed or anything in between the two). The actuator valves can provide a dual directional control. As well as controlling the actuator valves once installed, the location of the actuator valves (i.e. inlet control valves and boundary control valves which may include pressure reducing and/or sustaining valves) may also be optimised. The locations and the optimal pressure settings may be considered simultaneously. This establishes a co-design optimization problem. Once the co-design multi-objective optimization problem is solved, different control schemes can be applied for deriving the operational control settings of the optimal placed control valves in different operational scenarios.

The mathematical formulation for the simultaneous optimal design and control problem presents significant challenges. It requires the solution of a complex nonlinear optimization problem with both continuous and discrete variables. In the present embodiment, nodal pressures, pipe flows and valve locations are considered as decision variables, where binary variables are needed to model the placement of valves. In addition, the inclusion of hydraulic conservation laws results in non-convex, nonlinear constraints.

The co-design problem is addressed with respect to two objectives, i.e. the minimization of average zone pressure (AZP) and pressure variability (PVZ) via the optimal placement and control of valves. Therefore, a multi-objective optimization problem is formulated to establish the trade-offs between these two criteria. The resulting problem is a non-convex multi-objective mixed integer nonlinear program (MINLP).

Although the natural aspiration is to find a feasible optimal configuration for all the objectives, typically no single solution that simultaneously optimises all conflicting objectives exists. Therefore, the mathematical notion of Pareto optimality is adopted to characterise the best compromises between conflicting objectives. A feasible point for the multiobjective optimization analysis is called "not-dominated" (and consequently do lie on the Pareto Frontier) if there is no other feasible point that improves one or more of the objectives without making any other worse. The not-dominated feasible points define the Pareto Frontier.

While multiobjective co-design optimization problems for water distribution networks may be studied using evolutionary algorithms (EAs), which can be easy to implement by coupling them with simulation software, EAs have disadvantages. Firstly, they do not guarantee optimality of the solutions, not even local optimality. Secondly, a large number of function evaluations and simulations of the hydraulic model are required in order to generate useful solutions, which do not scale well with the size of the search space. Unlike mathematical optimization approaches, evolutionary methods do not explicitly or accurately handle nonlinear constraints. Since optimal control profiles for pressure management need to satisfy strict physical, quality and economic constraints, it is desirable to obtain a method that guarantees an accurate handling of constraints. Moreover, it is preferable to adopt a scalable approach for large scale water systems, which may exploit the particular sparse structure of water distribution networks.

Therefore, the preferred embodiment adopts a different approach, focusing on techniques which approximate the image of the Pareto set through the objective functions, called Pareto frontier or front. The aim is to provide a wide and uniform distribution of points of the front so that the decision maker (network operator) can choose the Pareto optimal solution which better fits its needs. This is an a posteriori approach since the articulation of the preferences by the decision maker is made after the generation of the approximated Pareto front. Common a posteriori approaches include scalarization methods (such as those described in, for example: Marler, R.T. & Arora, J.S., 2004. Survey of multi-objective optimization methods for engineering. Structural and Multidisciplinary Optimization, 26(6), pp.369-395; Das, I. & Dennis, J.E., 1998. Normal-Boundary Intersection: A New Method for Generating the Pareto Surface in Nonlinear Multicriteria Optimization Problems. SIAM Journal on Optimization, 8(3), pp.631-657; Messac, A., Ismail-Yahaya, A. & Mattson, C.A., 2003. The normalized normal constraint method for generating the Pareto frontier. Structural and Multidisciplinary Optimization, 25(2), pp.86-98; and Kim, I.Y. & De Week, O.L., 2005. Adaptive weighted-sum method for bi-objective optimization: Pareto front generation. Structural and Multidisciplinary Optimization, 29(2), pp.149-158).

These methods parameterize the multi-objective problem into a series of single-objective optimization problems that can be solved using standard nonlinear programming techniques. A popular scalarization method is the weighted sum (WS) of the objectives . Despite its simplicity, some disadvantages of the WS approach include its inability to generate Pareto points in the non-convex parts of the Pareto front and the fact that an even spread of weights may not correspond an uniform distribution of points of the front (see, for example, Das, I. & Dennis, J.E., 1997. A closer look at drawbacks of minimizing weighted sums of objectives for Pareto set generation in multicriteria optimization problems; Structural Optimization, 14(1), pp.63-69). In the last 20 years new methods have been proposed to deal with these drawbacks, in particular, the Normal Boundary Intersection (NBI) method (see, for example Das, I. & Dennis, J.E., 1998. Normal- Boundary Intersection: A New Method for Generating the Pareto Surface in Nonlinear Multicriteria Optimization Problems. SIAM Journal on Optimization, 8(3), pp.631-657) and the Normalized Normal Constraint (NNC) method (see, for example, Messac, A., Ismail-Yahaya, A. & Mattson, C.A., 2003. The normalized normal constraint method for generating the Pareto frontier. Structural and Multidisciplinary Optimization, 25(2), pp.86- 98). These methods can be combined with efficient gradient-based algorithms to find optimal solutions even for large-scale and highly constrained optimization problems.

So, the preferred embodiment adopts an approach designed to minimise average zone pressure (AZP) and pressure variability (PVZ) within the zones 10 and 15 (Figure 2). The approach adopted is described in more detail below.

1 Average Zone Pressure

Average zone pressure (AZP) according to the preferred embodiment is formulated as the weighted sum of nodal pressure, where weights represent the lengths of pipes connected to the corresponding nodes, and can be expressed using the weighted sum: k

(1) k=l i=l where a _t ^■ ■=∑ _je i ₍ i ₎ - ^an d K is ^tne set of indices for pipes incident at node i counted only once.

The normalization factor A = is chosen so that, for each time step, ^ ^ j is a convex combination of pressures.

2 Pressure variability

Pressure varia

This function encodes pressure variability indicators that have been linked with pipe failures and cumulative pressure induced stress. The function in (2) is a quadratic form and can be written as x ^T Qx, with Q e R ^NxN symmetric positive semidefinite.

Moreover, since the quadratic form only depends on pressure head values, we have that x ^T Qx = p ^T Q _p p with p = [p ¹ , ... , p ⁿⁱ ] ^T and Q _p e R ⁿ i ⁿ n ^xn i ⁿ n symmetric positive semidefinite. A preliminary numerical analysis on the considered optimization problem has shown that gradient-based nonlinear programming solvers perform better if at least matrix Q _p is positive definite, hence we add a small regularization term δ > 0 to the diagonals elements of Q _p so that it becomes diagonal dominant and then positive defined. Let Qp '■= Qp + ^SI ninn ^an d Q ^tne matrix such that x ^T Qx = p ^T Q _p p. If the perturbation term is sufficiently small it will not significantly affect the optimization process. Finally, the second objective function used to account or pressure variability can be defined as μ ₂ ( ) == x ^T Qx.

3 Risk based selection of weights

In minimizing pressure variability across the network, it may be beneficial to consider both the relative criticality and vulnerability of pipes. This is because pipes that are vital in the provision of supply (e.g. pipes connecting important users like hospitals, schools etc. to water sources or have high-consequence of failure) and pipes which are vulnerable to pressure-induced failure (e.g. pipes that are relatively older and pipes with high likelihood of failure) can be prioritized to have lower pressure variations on them (see 60 in Figure 2).

The framework in equations (1) and (2) above can treat the optimization of pressure variability across all pipes of a water distribution network equally by using weights all equal to 1 in (2). However, by appropriately choosing the weights based on pipe criticality (i.e. consequence of failure) and vulnerability (i.e. according to their likelihood of failure), the same proposed framework allows the decision maker to direct the optimization process more towards high risk pipes.

3.1 Critical pipes

The criticality of a pipe can be defined based on the relative impact on the network if that pipe fails. One aspect of the impact relates to the number of customers affected by the interruption of supply in the failure (see 60 in Figure 2). A second aspect relates to the type of customer affected; for example, an interruption of supply for a hospital may be less tolerated compared to an interruption of supply for an industrial user).

For a small network, critical pipes and their impact is measured through hydraulically simulating failures, i.e. using a critical link analysis. Such methods have combinatorial computational complexity, i.e. simulating cascading effects through different combinations of failures takes an infeasibly large amount of time for large and medium size models. In order to address this, preferred embodiments of this disclosure use scalable graph theory tools for automatically identifying most critical and high-risk pipes in the networks

To identify pipes critical to the supply of important users (such as hospitals or schools, or other users important with respect to other social and economic criteria), the preferred embodiment adopts a hydraulic capacity (see, for example, Park, J. I., Lambert, J.H. & Haimes, Y.Y., 1998. Hydraulic power capacity of water distribution networks of deterioration. Water Resources Research, 34(12), pp.3605-3614) weighted K-shortest path algorithm (see, for example, Eppstein, D., 1998. Finding the k shortest paths. SIAM Journal on Optimization, 28(2), pp.652-673). In this approach, by restricting a weighted betweenness computation to paths connecting important nodes (e.g. hospitals) to water sources, the k-shortest paths between each such pair-wise nodes can be computed, and the number of times each pipe appears in such paths can be identified. The paths can be weighted by their hydraulic power capacities in the k-shortest path computations to introduce hydraulic meaning to each pipe's importance. The pipes that appear in all the paths with the highest hydraulic capacity can be deemed critical/important to supply. This may be done for all pairwise combination of sources and important nodes. The pipes that do not appear in the pairwise paths may not be considered critical. A set of weights can also be generated for each pipe based on these computations.

3.2 Vulnerable pipes

The weighting may also take account of the vulnerability of the pipes. For example, special attention can be given to the pipes that are more vulnerable to failure under cumulative pressure variations (cumulative pressure induced stress). By automatically analysing the network geographical information system (GIS) data, pipes can be categorized by age, material type, location, and other factors. By assigning weights based on this, the optimization process can made to prevent large pressure variations on the most vulnerable of the critical pipes while still reducing the average zone pressure.

3.3 Risk Ranking

Using criticality and vulnerability rankings as described above, and overall risk ranking can be defined to establish the weighting for pipes within the network.

4 Solution methods

Having defined the functions to be minimised, the referred embodiment adopts a scalable solution method for the considered non-convex multiobjective mixed-integer nonlinear problem thereby established. In particular, scalarization methods can be used to obtain the Pareto fronts. Techniques that can be adopted to generate the Pareto fronts include the methods mentioned above such as the Normal Boundary Intersection (NBI) method, the Normalized Normal Constraint (NNC) method and the Weighted Sum (WS) method. Each scalarization method addresses the generation of points on the Pareto front by solving a sequence of single-objective M IN LPs.

4.1 Solution of single-objective MINLPs

A possible approach for non-convex MINLP is the application of a branch and bound algorithm to find at least local solutions. However, this would result in infeasible computational time for large scale water distribution networks, since it is necessary to solve a series of MINLPs to generate the Pareto front. Accordingly, in the preferred embodiment, an alternative approach is adopted. Since the integer variables involved in our optimization problem are binaries, they can be reformulated as complementarity constraints, enforcing the variables to take only one of the two complementary values 0 and 1. The feasible set of mathematical programs with complementarity constraints (MPCCs) has a pathological structure that prevents the solution with direct application of standard nonlinear programming solvers. However, MPCCs have been successfully solved using specialized algorithms in other engineering frameworks.

The preferred embodiment may implement relaxation (such as that described in Raghunathan, A.U. & Biegler, L.T., 2005. An Interior Point Method for Mathematical Programs with Complementarity Constraints (MPCCs). SIAM Journal on Optimization, 15(3), pp.720-750; or Hoheisel, T., Kanzow, C. & Schwartz, A., 2011. Theoretical and numerical comparison of relaxation methods for mathematical programs with complementarity constraints. Mathematical Programming, 137(1-2), pp.257-288) and penalization methods (such as those described in: Hu, X.M. & Ralph, D., 2004. Convergence of a penalty method for mathematical programming with complementarity constraints. Journal of Optimization Theory and Applications, 123(2), pp.365-390; Leyffer, S., Lopez-Calva, G. & Nocedal, J., 2006. Interior Methods for Mathematical Programs with Complementarity Constraints. SIAM Journal on Optimization, 17(1), pp.52-77; or Ralph, D. & Wright, S.J., 2004. Some properties of regularization and penalization schemes for MPECs. Optimization Methods and Software, 19(5), pp.527-556.) for the solution of the MPCCs related to the considered co-design optimization problem. The former relies on the solution of a sequence of relaxed problems with 'regular' approximations for the complementarity constraints; these approximate feasibility sets are described using a relaxation parameter. The relaxed problems can be solved with standard NLP solvers. The relaxed feasible sets converge to the pathological feasibility set of the MPCC, hence the sequence of solutions will converge to a stationary solution of the original problem.

Alternatively to relaxation methods, the implementation of penalization methods may be adopted, where the optimization is performed in a feasible set obtained by dropping the complementarity constraints, which are then embedded into the objective function via a function that penalizes complementarity violations. A sequence of penalty problems can then be solved using standard NLP solvers and the solutions will converge to a stationary point of the original MPCC problem (Hu, X.M. & Ralph, D., 2004. Convergence of a penalty method for mathematical programming with complementarity constraints. Journal of Optimization Theory and Applications, 123(2), pp.365-390). Moreover, sparse techniques can be used within the NLP solvers to improve the performance of relaxation and penalization methods since the large scale mixed integer nonlinear programs that arise in the framework of optimal co-design for water distribution networks have sparse constraints. The Jacobians and Hessians of relaxed and penalized NLPs retain this sparsity, which can be taken advantage of.

4.2 Overview of the solution method

With reference to Figure 3, the overall solution method can be understood as follows.

At step 7a, a multi-start technique is adopted. This allows the process to obtain, at step 7b, multiple different individual local minima (anchor points) for the two objectives (AZP and PVZ). In the next step 7c, a parallelization process is operated on multicore processors of multiple different runs of the scalarization methods, each of which is initialized by the plurality of anchor points generated at step 7b.

Step 7d then operates in parallel and allows the application of the user preferred or all scalarization methods mentioned above - NBI, NNC, WS. Subsequently, at step 7e the solution of the single-objective MINLPs generated by each instance of the scalarization method is computed. At step 7f, the Pareto front is composed from the identified instances. At step 7g a Pareto filter is applied to the different Pareto fronts generated by the previous steps in order to eliminate dominated points (step 7h); as a result at step 7i a robust non-dominated front can be generated.

The performance of the scalarization methods can be improved by using the information obtained from the solution of a single-objective problem to initialize the subsequent iteration. Relaxation and penalty methods require the solution of a series of sparse nonlinear programs (NLPs) which can be efficiently solved by a plurality of state-of-the-art sparse NLP solvers, implementing active-set or interior point algorithms.

Variations and modifications will be apparent to the skilled person. Such variations and modifications may involve equivalent and other features which are already known and which may be used instead of, or in addition to, features described herein. Features that are described in the context of separate embodiments may be provided in combination in a single embodiment. Conversely, features which are described in the context of a single embodiment may also be provided separately or in any suitable sub-combination. It should be noted that the term "comprising" does not exclude other elements or steps, the term "a" or "an" does not exclude a plurality, a single feature may fulfil the functions of several features recited in the claims and reference signs in the claims shall not be construed as limiting the scope of the claims. It should also be noted that the Figures are not necessarily to scale; emphasis instead generally being placed upon illustrating the principles of the present invention.