Next Article in Journal
Optimal Design of a Band Pass Filter and an Algorithm for Series Arc Detection
Previous Article in Journal
A Hydropower Biological Evaluation Toolset (HBET) for Characterizing Hydraulic Conditions and Impacts of Hydro-Structures on Fish
Previous Article in Special Issue
Multi-Model Prediction for Demand Forecast in Water Distribution Networks
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Advanced Control Systems (ACS) Research Group at Institut de Robòtica i Informàtica Industrial (IRI), CSIC-UPC, Universitat Politècnica de Catalunya-BarcelonaTech (UPC), C/. Llorens i Artigas 4-6, 08028 Barcelona, Spain
2
Departamento de Ingeniería de Sistemas y Automática, Universidad de Sevilla, Avenida de los Descubrimientos, S/N, 41092 Sevilla, Spain
3
Cetaqua, Water Technology Centre, Ctra. d’Esplugues 75, Cornellà de Llobregat, 08940 Barcelona, Spain
*
Author to whom correspondence should be addressed.
Energies 2018, 11(4), 991; https://doi.org/10.3390/en11040991
Submission received: 28 February 2018 / Revised: 13 April 2018 / Accepted: 16 April 2018 / Published: 19 April 2018
(This article belongs to the Special Issue Smart Water Networks in Urban Environments)

Abstract

:
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 efficiency 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 flow in interconnected pipes is nonlinear. Specifically, the considered WDN model includes two categories of nonlinear algebraic equations for unidirectional and bidirectional flows 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 efficiency 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.

1. 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.

2. Economic Model Predictive Control of WDNs

2.1. 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
x ( k + 1 ) = A x ( k ) + B u u ( k ) + B v v ( k ) + B d d ( k ) ,
0 = E u u ( k ) + E v v ( k ) + E d d ( k ) ,
0 = P x x ( k ) + P z z ( k ) + ψ v ( k ) ,
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, A R n x × n x , B u R n x × n u , B v R n x × n v , B d R n x × n d , E u R m l × n u , E v R m l × n v , E d R m l × n d , P 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:
x ̲ x ( k ) x ¯ , k N ,
u ̲ u ( k ) u ¯ , k N ,
where x ̲ and u ̲ are lower bounds of x and u while x ¯ and u ¯ are upper bounds of x and u.

2.2. 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
z ( k ) z ̲ , k N ,
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
e u ( k ) , p i = p i T u ( k ) , i = mod ( k , T ) .
where p is modeling a periodic variation of the electricity price with the period T as
p = p 1 , , p T .
Besides, the water levels at storage tanks are penalized by
x x ( k ) = x ( k ) 1 ,
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
J T ( x , u , p ) = i = 0 T 1 λ 1 e u ( i ) , p i + λ 2 x x ( i ) ,
where λ 1 and λ 2 are prioritization weights. The general optimization problem for the operational management of WDNs is formulated as follows:
minimize x ( 0 ) , , x ( T ) u ( 0 ) , , u ( T 1 ) J T x , u , p ,
subject to
1 1 x ( i + 1 ) = A x ( i ) + B u u ( i ) + B v v ( i ) + B d d ( i ) ,
0 = E u u ( i ) + E v v ( i ) + E d d ( i ) ,
0 = P x x ( i ) + P z z ( i ) + ψ v ( i ) ,
x ̲ x ( i ) x ¯ ,
u ̲ u ( i ) u ¯ ,
x ( 0 ) = x ( T ) ,
x ( j ) = x ^ ( k ) , j = mod k , T .
where i = 0 , , T 1 and J T x , u denotes the economic cost function for management of WDNs, and x ^ ( 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.

3. 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], ψ v i = α i v i v i β 1 . Therefore, the head-flow relation in the nonlinear algebraic equation (1c) can be explicitly written as
0 = α i v i v i β 1 + Δ h i ,
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
0 = α i v i β + Δ h i ,
with
0 v i v ¯ i ,
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.

3.1. Nonlinear Constraint Relaxation for Unidirectional Pipes

