Economic Model Predictive Control with Nonlinear Constraint Relaxation for the Operational Management of Water Distribution Networks

: This paper presents the application of an economic model predictive control (MPC) for the operational management of water distribution networks (WDNs) with periodic operation and nonlinear constraint relaxation. In addition to minimizing operational costs, the proposed approach aims to reduce the computational load and to improve the implementation efﬁciency associated with the nonlinear nature of the MPC problem. The behavior of the WDN is characterized by a set of difference-algebraic equations, where the relation of hydraulic pressure/head and ﬂow in interconnected pipes is nonlinear. Speciﬁcally, the considered WDN model includes two categories of nonlinear algebraic equations for unidirectional and bidirectional ﬂows in pipes, respectively. In this paper, we propose an iterative algorithm to relax these nonlinear algebraic equations into a set of linear inequality constraints that will be implemented in the economic MPC design, which improves the implementation efﬁciency and meanwhile optimizes the economic performance. Finally, the proposed strategy is applied to a well-known benchmark of the Richmond WDN. The closed-loop simulation results are shown and the proposed strategy is also compared with a nonlinear economic MPC using several key performance indexes.


Introduction
Water distribution networks (WDNs) are complex and large-scale interconnected systems, which are designed to supply water to consumers in cities [1]. They are composed of a number of storage tanks located across urban cities in different elevations, pumping stations including parallel booster pumps, pressure reducing valves, water demand nodes as well as interconnected pipes. The network configurations may have a complex mesh layout to achieve the necessary services with a convenient geographical and topological distribution [2,3]. Moreover, these networks have been built incrementally throughout the urban development and increasing population, and contain pipes of different materials, lengths and diameters. By revising the literature, a short overview of water systems and some current challenging water issues can be found in [4].
The dynamical and static behaviors of a WDN can be mathematically described by a set of nonlinear difference-algebraic equations (DAEs) in discrete-time [5][6][7]. The system dynamics are described by difference equations while static relations involving flows and pressures are characterized by algebraic equations. Following the modeling approach presented in [6], difference equations are obtained from the mass balance at the storage tanks and algebraic equations come from (i) linear algebraic equations obtained from the mass balance at the non-storage nodes; (ii) nonlinear algebraic equations describing the relation of hydraulic pressure/head and flow in interconnected pipes. Depending on the pipe usage, two categories of nonlinear algebraic equations are considered: (i) (in unidirectional pipes) since the water is transported only in one direction, the water flow can be always positive by choosing a direction; (ii) (in bidirectional pipes) since the direction may be reversed, the water flow can be positive and negative.
Model predictive control (MPC), as one of the most popular optimal control strategies, has attracted a lot of attention in a large number of industrial applications [8,9]. Recently, economic MPC appears [10][11][12][13][14][15], where the optimal control actions are computed by optimizing a purely economic cost function that directly measures the economic performance of the considered system. In contrast with the tracking MPC strategy designed to track a given trajectory or reference, general economic cost functions are used and these may not necessarily be positive-definite with respect to the optimal steady states.
For the operational management of WDNs, different optimal control strategies have been investigated (see, e.g., [16][17][18][19][20][21]). For achieving multiple objectives, MPC strategies, and in particular economic MPC, have already been applied to WDNs (see, e.g., [5,6,[22][23][24][25][26][27]). The main challenge of applying MPC strategy to the WDN is to solve a nonlinear optimization problem taking into account the full system dynamics and a large amount of decision variables within a sampling time. Although for WDNs the sampling time is usually chosen to be one hour, suboptimal solutions may be obtained after a huge number of computational iterations.
In terms of the operational behavior of the WDN, some endogenous and exogenous variables usually present daily patterns. For instance, the urban water demand presents a daily cycle and the periodically time-varying electricity price is usually selected as an exogenous signal for evaluating the energy cost of pumping in a WDN. The electricity tariff also presents a periodicity of 24 h. Thus, this naturally induces a periodic operation of the closed-loop system; that is, the optimal system trajectories follow a cyclic pattern. A single-layer economic MPC for periodic operation was proposed in [28] for a WDN case study.
The main contribution of this paper is to present an application of economic MPC strategy with nonlinear constraint relaxation for the operational management of WDNs. In this paper, the control-oriented model of WDNs uses a set of difference and algebraic equations described in [6]. Besides reducing operational costs, the proposed strategy is able to deal with nonlinear algebraic equations in order to reduce the computational load and improve the implementation efficiency associated. To this aim, firstly, we propose an iterative algorithm to relax nonlinear algebraic equations into a set of linear inequality constraints. Then, these obtained linear inequality constraints are integrated into the economic MPC design for periodic operation. Finally, we apply the proposed strategy to a real benchmark, namely the Richmond WDN, which is a medium-size network in the UK. Furthermore, we compare the proposed strategy with a nonlinear economic MPC and we also assess the performance with relaxed constraints using defined key performance indicators (KPIs).
The remainder of the paper is structured as follows. Economic MPC of WDNs is presented in Section 2. An iterative algorithm for nonlinear constraint relaxation is introduced in Section 3 as well as the economic MPC strategy with relaxed constraints for periodic operation. The Richmond WDN case study is chosen to test the proposed economic MPC strategy with nonlinear constraint relaxation in Section 4. Finally, conclusions are addressed in Section 5.

Control-Oriented Model of WDNs
Consider that a WDN can be formulated as a discrete-time DAE model (see [6] for more details) for control purposes where x ∈ R n x denotes the difference state vector of hydraulic heads (the sum of elevation and water level) for water storage tanks, u ∈ R n u denotes the control input vector of water flows through actuators (pumps and valves), v ∈ R n v denotes the input vector of non-manipulated water flows through pipes, d ∈ R n d denotes the known input vector of water demands, z ∈ R n z denotes the algebraic state vector of hydraulic heads of the non-storage nodes, and k ∈ N is the sampling time instant. Moreover, x ∈ R m n ×n x and P z ∈ R m n ×n z are the system matrices. ψ(v) ∈ R n e is a collection of nonlinear functions that describes the static relation of the hydraulic head and the water flow though interconnected pipes in (1c). We assume that the state variables x(k) are fully measurable, ∀k ∈ N and the water demand d(k) can be known in a short-term prediction horizon using a certain demand forecasting method, for instance [29]. Besides, the water demand d(k) also presents a periodic behavior, that is d(k) = d(k + T), ∀k ∈ N, where T is a period. In the WDN, the system variables in (1) are constrained according to the physical limitations as follows: where x and u are lower bounds of x and u while x and u are upper bounds of x and u.

Economic MPC of WDNs
In this section, we present an economic MPC for the optimal operational management of WDNs. Firstly, the constraints and economic cost functions are defined based on the control objectives for the WDN management. Then, we formulate the optimization problem to implement the proposed economic MPC. Considering the periodicity in the WDN operation, the MPC prediction horizon H p is chosen to be H p = T.
The main control objectives for the management of a WDN, considered in this paper, are summarized as follows: • Economic: To minimize the operational costs while supplying enough water to satisfy all the demands with suitable pressure; • Safety: To reserve suitable volumes of water in storage tanks for satisfying all the underlying uncertain demands.
To achieve these objectives, system constraints and the individual cost function with respect to each objective are formulated in the following. Firstly, we discuss the constraints: In terms of a MPC strategy, the prediction model is chosen as (1a) and (1b), and the nonlinear algebraic equation (1c). Besides, system variables are constrained by (2a) and (2b). To guarantee the safety objective, the minimal pressure at the demand sectors is required, which can be achieved by setting the following constraint on the hydraulic head z as where z guarantees the minimal required pressure for the nodes inside the WDN. Secondly, we discuss the cost function settings for the management of WDNs. We would like to penalize the economic cost with respect to the periodic price signal sequence p along the MPC prediction horizon to achieve the economic performance. Specifically, we define the economic cost function as where p is modeling a periodic variation of the electricity price with the period T as Besides, the water levels at storage tanks are penalized by which means that we would like to reserve the minimal required water volume inside the storage tanks over the prediction horizon. In general, the economic cost function along the MPC prediction horizon is defined as where λ 1 and λ 2 are prioritization weights. The general optimization problem for the operational management of WDNs is formulated as follows: subject to where i = 0, . . . , T − 1 and J T (x, u) denotes the economic cost function for management of WDNs, andx(k) is the measured state at time instant k and mod(a, b) is the modulo operation of two numbers a and b.
As formulated in (7), this is a nonlinear optimization problem due to the nonlinear algebraic equation (1c). The computational cost of solving a nonlinear optimization problem is high and leads to a suboptimal solution. To avoid solving a nonlinear optimization problem, a suitable algorithm to deal with the nonlinear algebraic constraints of WDNs is proposed in the next section.
Besides, taking into account the periodicity presented in the WDN (induced by daily periodic water demands and electricity prices), we expect to obtain an optimal periodic (or quasi-periodic) closed-loop trajectory with the designed economic MPC controller for WDNs.

Economic Model Predictive Control of WDNs with Nonlinear Constraint Relaxation
In this section, we propose an iterative algorithm to relax the nonlinear algebraic equation (1c) to obtain linear inequality constraints for bounding (1c) that will be plugged in the economic MPC design. For the relaxation of nonlinear constraints (1c), two cases are considered as follows.
Let us denote v i , i = 1, . . . , n v as the water flow for the pipe i and ∆h i , i = 1, . . . , n v as the i-th row of P x x + P z z. According to [30], Therefore, the head-flow relation in the nonlinear algebraic equation (1c) can be explicitly written as where α i ∈ R + is the pipe resistance coefficient for the i-th nonlinear equation due to friction with the pipe, and β is flow exponent that depends on the particular approximation, such as in Hazen-Williams, Darcy-Weisbach and Chezy-Manning formulas but, in all cases β > 1 according to ([30], Table 3.1).
The interconnected pipes in WDNs may be unidirectional or bidirectional. For the unidirectional pipe with a chosen positive direction, (8) becomes where v i denotes the upper bound of the i-th flow.
The goal of dealing with these nonlinear algebraic equations in (8) is to relax them obtaining a set of linear inequality constraints using an iterative over-bounding algorithm. Note that finding these linear constraints with a proper constraint relaxation method in this work is different than the traditional linearization method with a chosen operating point.

Nonlinear Constraint Relaxation for Unidirectional Pipes
We first discuss the relaxation for unidirectional pipes. By choosing a positive direction of the Figure 1, the objective is to find a set of upper and lower linear bounds for over-bounding this term (shown in blue solid line).
The nonlinear algebraic equation (9) is equivalent to the satisfaction of the following inequalities: Therefore, the inequality (11a) can be relaxed considering (10) as On the other hand, from the convex nature of v β i , we have that every linearization constitutes a lower bound (dashed dotted lines in Figure 1). The constraint (11b) can be approximated by considering N a sampled operating points v i,j for j = 1, 2, . . . , N a as in which parameters a j and b j are given by In general, for a unidirectional pipe, the nonlinear algebraic equation (9) can be relaxed by using N a + 1 inequality constraints as presented in (12) and (13). Figure 1 shows graphically the obtained relaxation. As a potential improvement, this relaxation can be refined iteratively. The iterative algorithm of nonlinear constraint relaxation can be improved by adding a penalty term in order to refine the region of v i . As shown in Figure 2, the upper bound can be moved by a scalar (12) can be replaced by where a small positive τ i can be found in the MPC optimization loop. Hence, the cost function for the scalar τ i (j) varying in the MPC prediction horizon H p = T can be penalized as where λ e (j), j = 1, 2, . . . , T is a weight that can be set as a forgetting (monotonically decreasing) factor along the MPC prediction horizon T.

Nonlinear Constraint Relaxation for Bidirectional Pipes
As shown in Figure 3, the goal is to relax the nonlinear algebraic equation for bidirectional pipes also by linear inequality constraints. As in the unidirectional case, the nonlinear algebraic equation (8) is equivalent to From (17a) and (17b), we can see that these two inequality constraints are not convex along v i ≤ v i ≤ v i . In order to obtain a convex relaxation for (17a), we consider lower bounds for |v i | β−1 with v i ≤ v i ≤ v i in the following form: where a l j and b l j for j = 1, . . . , N b + 1 are two scalars. With a given a j , the condition for the parameter b j should be satisfied: and let us consider the right side of the previous inequality: where We now summarize the way to find a l j and b l j for j = 1, . . . , where v l, i can be obtained by satisfying the following condition: Denote for j = 2, . . . , N b + 1. Therefore, the parameter b l j can be obtained by Furthermore, the upper and lower bounds are symmetric as shown in Figure 3. Therefore, we can find the upper bounds in a similar way. Let us consider upper bounds of (17b) in the following form: Because of symmetry, a r 1 can be determined by where v r, i can be obtained by satisfying the following condition: Denote for j = 2, . . . , N b + 1. Therefore, the parameter b r j can be computed as As a result, (8) for bidirectional pipes can be relaxed as 2N b + 2 linear inequalities. From (17a) and (17b), we obtain the relaxed linear inequality constraints: for j = 1, . . . , N b + 1.

Economic MPC with Nonlinear Constraint Relaxation
From the above results, we can obtain the relaxed linear inequality constraints in the MPC prediction horizon H p = T as follows: for j = 1, . . . , T. Taking into account the proposed iterative algorithm for the nonlinear constraint relaxation, the nonlinear algebraic constraint (1c) can be replaced by (32a) and (32b) along the MPC prediction horizon. Besides, since the optimal states are expected to be periodic, the current measured statex(k) is inserted into the shifted position by introducing the modulo operator mod(·). We now formulate the optimization problem for implementing the economic MPC with nonlinear constraint relaxation as follows minimize x(0),...,x(T) u(0),...,u(T−1) subject to with i = 0, . . . , T − 1.
Note that the optimization problem (33) is solved with a given prediction horizon T. Let us denote the optimal solution of the optimization problem (33) as u * (j). Based on the receding horizon strategy, the optimal control action u(k) at time step k is chosen as

Application: the Richmond Water Distribution Network
In this section, we apply the proposed relaxation algorithm and control strategy to a realistic benchmark network called the Richmond WDN (This benchmark is available from the link: http: //emps.exeter.ac.uk/engineering/research/cws/resources/benchmarks), which is a medium-size network and the control-oriented model used in EPANET is formulated in the form of (1). In particular, nonlinear algebraic equations are included.

System Description
The topology and layout of the Richmond WDN is shown in Figure 4. The Richmond WDN has 6 water storage tanks, 7 booster pumps and 11 water demand sectors. Besides, there are 41 non-storage nodes and 41 pressurized pipes connected in this network. Following the modeling methodology presented in [6], the control-oriented model of the Richmond WDN is in the DAE form (1). By using the mass balance at storage tanks, we obtain the system dynamics in (1a). The demand pattern is also given for a 24-h period, that is T = 24. By also using the mass balance at non-storage nodes, we obtain the linear algebraic equation in (1b). Besides, we use the Chezy-Manning head-flow formula to obtain (1c) as follows [30]: where z i and z j correspond to the hydraulic heads at any two adjacent nodes, and v i,j is the corresponding water flow. The parameter R i,j in the Chezy-Manny formula is given by where L i,j , D i,j and C i,j are the length, diameter and roughness coefficient of the corresponding pipe, respectively. L i,j and D i,j are given in the EPANET model of the Richmond WDN. The values of C i,j in Chezy-Manny formula are provided for each pipe in the Richmond network in Table 1. As shown in Figure 4, this WDN has two bidirectional pipes in (8) and 39 unidirectional pipes in (9). In addition to the economic cost function defined in (6), for the relaxed linear constraints (32a) and (32b) with the setting in (15), a penalty term λ e (j) is set to be a forgetting factor as where denotes the relaxed step and λ e is the initial value of this weight. Therefore, the total MPC cost function is set to bē where τ denotes a slack variable for all the constraints in (15). In this simulation, we select the weights as λ 1 = 10, λ 2 = 1, λ e = 0.1 and = 0.01 obtained from the tuning procedure presented in [31]. For the Richmond network, the period T is considered to be T = 24h with the sampling time ∆t = 1 h because of the periodicities of the water demand and electricity price considering the variations in the daily tariff. Hence, the prediction horizon of the proposed economic MPC strategy is also chosen to be T = 24 h. The minimal pressure at all the demand sectors is set to be 10 m. Furthermore, for the implementation of the proposed nonlinear constraint relaxation, we choose N a = N b = 10. Therefore, there are 11 relaxed constraints replacing (9) and 22 relaxed constraints replacing (8) for each pipe.
To assess the performance of the proposed iterative algorithm in the economic MPC, a nonlinear periodic economic MPC strategy can be implemented by solving the following optimization problem: subject to x(j) =x(k), j = mod (k, T) .
Besides, the optimal periodic steady trajectory corresponding to this nonlinear periodic economic MPC can be found. Following [28], we present a finite-horizon optimization problem with the complete model of WDNs in (1) to find the optimal periodic steady trajectory, known as the planner optimization problem. This optimization problem yields the same solution if the time frame to be considered is any period, that is i ∈ [k, k + T], ∀k ∈ N.
subject to Note that the planner optimization problem (40) is set without any initial state. Hence, we only need to solve it once to obtain the optimal solutions as the optimal periodic steady trajectory.
The simulations have been carried out with the MATLAB R2015a and the EPANET simulator [30] for seven days (168 h) in a PC of Intel i7-5500U CPU and 12GB RAM. The linear optimization problems are solved using the Yalmip toolbox [32] and the Mosek solver [33]. The nonlinear optimization problems are solved using the nonlinear programming through the Yalmip toolbox and the IPOPT solver available in the OPTI toolbox [34]. The Richmond network is given in the EPANET simulator as the simulation model.

Key Performance Indicators
To compare the performance of the proposed economic MPC with the nonlinear economic MPC, we define the following KPIs: where KPI E is the economic KPI that measures the operational costs at each time step. KPI S is the safety KPI that computes the average differences of the water storage that are lower than safety hydraulic head x s i given in Table 2. KPI M is the measurement KPI that represents the additional water reserved in storage tanks. Based on the original benchmark available online, all the tanks are cylindrical and the relationship between water level and stored volume is considered to be linear and constant.
On the other hand, with the optimal solutions of the optimization problem (33), we would like to check whether all the nonlinear algebraic equations in (1c) are satisfied. To assess the relaxation algorithm for 40 nonlinear algebraic equations in the Richmond WDN, the error measurements for (1c) including mean squared error (MSE), mean absolute error (MAE) and symmetric mean absolute percentage error (SMAPE) are introduced as follows: where P j x , P j z and ψ j (·) denote the j-th row of P x , P z and ψ(·), respectively. n e denotes the total number of nonlinear algebraic equations. In terms of MSE and MAE, they represent the violation of nonlinear algebraic equations. SMAPE is an indicator based on percentage errors.

Simulation Results
For the notation simplicity, we denote the simulation results of applying the proposed economic MPC with nonlinear constraint relaxation as EMPC-NCR, while the comparison results with the solutions of optimization problems (39) and (40) are denoted by EMPC-Planner and NEMPC, respectively. The closed-loop simulation results of system states and control inputs are shown in Figures 5-8. In Figures 5 and 6, the state trajectories obtained from applying the proposed EMPC-NCT are in solid lines with circles. Due to the convexity, the steady states can be obtained from the solution of the optimization problem (40) shown in dashed line. As a comparison, the state trajectories of NEMPC are also shown in solid lines with cross marks. From these results, we can see that the closed-loop trajectories obtained using the EMPC-NCR and NEMPC strategies are similar to those of the optimal planner trajectories (both states and control inputs). The NEMPC results are smoother and closer to the planner trajectories since a more accurate nonlinear model is used in the NEMPC optimization problem. Similarly, in terms of control inputs, as shown in Figures 7 and 8, three trajectories of EMPC-NCR, NEMPC and EMPC-Planner are plotted. The input trajectories of EMPC-NCR are approaching the ones of EMPC-Planner.
To assess the performance of different control strategies, the comparison is also provided based on the defined KPIs. The computation results using the defined KPIs are shown in Table 3. In general, the performances of both MPC strategies are similar. Specifically, from the KPI E results, the pure economic cost of EMPC-NCR is slightly cheaper than the one of NEMPC. According to KPI S and KPI M results, small differences between the reserved water in the storage tanks can be seen for both MPC strategies. This is because in the EMPC-NCR we use the pressure constraint on the variable z to guarantee the safety objective, which implies the water levels in the storage tanks should be greater than some certain values.   The results of error indicators for the EMPC-NCR and NEMPC strategies are shown in Figure 9. Through the MSE and MAE results, it is obvious that the result of NEMPC is similar to the one of EMPC-NCR, although none of them are identically equal to zero. This is because the tolerance of the nonlinear solver is chosen as 10 −5 . The SMAPE of NEMPC is smaller and closer to zero than the one of EMPC-NCR, which means that nonlinear algebraic equations are satisfied by NEMPC better than EMPC-NCR since the nonlinear programming technique is able to solve hard constraints. However, the EMPC-NCR strategy is able to produce a similar performance according to three error measurement results.
For the comparison of simulation time in a scenario of 168 h, it takes 62.86 min for NEMPC while 1.43 min for EMPC-NCR. Hence, the EMPC-NCR strategy has a significant improvement in the reduction of computational load and meanwhile based on the above comparison result, the performance of the EMPC-NCR strategy is similar to the NEMPC strategy. This reduction in the computation time would be more relevant in larger networks.

Conclusions
In this paper, we have presented the application of an economic MPC strategy with nonlinear algebraic constraint relaxation for the operational management of WDNs. To reduce the computational load and to improve the implementation efficiency, we propose an iterative algorithm to handle nonlinear algebraic equations in the control-oriented model of WDNs. When the size of the WDN increases, the proposed relaxation algorithm can significantly reduce the computational load and transform the nonlinear optimization problem into a linear one. From the application results of the Richmond WDN, the proposed strategy is tested and its performance is similar to a corresponding nonlinear economic MPC. From the comparison of the computation time, the proposed strategy is significantly faster than using a nonlinear economic MPC. Hence, the proposed strategy could be applied to a large-scale WDN.