Economic model predictive control based on a periodicity constraint

Economic model predictive control based on a periodicity constraint Ye Wanga, José Ramón Salvadorb, David Mu˜ noz de la Pe˜ nac, Vicenc¸ Puiga,∗, Gabriela Cembranoa,c aAdvanced Control Systems (SAC) Research Group at Institut de Robòtica i Informàtica Industrial (IRI), CSIC-UPC, Automatic Control Department, Universitat Politècnica de Catalunya-BarcelonaTech (UPC), C/ Llorens i Artigas 4-6, 08028 Barcelona, Spain bDepartamento de Ingeniería de Sistemas y Automática, Escuela Superior de Ingenieros, Universidad de Sevilla, Avenida de los Descubrimientos, S/N, 41092 Sevilla, Spain cCetaqua, Water Technology Centre, Ctra. d’Esplugues 75, Cornellà de Llobregat, 08940 Barcelona, Spain a r t i c l e i n f o Article history: Received 25 October 2017 Received in revised form 14 June 2018 Accepted 15 June 2018 Available online 30 June 2018 Keywords: Economic model predictive control Periodicity constraint Recursive feasibility Convergence analysis Water distribution networks a b s t r a c t This paper addresses a novel economic model predictive control (MPC) formulation based on a periodicity constraint to achieve an optimal periodic operation for discrete-time linear systems. The proposed control strategy does not rely on forcing the terminal state by means of a terminal equality constraint and hence it does not require a priori knowledge of a periodic steady trajectory. Instead, at each sampling time step the economic cost function is optimized based on a periodicity constraint over all the periodic trajectories that include the current state. The recursive feasibility and the closed-loop convergence to a periodic steady trajectory are discussed. Moreover, an optimality certificate of this steady trajectory is provided based on the Karush–Kuhn–Tucker (KKT) optimality conditions. Finally, an application to a well-known water distribution network benchmark is presented to demonstrate the proposed economic MPC in which the closed-loop simulation results obtained with a linear model and a virtual–reality simulator are both provided. 1. Introduction Economic model predictive control (MPC) has attracted an increasing attention during the past decade [1–3]. Unlike the conventional MPC formulations [4,5], the main control objective of economic MPC is to optimize an economic performance index without regulating the system to a given steady state or trajectory. Economic cost functions are not necessarily quadratic or positive definite with respect to the given trajectories or references as tracking MPC. Economic MPC has been applied to a variety of industrial applications, such as water distribution networks [6–8], wastewater treatment processes [9], smart grids [10,11] and chemical processes [12,13]. Recently, the closed-loop stability of economic MPC has been widely investigated. Unlike the conventional MPC, the cost function may not be regarded as a quadratic function with respect to a given reference. Hence, the standard stability analysis cannot be directly applied. In [3,14,15], stability analysis of economic MPC has been established under the strong duality and the dissipativity assumptions. Terminal cost and constraint around the optimal steady state ∗Corresponding author. E-mail addresses: [email protected] (Y. Wang), [email protected] (J.R. Salvador), [email protected] (V. Puig), [email protected] (G. Cembrano). are used. In [16], a review is presented for discussing the role of constraints in economic MPC, where the convergence of economic MPC can be enforced by adding terminal constraints. Besides, economic MPC without terminal constraints is studied in [17,18]. Based on the turnpike and controllability properties, the closed-loop convergence is proved. In [19], economic MPC with extended prediction horizon is designed based on an auxiliary controller. An additional term with the auxiliary control law is included in the cost function in order to guarantee the closed-loop convergence. From the application point of view, periodic system behavior appears in some specific systems, such as water distribution networks (WDNs) [20,21] and electrical networks [10]. One specific example stems from the periodic behavior of customer demands in water distribution networks. A WDN generally consists of a large number of hydraulic elements, such as storage tanks, pressurized pipelines, pumping stations (including several parallel pumps) and valves. Economic MPC is suitable for optimizing the economic performance of operations in WDNs, as shown in [6,7,22], but these methods do not take specific advantage of the periodic nature of the consumer demands and energy costs. Taking into account the daily water demand patterns and periodic electricity prices, periodic operations can also be considered in the economic MPC design. A study of periodic economic MPC for the management of a WDN has been reported in [23], where the nonlinear model of a WDN is used with differential and algebraic equations. The closed-loop simulations show that optimal periodic operation achieved by an economic MPC controller is also appearing in the operational management of WDNs. In the literature, the Lyapunov stability of economic MPC with periodic operation is discussed for nonlinear processes in [24]. In [20], a single-layer economic MPC is proposed for periodic operation. An online planner is included in the corresponding optimization problem in order to find the open-loop optimal periodic steady states that are subsequently set as the terminal states for economic MPC. This single-layer strategy for periodic operation is used in [25], where economic MPC without terminal constraints is designed for optimal periodic behavior under the dissipativity and controllability conditions. The closed-loop trajectory converges to the optimal periodic orbit. The main contribution of this paper is to propose a novel economic MPC based on a periodicity constraint for discrete-time linear systems. We formulate an economic MPC optimization problem without setting a terminal state. Hence, it does not need to know a periodic steady trajectory as a priori knowledge. Therefore, the economic cost function is optimized with a periodicity constraint considering all the periodic trajectories including the current state along the prediction horizon. Different from the conventional MPC optimization formulation, the current state is set as a shifted position and not necessarily being the first prediction state. The recursive feasibility and the closed-loop convergence to a periodic steady trajectory are discussed. In order to investigate the optimality of this periodic steady trajectory, an optimal periodic steady trajectory can be obtained by solving an open-loop finite-horizon optimization problem. And an optimality certificate is provided based on the Karush–Kuhn–Tucker (KKT) optimality conditions [26,Section 5.5.3]. Furthermore, an application to the Richmond WDN is presented with the proposed economic MPC controller. Two closed-loop simulation results are provided with linear control-oriented model and with a virtual–reality hydraulic simulator, EPANET. In the case of the linear control-oriented model, no modeling errors and system uncertainties are taken into account. The expected result is that the closed-loop trajectory can converge to the planner trajectory as well as the MPC cost and the closed-loop operational cost. In the case of EPANET, because of the mismatches between the prediction model in economic MPC controller and the simulation model in EPANET, the closed-loop trajectory could not be perfectly periodic but the cost converges to the optimal one provided by the planner. The remainder of this paper is organized as follows. In Section 2, the problem statement of economic MPC is expressed. The economic MPC controller is proposed in Section 3and the recursive feasibility and convergence analysis of the proposed economic MPC are discussed in Section 4. The simulation results of applying the proposed economic MPC into the Richmond WDN are shown in Section 5. Finally, some conclusions are highlighted in Section 6. 2. Problem statement Consider the class of discrete linear time-invariant systems xk+1=Axk+Buk,(1) where x∈Rnxand u∈Rnudenote the system state vector and the control input vector, respectively. Moreover, A∈Rnx×nxand B∈ RRnu×nuare system matrices. For the system (1), system states and control inputs are limited by the following constraints: xk∈X, uk∈U, ∀k∈N,(2) where Xand Uare strictly convex sets of states and inputs. The economic performance of the system (1) is measured by a time-varying economic cost function k(xk, uk, pi), i =mod(k, T)(3) where T∈Z+is a period index and piis a time-varying exogenous signal usually indicating the unit prices, which is stored in a known sequence pas p=pi, i =1, . . ., T, and exhibiting a periodic behavior is implemented using the modulo operator mod(k,T). It is worth mentioning that k(xk, uk, pi)is not necessarily a quadratic function that depends on a sequence of references for tracking. The main control objective is to minimize the closed-loop economic cost measured by k(xk, uk, pi)that is a strictly convex function, ∀k∈Nand the periodicity of this economic stage cost function is given by k(xk,uk,pi) = k+T(xk+T,uk+T, pi) with i= mod(k,T). In this paper, we propose an economic MPC formulation with guaranteeing the closed-loop system convergence to a periodic steady trajectory that minimizes the economic cost while satisfying all the constraints. A procedure to certify that the reached trajectory is optimal with respect to the economic cost provided. In addition, the proposed controller does not lose feasibility even in the presence of sudden changes in the economic cost. 3. Economic model predictive control based on a periodicity constraint In principle, MPC controllers are based on solving a finite horizon optimization problem. If a steady-state trajectory is known, a terminal constraint is included forcing the predictions to reach this steady trajectory at the end of the MPC prediction horizon. While several controllers proposed in the literature are based on a standard terminal region/constraint approach, in this paper a periodic steady-state trajectory is assumed to be unknown in the economic MPC design. We propose a different approach in which the MPC controller seeks to minimize the economic cost function over a single period that includes the current state. Besides, an open-loop finite-time optimization problem is also proposed to find an optimal periodic steady trajectory that will be used for the analysis of the closed-loop convergence. The proposed controller guarantees recursive feasibility and hence the closed-loop convergence even in the presence of sudden changes in the economic cost function, because the constraints of the optimization problem are independent of this cost function. Note that standard approaches that depend on terminal constraints often lead to optimization problems that have to be modified if the economic cost function changes, which in general lead to a more complex control scheme and even to a possible loss of feasibility issues [1,25]. 3.1. Optimal periodic steady trajectory To find the optimal periodic steady trajectory, an open-loop finite-horizon optimization problem, the so-called planner, is presented. Because of the periodic nature discussed above, it can be proved that the infinite horizon problem is equivalent to the following finite horizon optimization problem in which a single period is taken into account [20]. This optimization problem yields the same solution if the time frame to be considered is any period. minimize x0, . . ., xT, u0, . . ., uT−1 JT(x, u)= T−1  i=0 i(xi, ui, pi), (4a) subject to xi+1=Axi+Bui,(4b) xi∈X, (4c) ui∈U, (4d) x0=xT.(4e) Remark 1. Note that in formulation above, the time step i= 0 is chosen as the first step of one period. If a different initial step is chosen,thefunctionswouldbedifferentbutwould lead to an equivalent problem. This choice will affect the proposed economic MPC optimization problem as it will be based on solving a finite horizon optimization problem in a period that starts at some multiple of T, that is, at the same time step used to define the planner. 3.2. Economic MPC formulation The economic MPC strategy is proposed by implementing the following optimization problem. Considering the periodicity, the current state xkis inserted into the shifted position by using the modulo operator mod(k,T). minimize x0, . . ., xT, u0, . . ., uT−1 JT(x, u), (5a) subject to xi+1=Axi+Bui,(5b) xi∈X, (5c) ui∈U, (5d) x0=xT,(5e) xj=xk, j =mod (k, T).(5f) Due to the periodic system behavior, the optimization problem (5) is always initialized from time step i= 0. At each time step, this optimization problem is solved with a fixed prediction horizon of T. Note that the current state xkis not always set as the first state prediction. Let u*be the optimal solution of the optimization problem (5) with the initialization of xk. According to the receding horizon strategy, the optimal control action ukapplied to the closed-loop system at time step∼kis chosen by uk=u∗ j, j =mod (k, T).(6) Remark 2. In the formulations of the optimization problems (4) and (5), the subscript k,∀k∈Ncorresponds to a time instant while the index iwith i= 0, 1, . . .,T−1 refers to a prediction step in the optimization problem. The following two remarks are given including the properties of the proposed economic MPC controller. The detailed discussion and proof will be presented in the next section. Remark 3. Note that the constraints of the optimization problem (5) do not depend on the economic cost function, so recursive feasibility is guaranteed even in the presence of a sudden change. Remark 4. The optimization problem (5) is feasible if there exists a feasible periodic trajectory over a length of Tthat includes the current state xk. This implies that the domain of attraction, that is, the feasibility region of (5) is in general very large, as it is not constrained to reach a specific target in a fixed time as in standard MPC formulations with terminal regions. 4. The closed-loop properties of the proposed controller In this section, the closed-loop properties of the proposed economic MPC are discussed. Recursive feasibility and convergence analysis are standard notions in MPC designs [1]. In the following, these properties are discussed for the proposed economic MPC controller. In particular, under a certain assumption, the closedloop system trajectory converges to the optimal periodic steady trajectory obtained by the planner (4). Theorem 1. The system (1) in closed-loop with the economic MPC implemented by the optimization problem (5) is stable and converges to a periodic steady trajectory. This trajectory is equal to the optimal trajectory obtained from the optimization problem (4), if there exists a time step M > 0 such that for any k ≥M, the dual variables corresponding to the equality constraints (5f) in KKT optimality conditions are zero. Proof. We first prove the recursive feasibility of the closed-loop control system. (Recursive feasibility) If the optimization problem (5) is feasible at time step k∈N, then it is also feasible at time step k+ 1. Let us denote xjand ujwith j= mod(k,T) be feasible solutions of the optimization problem (5) at time step k∈N. All the constraints (5b),(5c),(5d),(5e),(5f) are satisfied at time step k∈N. Thanks to the formulation in (5), as discussed in Remark 2, constraints (5b), (5c),(5d),(5e) are also satisfied at time step k+ 1 since they do not depend on the time step k+ 1. From (5f), we have xj=xk+1, j =mod (k+1, T), which is equivalent to xj+1=xk+1, j =mod (k, T).(7) If the constraint (7) holds, then the optimization problem (5) is feasible at time step k+ 1. Recall the feasible solutions xjand ujwith j= mod(k,T) and the conditions xj=xkand uj=ukhold. From (5b), we can derive xj+1=Axj+Buj, and with the control action ukchosen in (6) and the system (1), the constraint (7) is satisfied. Thus, the optimization problem (5) is also feasible at time step k+ 1. Since the optimization problem (5) is recursively feasible, let us denote its optimal MPC cost as JT(k,x,u) at time step k. By optimality [26],JT(k+ 1, x,u)≤JT(k,x) holds, which implies the cost of the optimization problem (5) is a non-increasing sequence and therefore the closed-loop system is stable. Next, we discuss the convergence of the closed-loop control system. (Convergence analysis) To use the convex analysis for discussions of convergence, the optimization problems (4) and (5) can be reformulated into standard convex formulations. Define the vector z=[x0· · · xTu0· · · uT−1],(8) and the cost function JT(x,u) becomes JT(z). Then, the optimization problem (4) is rewritten to be a standard convex form as follows: minimize zJT(z),(9a) subject to hr(z)≤0, r =1, . . ., m, (9b) gi(z)=0, i =1, . . ., n. (9c) where the functions hrfor r= 1, . . .,mcorrespond to the constraints (4c) and (4d), and the functions gifor i= 1, . . .,nrepresent the prediction model (4b) and the periodicity constraint (4e). Besides, let us denote the optimal solution of the optimization problem (9) (the optimization problem (4)) as zp, that is, the optimal periodic steady trajectory. By the convexity assumption, it follows that there exist dual variables p=⎡ ⎢ ⎢ ⎣ p 1 . . . p m ⎤ ⎥ ⎥ ⎦ , p=⎡ ⎢ ⎢ ⎣ p 1 . . . p n ⎤ ⎥ ⎥ ⎦ , such that the following KKT optimality conditions of (9) hold ∇JTzp+ m  r=1 p r∇hrzp+ n  i=1 p i∇gizp=0,(10a) hrzp≤0, r =1, . . ., m, (10b) gizp=0, i =1, . . ., n, (10c) p r≥0, r =1, . . ., m, (10d) p rhrzp=0, r =1, . . ., m. (10e) Similarly, with the vector zdefined in (8), the optimization problem (5) is also rewritten to be a standard convex form as follows: minimize zJT(z),(11a) subject to hr(z)≤0, r =1, . . ., m, (11b) gi(z)=0, i =1, . . ., n, (11c) Qjz=xk, j =mod (k, T),(11d) where Qjwith j=mod (k, T)is defined based on (5f). Denote the optimal solution of the optimization problem (11) at time step kas zk. Therefore, there exist dual variables k=⎡ ⎢ ⎢ ⎣ 1,k . . . m,k ⎤ ⎥ ⎥ ⎦ , k=⎡ ⎢ ⎢ ⎣ 1,k . . . n,k ⎤ ⎥ ⎥ ⎦ , k=⎡ ⎢ ⎢ ⎣ 1,k . . . nx,k ⎤ ⎥ ⎥ ⎦ , such that the following KKT optimality conditions of (11) hold ∇JT(zk)+ m  r=1 r,k∇hr(zk)+ n  i=1 i,k∇gi(zk)+ nx  l=1 l,kQl mod(k,T)=0,(12a) hr(zk)≤0, r =1, . . ., m, (12b) gi(zk)=0, i =1, . . ., n, (12c) Qmod(k,T)zk=xmod(k,T),(12d) r,k ≥0, r =1, . . ., m, (12e) r,khr(zk)=0, r =1, . . ., m, (12f) where Ql mod(k,T)denotes the lth row of Qmod(k,T)transposed. As discussed above, the cost of the optimization problem (5) is a non-increasing sequence. Taking into account that by assumption that the cost function JT(z)is strictly convex, it is not possible that there exist two consecutive time steps kand k+ 1 such that the costs JT(zk)=JT(zk+1)with zk/=zk+1. Hence, if JT(zk+1)=JT(zk),∀k≥M, the system (1) in closed-loop reaches a periodic steady trajectory, that is zM=zM+1 =· · ·. Without loss of generality, we assume that mod (M, T)=0. Denote zsas this periodic steady trajectory. The solution zsis also feasible for the optimization problem (4). On the one hand, the closed-loop solution zkshould be equal to the optimal solution zs, that is zs=zk,∀k≥M. On the other hand, zsis an optimal solution of the optimization problem (5) such that zs,∀k≥Msatisfies the KKT conditions in (12a),(12b),(12c),(12d),(12e),(12f). If the dual variables in kare zero, then (12d) can be disabled. As a result, zs also satisfies the KKT conditions in (10). Hence, it can be concluded that zs=zp, which means that the closed-loop trajectory zkconverges to the optimal periodic steady trajectory zp,∀k≥M.䊐 Remark 5. For a periodic steady trajectory zs, all the constraints of (11) must be satisfied for all k= 1, . . .,Twith zk=zsand xk=Qjzs with j= mod(k,T). The solution provided by the planner (9) satisfies this condition by definition with dual variables of constraint (11d) equal to zero. Although it is rarely found because the number of constraints is larger than the number of variables, other trajectories may also satisfy the condition. In this case, the closed-loop system may converge to a periodic trajectory different from the one of the planner. To better illustrate Remark 5, the following example is presented, in which the closed-loop system converges to a periodic trajectory different from the planner. Example 6. Consider the following system subject to additive disturbances xk+1=Axk+Buk+Bddk, with system matrices A=0.5 0.5 1 0.25 , B =1 1, Bd=1 0, where dis a periodic known disturbance signal with a period T= 3. The values of these periodic signals are given by dk=diwith i= mod(k,T), where d1=−0.1, d2=−0.2 and d3=−0.1 for i= 1, 2, 3. This system is controlled by the proposed economic MPC. In this example, consider the formulations in (9) and (11) and the quadratic cost function JT(z)=1 2zHz +fzwith H=diag( [111011201110]), f=[0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1], where diag(·) returns a diagonal matrix with diagonal elements defined by its argument. The input must belong to the set U= u∈R:−0.1≤u≤0.1and no constraints are considered for the states. The simulation with this example has been carried out for 60 sampling steps. As shown in Fig. 1, the closed-loop trajectories of both entries x1,kand x2,kof the state states converge to a periodic Fig. 1. The closed-loop state trajectory of the example. Fig. 2. The closed-loop MPC cost of the example. trajectory that is different from the one corresponding to the planner. Fig. 2 shows that the cost of the MPC optimization problem is a non-increasing sequence that reaches a constant value when the system converges to the periodic trajectory. In this case, because this trajectory is different from the optimal one computed by the planner, its corresponding cost is higher than the planner cost. Besides, a measurement (defined by the 2-norm) of the dual variable kis shown in Fig. 3. From this figure, it can be seen that dual variables corresponding to (11d) are not zero at any time since the closed-loop trajectory cannot reach the planner one. 5. Application to the Richmond water distribution network In this section, the proposed economic MPC controller is applied to the Richmond WDN, which is a well-known benchmark available in the link.1 1http://emps.exeter.ac.uk/engineering/research/cws/resources/benchmarks/ operation/richmond.php. Fig. 3. The measure of dual variables kof the example. 5.1. System description The topology of the Richmond WDN is shown in Fig. 4. In this network, there are 6 tanks, 7 actuators (pumps), 11 water demand sectors, 41 non-storage nodes and 41 interconnected pipes. The main objective of operational management of this network is to supply enough water satisfying all the demand sectors with appropriate pressure. From industrial experience, the possible control actions for management of a WDN could be obtained with a periodic behavior because water demands and electricity price usually present daily repetitive patterns. Hence, a periodic economic MPC is suitable for optimal management of a WDN in order to guarantee the multiple control objectives. The common operational objectives for the operational management of WDNs can be summarized as follows [7] •Economic: To provide a reliable water supply with required pressures minimizing operational costs. •Safety: To guarantee the availability of enough water and hydraulic head to satisfy its underlying water demands with enough pressure. •Smoothness: To operate actuators (pumps and valves) of the WDN under smooth control actions. To this end, actuators in the network have to be operated so that these objectives are satisfied taking into account the available information on the future demands and other relevant parameters such as electricity costs. The actuators (pumps and valves) in most WDNs are usually operated in an ON/OFF fashion, defining the control laws as simple switching logic which in general are difficult to tune [27]. Switching actuators lead to hybrid models, where hybrid MPC could be considered as an alternative possibility to address this issue. However, it may be difficult to consider hybrid models in an economic MPC setting. In this case study, to avoid dealing with binary variables, we follow the two-layer control scheme presented in [8], in which a periodic economic MPC controller will decide the mean flow in each actuator of the WDN, and a lower-level layer will implement the switching sequence to guarantee that mean flow. In the case that the pumps and valves can be regulated by controlling directly the frequency or the degree of opening, PID controllers can be used. If only ON/OFF sequences can be implemented, which is a common constraint in real systems, an appropriate open-loop or closed-loop strategy has to be used. In particular, we propose a control strategy in which at the beginning of the sampling period, all pumps and valves are switched on, and are switched off only when the amount of water needed in the period is reached. In this case, it is important to note that the hydraulic head at the nodes depends on the state of the actuators, leading to sudden pressure changes. For this reason, the output variable is defined as the mean hydraulic head at the demand nodes. According to [7], the linear control-oriented model of the Richmond WDN is chosen as follows: xk+1=Axk+Buuk+Bvvk+Bddk,(13a) 0=Euuk+Evvk+Eddk,(13b) where xk∈Rnxdenotes the vector of hydraulic heads at the storage nodes (reservoirs/tanks) as differential states, uk∈Rnudenotes the vector of mean manipulated flows of the actuators (pumps and valves), vk∈Rnvdenotes the vector of mean non-manipulated flows through the interconnected pipes and dk∈Rnddenotes the vector of mean water demands, which are regarded as the measured system disturbances. The water demand dk+ifor i= 1, 2, . . ., Hpalong the prediction horizon Hpis assumed to be predicted as ˆ dk+iby using a suitable short-term forecasting method, such as in [28,29]. The demand forecasts are periodic with a period T. Moreover, A,Bu,Bv,Bd,Eu,Evand Edare system matrices of appropriate dimensions based on the system topology. Eq. (13a) describes the system dynamics, (13b) presents the physical and static relations in the water system by means of mass balance at non-storage nodes. To obtain this control-oriented model (13a) and (13b) of the Richmond network, several simplifications were made in the Richmond WDN. In particular, the flows through the pipes (ID: 1783) and (ID: 1793) were assumed to be equal and tank E was considered to be full at all times, acting as a pipe with a fixed output hydraulic pressure. These simplifications introduce modeling errors between the control model and EPANET, in particular, the errors will be high if tank E is not full. By simulation, it can be demonstrated that because of the head and size difference between tanks D and E, tank E is always full in normal operation. Note also that both tanks are connected by a set of passive links without pumps. Regarding the information available in the link, the sign of the demand at Junction (ID: 777) is modified so that it consumes water. In addition, both the electrical price pattern and the demands are supposed to start at 07:00 when the electrical price pattern starts from the peaks. At each hour, the tank levels are read from the Richmond EPANET model (see Fig. 5), the optimization problem is solved to obtain the desired mean flows, and then a series of EPANET simulations are carried out to decide at which minute, pumps should be turned off depending on the amount of water delivered. Fig. 4. The topology of the Richmond water distribution network. Fig. 5. The Richmond EPANET model. Table 1 The elevations and capacities of storage tanks in the Richmond benchmark. Tank State Elevation [m] Tank capacity [m] A x1,k184.13 3.37 B x2,k216 3.65 C x3,k258.90 2 D x4,k241.18 2.11 F x5,k235.71 2.19 5.2. Economic MPC setup The economic MPC for management of the Richmond network is set up as follows. First, system variables xand uhave the following physical constraints: xk∈X, uk∈U, where Xis defined by the sum of elevations and capacities of the storage tanks, and Ucorresponds to the maximum and minimum pumping flows that can be provided by actuators in the Richmond benchmark. The data of these two sets are given in Tables 1 and 2. Table 2 The maximum and minimum pumping flows through actuators in the Richmond benchmark. Pump Input Minimum flow [m3/s] Maximum flow [m3/s] 1A u1,k0 0.031 2A u2,k0 0.031 3A u3,k0 0.057 4B u4,k0 0.025 5C u5,k0 0.004 6D u6,k0 0.010 7F u7,k0 0.001 Considering a relation between the flow and the energy used by a pump, the economic performance can be estimated by e(uk, pi)=p iuk, i =mod(k, T), where e(uk, pi)directly measures the economical cost of pump operations. p=pi,i= 1, . . .,Tis a sequence of electricity price signals. Fig. 6. The system trajectories of Tank F: blue crossed line is the closed-loop trajectory; green dash line is the planner trajectory; black dash lines are predictions of the economic MPC. (For interpretation of the references to color in this legend, the reader is referred to the web version of the article.) Fig. 7. The optimality certificate: a measure of dual variables corresponding to the equality constraint (11d). Meanwhile, the water levels in the storage tanks are penalized to guarantee safety water volumes by x(xk)= xk 1, where ·1is the 1-norm operator. The cost function x(xk) is used to guarantee the strict convexity of the economic cost function with respect to states, which is not set using the 2-norm due to that the objective is to find periodic water volume evolutions in storage tanks based on the economic cost e(uk, pi)rather than to penalize the storage tanks to be empty. The smoothness control actions can be obtained by optimizing the following cost function s(uk)= uk  2 2, where ukuk−uk−1and ·2is the 2-norm operator. The cost function s(uk)is set to penalize the slew rate of control inputs for the purpose of protecting actuators and maximizing their working lives by avoiding unnecessary ON/OFF switches. In total, the economic stage cost function can be defined as T(xi, ui, pi)=1e(ui, pi)+2x(xi)+3s(ui), where 1,2and 3are the prioritization weights. Therefore, over the prediction horizon Hp=T, the economic cost function ˜ JT(x, u, p) is defined as ˜ JT(x, u, p)= T−1  i=0 T(xi, ui, pi).(14) Note that the parameter 2should be set to be sufficiently small compared to the other weights because the main objective for the management of WDNs is to minimize the economic cost depending on pumping with periodically time-varying electrical prices. Hence, the optimal strategy seeks to accumulate water in the storage tanks when the electricity price is low. For the Richmond WDN, the prioritization weights in the cost function are set as 1= 10, 2= 0.001 and 3= 0.1. The time-varying electricity prices pk, the physical constraint on xkare obtained from the original Richmond system given in EPANET. The physical constraints on ukand the minimum heads around pumps are obtained from running the simulation in EPANET. The short-term water demand forecasts are available on a daily basis. Hence, the period Tof the economic MPC control strategy is chosen as 24 h. Fig. 8. The average economic MPC cost and closed-loop operational cost. (For interpretation of the references to color in this legend, the reader is referred to the web version of the article.) 5.3. Simulation results with the linear control-oriented model To verify optimality certificate in Theorem 1, the first simulation has been carried out with the linear control-oriented model (13) for 5 days (120 h) in MATLAB using the CPLEX solver. The optimization problem (4) has been solved offline to find the optimal periodic steady trajectory. For notation simplicity, this trajectory is called the planner trajectory in the plot. In Fig. 6, the closed-loop trajectories of the simulation are shown (blue crossed line). It can be seen that the system reaches the optimal periodic steady trajectory (green dash line). In addition, at time steps 10, 20 and 30, the predictions of the periodic optimal solution of the economic MPC optimization problem (5) are also shown (black dash lines). It can be seen that the optimal solution of the economic MPC optimization problem (5) is a sequence of periodic trajectories that converge to the planner trajectory corresponding to (4) and that the closed-loop system does not remain in any of the periodic solutions until the optimal periodic trajectory is reached. Besides, for Tank A,B,Cand E, they all reach their planner trajectories earlier than Tank F. Hence, the closed-loop system is stable and converges to the planner trajectory. The values of dual variables k,∀k∈Ncorresponding to the equality constraint (11d) have been extracted online from the optimal solutions. The 2-norm measurements of kat each time step are plotted in Fig. 7. Before the closed-loop system reaches the periodic steady trajectory, all the values in the vector kare far away from zero. From the time step k= 47, all the values in kare approximately equal to zero when the closed-loop system state trajectory reaches the periodic steady trajectory of the planner as shown in Fig. 6. The average MPC cost and closed-loop sliding window cost are shown in Fig. 8. The average MPC cost is computed by ¯ JT= 1 T˜ JT(x∗, u∗, p), where x*,u*and ˜ JT(x∗, u∗, p) derive from the optimal solutions that can be obtained from the linear solver. For k≥T, the closed-loop cost is also computed in a sliding window by ˆ JT=1 T T−1  i=0 Txk−i, u∗ k−i, pj, k ≥T, j =mod(k−i, T), wherexk−iis the state obtained from the linearcontrolmodel.Based on the periodicity, the sliding window is chosen as T. By optimality, the MPC cost (solid line in Fig. 8a) is monotonically decreasing and at time 47 reaches the planner cost (dash line). The closed-loop cost is plotted from the time step k= 25 since the sliding window is chosen to be 24. From Fig. 8b, the closed-loop cost is growing from the beginning since the expense of pumping water to satisfy the demands is measured by the economic cost function. 5.4. Simulation results with EPANET The second simulation has been carried out with the EPANET simulator [30] following a two-layer simulation framework introduced in [8]. The proposed economic MPC controller is applied as the upper layer in MATLAB using the YALMIP toolbox [31]. The optimization problem (5) is solved to find the optimal flow setpoints for all the actuators using the IPOPT solver through OPTI toolbox