We first discuss the relaxation for unidirectional pipes. By choosing a positive direction of the flow v i , the nonlinear term v i v i β 1 becomes v i v i β 1 = v i β with v i 0 . As shown in 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:
α i v i β + Δ h i 0 ,
α i v i β + Δ h i 0 ,
in which v i β is a convex function due to β > 1 . From 0 v i v ¯ i , we have that v i β v ¯ i β 1 v i . Therefore, the inequality (11a) can be relaxed considering (10) as
α i v ¯ i β 1 v i + Δ h i 0 .
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
0 α i v i β + Δ h i α i ( a j v i + b j ) + Δ h i ,
in which parameters a j and b j are given by
a j = β v i , j β 1 ,
b j = ( 1 β ) v i , j β .
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 τ i > 0 . Therefore, the region of v i is refined to be v i a , v i b 0 , v ¯ i .
Considering a slack decision variable τ i , (12) can be replaced by
α i v ¯ i β 1 v i + Δ h i τ i 0 ,
τ i 0 ,
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
e ( τ i ( j ) ) λ e ( j ) τ i ( j ) ,
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.

3.2. 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
α i v i v i β 1 + Δ h i 0 ,
α i v i v i β 1 + Δ h i 0 .
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:
a j l v i + b j l v i v i β 1 , j = 1 , , N b + 1 ,
where a j l and b j l for j = 1 , , N b + 1 are two scalars. With a given a j , the condition for the parameter b j should be satisfied:
b j l min v ̲ i v i v ¯ i ( v i v i β 1 a j l v i ) ,
and let us consider the right side of the previous inequality:
M min v ̲ i v i v ¯ i ( v i v i β 1 a j l v i )
= min M + , M ,
where
M + = min 0 v i v ¯ i ( v i β a j l v i ) ,
M = min v ̲ i v i < 0 ( v i ( v i ) β 1 a j l v i ) .
We now summarize the way to find a j l and b j l for j = 1 , , N b + 1 . The minimum a 1 l along v ̲ i v i v ¯ i can be determined by
a 1 l = β v i l , β 1 ,
where v i l , can be obtained by satisfying the following condition:
β v i l , β 1 = v i l , β + v ̲ i β v i l , v ̲ i .
Denote v i , 1 l = v i l , . By choosing N b values of v i , j l in the interval v i l , v i , j l v ¯ i , we obtain a j l = β v i , j l β 1 for j = 2 , , N b + 1 . Therefore, the parameter b j l can be obtained by
b j l = v i , j l β a j l v i , j l , j = 1 , , N b + 1 .
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:
a j r v i + b j r v i v i β 1 , j = 1 , , N b + 1 .
Because of symmetry, a 1 r can be determined by
a 1 r = β v i r , β 1 ,
where v i r , can be obtained by satisfying the following condition:
β v i r , β 1 = v ¯ i β + v i r , β v ¯ i v i r , .
Denote v i , 1 r = v i r , . Similarly, by choosing N b values of v i , j r in the interval v ̲ i v i v i r , , we obtain a j r = β v i , j r β 1 for j = 2 , , N b + 1 . Therefore, the parameter b j r can be computed as
b j r = v i , j r ( v i , j r ) β 1 + a j r v i , j r , j = 1 , , N b + 1 .
As a result, (8) for bidirectional pipes can be relaxed as 2 N b + 2 linear inequalities. From (17a) and (17b), we obtain the relaxed linear inequality constraints:
0 α i v i v i β 1 + Δ h i α i a j l v i + b j l + Δ h i ,
0 α i v i v i β 1 + Δ h i α i a j r v i + b j r + Δ h i ,
for j = 1 , , N b + 1 .

3.3. 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:
P ˜ x l ( j ) x ( k + j ) + P ˜ z l ( j ) z ( k + j ) + P ˜ v l ( j ) v ( k + j ) + P ˜ b l ( j ) 0 ,
P ˜ x r ( j ) x ( k + j ) + P ˜ z r ( j ) z ( k + j ) + P ˜ v r ( j ) v ( k + j ) + P ˜ b r ( j ) 0 ,
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 state x ^ ( 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 ) J T ( x , u , p ) ,
subject to
1 1 x ( i + 1 ) = A x ( i ) + B u u ( i ) + B v v ( i ) + B d d ( i ) ,
0 = E u u ( i ) + E v v ( i ) + E d d ( i ) ,
P ˜ x l ( i ) x ( i ) + P ˜ z l ( i ) z ( i ) + P ˜ v l ( i ) v ( i ) + P ˜ b l ( i ) 0 ,
P ˜ x r ( i ) x ( i ) + P ˜ z r ( i ) z ( i ) + P ˜ v r ( i ) v ( i ) + P ˜ b r ( i ) 0 ,
x ̲ x ( i ) x ¯ ,
u ̲ u ( i ) u ¯ ,
x ( 0 ) = x ( T ) ,
x ( j ) = x ^ ( k ) , j = mod k , T ,
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
u ( k ) = u * ( j ) , j = mod k , T .

4. 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.

4.1. 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]:
z i z j = R i , j v i , j v i , j ,
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
R i , j = 10.29 L i , j C i , j 2 D i , j 5.33 ,
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
λ e ( j ) = max λ e ( j 1 ) ϵ , 0 ,
λ e ( 0 ) = λ e ,
where ϵ denotes the relaxed step and λ e is the initial value of this weight.
Therefore, the total MPC cost function is set to be
J ¯ T ( x , u , p ) = i = 0 T 1 λ 1 e u ( i ) , p i + λ 2 x x ( i ) + λ e ( i ) τ ( i ) ,
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 = 24 h 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:
minimize u ( 0 ) , , u ( T 1 ) J ¯ T ( x , u , p ) ,
subject to
1 1 x ( i + 1 ) = A x ( i ) + B u u ( i ) + B v v ( i ) + B d d ( i ) ,
0 = E u u ( i ) + E v v ( i ) + E d d ( i ) ,
0 = P x x ( i ) + P z z ( i ) + ψ v ( i ) ,
x ̲ x ( i ) x ¯ ,
u ̲ u ( i ) u ¯ ,
x ( 0 ) = x ( T ) ,
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 .
minimize x ( 0 ) , , x ( T ) , u ( 0 ) , , u ( T 1 ) J ¯ T ( x , u , p ) ,
subject to
1 1 x ( i + 1 ) = A x ( i ) + B u u ( i ) + B v v ( i ) + B d d ( i ) ,
0 = E u u ( i ) + E v v ( i ) + E d d ( i ) ,
0 = P x x ( i ) + P z z ( i ) + ψ v ( i ) ,
x ̲ x ( i ) x ¯ ,
u ̲ u ( i ) u ¯ ,
x ( 0 ) = x ( T ) .
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.

4.2. Key Performance Indicators

To compare the performance of the proposed economic MPC with the nonlinear economic MPC, we define the following KPIs:
K P I E 1 n s k = 1 n s p j T u ( k ) , j = mod ( k , T ) ,
K P I S 1 n s k = 1 n s i = 1 n x max 0 , x i s x i ( k ) ,
K P I M 1 n s k = 1 n s i = 1 n x x i ( k ) x i s ,
where K P I E is the economic KPI that measures the operational costs at each time step. K P I S is the safety KPI that computes the average differences of the water storage that are lower than safety hydraulic head x i s given in Table 2. K P I 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:
M S E ( k ) = 1 n e j = 1 n e ( P x j x ( k ) + P z j z ( k ) + ψ j v ( k ) ) 2 ,
M A E ( k ) = 1 n e j = 1 n e P x j x ( k ) + P z j z ( k ) + ψ j v ( k ) ,
S M A P E ( k ) = 100 n e j = 1 n e P x j x ( k ) + P z j z ( k ) + ψ j v ( k ) P x j x ( k ) + P z j z ( k ) + ψ j v ( k ) ,
where P x j , P z j 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.

4.3. 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 Figure 5, Figure 6, Figure 7 and Figure 8. In Figure 5 and Figure 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 Figure 7 and Figure 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 K P I E results, the pure economic cost of EMPC-NCR is slightly cheaper than the one of NEMPC. According to K P I S and K P I 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.

5. 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.

Acknowledgments

This work was partially funded by the Spanish State Research Agency (AEI) and the European Regional Development Fund (ERFD) through the project DEOCS (ref. DPI2016-76493), the María de Maeztu Seal of Excellence to IRI (ref. MDM-2016-0656) and the FPI grant (ref. BES-2014-068319).

Author Contributions

All the authors have participated in the discussions and preparation of this work.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Mays, L. Water Distribution Systems Handbook; McGraw-Hill: New York, NY, USA, 2000. [Google Scholar]
  2. Yazdani, A.; Jeffrey, P. A complex network approach to robustness and vulnerability of spatially organized water distribution networks. Phys. Soc. 2010, 15, 1–18. [Google Scholar]
  3. Di Nardo, A.; Natale, M.D.; Giudicianni, C.; Musmarra, D.; Varela, J.R.; Santonastaso, G.; Simone, A.; Tzatchkov, V. Redundancy Features of Water Distribution Systems. In Proceedings of the International Conference on Water Distribution Systems, Cartagena, DC, USA, 24–28 July 2016; pp. 412–419. [Google Scholar]
  4. Brown, C.; Lund, J.; Cai, X.; Reed, P.; Zagona, E.; Ostfeld, A.; Hall, J.; Characklis, G.; Yu, W.; Brekke, L. The future of water resources systems analysis: Toward a scientific framework for sustainable water management. Water Resour. Res. 2015, 51, 6110–6124. [Google Scholar] [CrossRef]
  5. Puig, V.; Ocampo-Martinez, C.; Perez, R.; Cembrano, G.; Quevedo, J.; Escobet, T. Real-Time Monitoring and Operational Control of Drinking-Water Systems; Springer: Cham, Switzerland, 2017. [Google Scholar]
  6. Wang, Y.; Puig, V.; Cembrano, G. Non-linear Economic Model Predictive Control of Water Distribution Networks. J. Process Control 2017, 56, 23–34. [Google Scholar] [CrossRef]
  7. Todini, E.; Pilati, S. A gradient method for the analysis of pipe networks. In Proceedings of the International Conference on Computer Applications for Water Supply and Distribution, Leicester, UK, 8–10 September 1987; pp. 1–20. [Google Scholar]
  8. Rawlings, J.; Mayne, D. Model Predictive Control: Theory and Design; Nob Hill Publishing: Madison, WI, USA, 2009. [Google Scholar]
  9. Mayne, D. Model predictive control: Recent developments and future promise. Automatica 2014, 50, 2967–2986. [Google Scholar] [CrossRef]
  10. Heidarinejad, M.; Liu, J.; Christofides, P. Economic model predictive control of nonlinear process systems using Lyapunov techniques. AIChE J. 2012, 58, 855–870. [Google Scholar] [CrossRef]
  11. Han, H.; Qiao, J. Nonlinear Model-Predictive Control for Industrial Processes: An Application to Wastewater Treatment Process. IEEE Trans. Ind. Electron. 2014, 61, 1970–1982. [Google Scholar] [CrossRef]
  12. Liu, S.; Zhang, J.; Liu, J. Economic MPC with terminal cost and application to an oilsand primary separation vessel. Chem. Eng. Sci. 2015, 136, 27–37. [Google Scholar] [CrossRef]
  13. Zeng, J.; Liu, J. Economic Model Predictive Control of Wastewater Treatment Processes. Ind. Eng. Chem. Res. 2015, 54, 5710–5721. [Google Scholar] [CrossRef]
  14. Pereira, M.; Limon, D.; Muñoz de la Peña, D.; Valverde, L.; Alamo, T. Periodic Economic Control of a Nonisolated Microgrid. IEEE Trans. Ind. Electron. 2015, 62, 5247–5255. [Google Scholar] [CrossRef]
  15. Ellis, M.; Liu, J.; Christofides, P. Economic Model Predictive Control: Theory, Formulations and Chemical Process Applications; Springer: Cham, Switzerland, 2017. [Google Scholar]
  16. Ohar, Z.; Ostfeld, A. Optimal design and operation of booster chlorination stations layout in water distribution systems. Water Res. 2014, 58, 209–220. [Google Scholar] [CrossRef] [PubMed]
  17. Creaco, E.; Alvisi, S.; Franchini, M. A Multi-step Approach for Optimal Design and Management of the C-Town Pipe Network Model. In Proceedings of the 16th Water Distribution System Analysis Conference, Bari, Italy, 14–17 July 2014; pp. 37–44. [Google Scholar]
  18. Schwartz, R.; Housh, M.; Ostfeld, A. Least-Cost Robust Design Optimization of Water Distribution Systems under Multiple Loading. J. Water Resour. Plan. Manag. 2016, 142, 04016031. [Google Scholar] [CrossRef]
  19. Price, E.; Ostfeld, A. Graph Theory Modeling Approach for Optimal Operation of Water Distribution Systems. J. Hydraul. Eng. 2016, 142, 04015061. [Google Scholar] [CrossRef]
  20. Price, E.; Ostfeld, A. Optimal Pump Scheduling in Water Distribution Systems Using Graph Theory under Hydraulic and Chlorine Constraints. J. Water Resour. Plan. Manag. 2016, 142, 04016037. [Google Scholar] [CrossRef]
  21. Di Nardo, A.; Di Natale, M.; Giudicianni, C.; Santonastaso, G.; Savic, D. Simplified Approach to Water Distribution System Management via Identification of a Primary Network. J. Water Resour. Plan. Manag. 2018, 144, 1–9. [Google Scholar] [CrossRef]
  22. Cembrano, G.; Wells, G.; Quevedo, J.; Pérez, R.; Argelaguet, R. Optimal control of a water distribution network in a supervisory control system. Control Eng. Pract. 2000, 8, 1177–1188. [Google Scholar] [CrossRef] [Green Version]
  23. Ocampo-Martinez, C.; Puig, V.; Cembrano, G.; Quevedo, J. Application of MPC strategies to the management of complex networks of the urban water cycle. IEEE Control Syst. Mag. 2013, 33, 15–41. [Google Scholar] [CrossRef]
  24. Wang, Y.; Ocampo-Martinez, C.; Puig, V. Stochastic Model Predictive Control based on Gaussian Processes applied to Drinking Water Networks. IET Control Theory Appl. 2016, 10, 947–955. [Google Scholar] [CrossRef] [Green Version]
  25. Pereira, M.; Muñoz de la Peña, D.; Limon, D.; Alvarado, I.; Alamo, T. Application to a drinking water network of robust periodic MPC. Control Eng. Pract. 2016, 57, 50–60. [Google Scholar] [CrossRef]
  26. Wang, Y.; Ramón Salvador, J.; Muñoz de la Peña, D.; Puig, V.; Cembrano, G. Periodic Nonlinear Economic Model Predictive Control with Changing Horizon for Water Distribution Networks. In Proceedings of the 20th IFAC World Congress, Toulouse, France, 9–14 July 2017; pp. 6588–6593. [Google Scholar]
  27. Wang, Y.; Mũnoz de la Peña, D.; Puig, V.; Cembrano, G. A novel formulation of economic model predictive control for periodic operations. In Proceedings of the European Control Conference, Limassol, Cyprus, 12–15 June 2018. [Google Scholar]
  28. Limon, D.; Pereira, M.; Muñoz de la Peña, D.; Alamo, T.; Grosso, J. Single-layer economic model predictive control for periodic operation. J. Process Control 2014, 24, 1207–1224. [Google Scholar] [CrossRef] [Green Version]
  29. Wang, Y.; Ocampo-Martinez, C.; Puig, V.; Quevedo, J. Gaussian-process-based demand forecasting for predictive control of drinking water networks. In Critical Information Infrastructures Security; Springer: Cham, Switzerland, 2016; pp. 69–80. [Google Scholar]
  30. Rossman, L.A. EPANET2 Users Manual. National Risk Management Research Laboratory: Cincinnati, OH, USA, 2000. [Google Scholar]
  31. Toro, R.; Ocampo-Martinez, C.; Logist, F.; Van Impe, J.; Puig, V. Tuning of Predictive Controllers for Drinking Water Networked Systems. In Proceedings of the 18th IFAC World Congress, Milan, Italy, 28 August–2 September 2011; pp. 14507–14512. [Google Scholar]
  32. Löfberg, J. YALMIP: A Toolbox for Modeling and Optimization in MATLAB. In Proceedings of the IEEE International Symposium on Computer Aided Control Systems Design, New Orleans, LA, USA, 2–4 September 2004. [Google Scholar]
  33. MOSEK ApS. The MOSEK Optimization Toolbox for MATLAB 8.1.0.49. Available online: https://docs.mosek.com/8.1/toolbox/index.html (accessed on 12 April 2018).
  34. Currie, J.; Wilson, D. OPTI: Lowering the Barrier Between Open Source Optimizers and the Industrial MATLAB User. In Proceedings of the Foundations of Computer-Aided Process Operations, Savannah, Georgia, 8–11 January 2012. [Google Scholar]
Figure 1. The relaxation for v i β in unidirectional pipes: original function v i β is plotted in blue bold line, its upper bound is in dashed line and its lower bounds are in dashed dotted lines.
Figure 1. The relaxation for v i β in unidirectional pipes: original function v i β is plotted in blue bold line, its upper bound is in dashed line and its lower bounds are in dashed dotted lines.
Energies 11 00991 g001
Figure 2. Improving nonlinear constraint relaxation for unidirectional pipes.
Figure 2. Improving nonlinear constraint relaxation for unidirectional pipes.
Energies 11 00991 g002
Figure 3. Relaxation for v i v i β 1 in bidirectional pipes: original constraint is plotted in blue bold line, upper bounds are shown in dashed line and lower bounds are shown in dashed dotted lines.
Figure 3. Relaxation for v i v i β 1 in bidirectional pipes: original constraint is plotted in blue bold line, upper bounds are shown in dashed line and lower bounds are shown in dashed dotted lines.
Energies 11 00991 g003
Figure 4. The topology of Richmond water distribution network.
Figure 4. The topology of Richmond water distribution network.
Energies 11 00991 g004
Figure 5. The simulation results of system states with the EMPC-NCR and NEMPC strategies.
Figure 5. The simulation results of system states with the EMPC-NCR and NEMPC strategies.
Energies 11 00991 g005
Figure 6. The simulation results of system states with the EMPC-NCR and NEMPC strategies.
Figure 6. The simulation results of system states with the EMPC-NCR and NEMPC strategies.
Energies 11 00991 g006
Figure 7. The simulation results of control inputs with the EMPC-NCR and NEMPC strategies.
Figure 7. The simulation results of control inputs with the EMPC-NCR and NEMPC strategies.
Energies 11 00991 g007
Figure 8. The simulation results of control inputs with the EMPC-NCR and NEMPC strategies.
Figure 8. The simulation results of control inputs with the EMPC-NCR and NEMPC strategies.
Energies 11 00991 g008
Figure 9. The comparison results of error measurements with the EMPC-NCR and NEMPC strategies.
Figure 9. The comparison results of error measurements with the EMPC-NCR and NEMPC strategies.
Energies 11 00991 g009
Table 1. The fitting values of coefficient C i , j of interconnected pipes.
Table 1. The fitting values of coefficient C i , j of interconnected pipes.
PipeCoeffi. C i , j PipeCoeffi. C i , j PipeCoeffi. C i , j PipeCoeffi. C i , j
7880.009810850.013213040.015618440.0092
7900.011211070.017516380.011818480.0092
7930.010311530.010116450.009918490.0114
7940.009811540.013016530.010718790.0098
8410.009411780.016517400.010019130.0113
9110.011511960.036717520.010319640.0119
9120.011812080.013117530.013019780.0111
9930.010112090.013717830.016820100.0120
10200.009412100.014317930.0189--
10330.010212780.010118320.0145--
10360.010113010.017218420.0097--
Table 2. The hydraulic heads at storage tanks to assess the safety constraints.
Table 2. The hydraulic heads at storage tanks to assess the safety constraints.
TankElevation (m)Volume (m)Hydraulic Head x s (m)
A184.131.02185.15
B2162.03218.03
C258.90.5259.40
D241.181.1242.28
E203.010.01203.03
F235.710.19235.90
Table 3. The KPI results with EMPC-NCR and NEMPC strategies.
Table 3. The KPI results with EMPC-NCR and NEMPC strategies.
MPC Strategy KPI E KPI S KPI M
EMPC-NCR0.69920.26046.7078
NEMPC0.70280.19146.5249

Share and Cite

MDPI and ACS Style

Wang, Y.; Alamo, T.; Puig, V.; Cembrano, G. Economic Model Predictive Control with Nonlinear Constraint Relaxation for the Operational Management of Water Distribution Networks. Energies 2018, 11, 991. https://doi.org/10.3390/en11040991

AMA Style

Wang Y, Alamo T, Puig V, Cembrano G. Economic Model Predictive Control with Nonlinear Constraint Relaxation for the Operational Management of Water Distribution Networks. Energies. 2018; 11(4):991. https://doi.org/10.3390/en11040991

Chicago/Turabian Style

Wang, Ye, Teodoro Alamo, Vicenç Puig, and Gabriela Cembrano. 2018. "Economic Model Predictive Control with Nonlinear Constraint Relaxation for the Operational Management of Water Distribution Networks" Energies 11, no. 4: 991. https://doi.org/10.3390/en11040991

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop