Multinode Real-Time Control of Pressure in Water Distribution Networks via Model Predictive Control

Leakage represents a crucial issue in the management of water distribution networks (WDNs). Due to leakage dependence on service pressure, the application of pressure real-time control (RTC) can effectively alleviate the problem. Current RTC implementations rely on a closed-loop control of pressure at a single (local or remote) node of the WDN. In this case, the regulation performance may be good at the selected node, but rather poor across the whole WDN. While conservative choices are usually carried out to mitigate it, this issue becomes extremely relevant in the case of multiple nodes reaching critical values of pressure during the day. This work proposes a novel multinode (MN) RTC approach, which explicitly considers closed-loop control of pressure at multiple WDN nodes. The control scheme is based on a Kalman Filter for state and disturbance estimation, a steady-state auxiliary target calculator, and a model predictive controller for regulation and disturbance rejection. A detailed pressure-driven, unsteady flow model is used to simulate a real WDN under different demand scenarios and assess the performances of the proposed approach, which delivered satisfactory results. Moreover, the MN-RTC approach discussed in this work is suitable for in situ implementation, due to its low computational complexity. Finally, as demonstrated in the simulated environment, the tuning of the control algorithm can be performed by relying on input–output data collected directly from the plant, with no need for a hydraulic simulator.


I. INTRODUCTION
R EAL-time control (RTC) of service pressure plays a fundamental role in the context of water distribution network (WDN) management, allowing for leakage reduction [1], [2], pipe burst abatement [2], [3], [4], and overall infrastructure life extension. These goals can be in fact achieved by suitably reducing pressure excess across the WDN. However, pressure should always be sufficient to fulfill the water demand, which is uncertain and time-varying. To ease control operations, the WDN is first subdivided into homogeneous pressure zones [5], and a pressure control valve (PCV) is installed in each zone [6]. At this point, a hierarchical control scheme is typically adopted, where high-level, model predictive controllers optimize the overall network operation and provide specific setpoints to the lower level RTC controllers [7], [8]. While the former aspect of the control scheme is widely addressed in the automatic control literature (see, e.g., [9], [10], [11], [12]), the design of lower level RTC is recently gaining attention due to the spread of smart WDNs, as promoted by the Water 4.0 approach [13]. In smart WDNs, sensors, actuators, and computing units are connected by wire (e.g., via optical fiber) or by high-end wireless networks (e.g., NB-IoT) to ensure fast and reliable communication. From an RTC perspective, this allows moving from local RTC, where the pressure is controlled right downstream of the control valve [14], [15], [16], to remote RTC, where the pressure can be controlled at any point of the WDN [17], [18], [19], [20]. This represents an important advantage with respect to local RTC, where only the pressure at the valve site is controlled in a closed loop, while the pressure at all the demanding nodes can only be controlled in an open loop [16]. In particular, in remote RTC, the WDN node characterized by the minimum daily pressure (denoted as critical node) is usually chosen as controlled node [6], with the aim of avoiding low-pressure values that may in turn result in unsatisfied users' demand of water. However, a WDN (or WDN pressure zone) may show multiple nodes associated with very low-pressure values (e.g., due to different elevations of the WDN areas). When the demand in the WDN is highly variable in the day, due, for instance, to the presence of commercial or industrial users, two or more nodes may switch the role of critical node, i.e., the node with the lowest service pressure in the WDN: as an example, one node could be the critical node in a certain time slot (e.g., at daytime) and another node could be the critical node in another timeslot (e.g., at nighttime). Under these unfavorable conditions, even remote RTC may fail in properly controlling the pressure and fulfill the users' demand at all times. Again, this is due to the open-loop nature of the control scheme for pressure associated with the nonmonitored nodes. To overcome this limitation, this work proposes a single-inputmultiple-output (SIMO) control design methodology that can effectively cope with the presence of multiple critical nodes; in the rest of this work, this novel approach will be denoted as multinode (MN) RTC, while the standard one will be denoted as single-node (SN) RTC. Specifically, a constrained model predictive control (MPC) scheme is proposed [21], based on a linear parameter-varying (LPV) description of the WDN dynamics [22]. The LPV nature of the system is handled in the MPC by considering a constant prediction for the scheduling variable along the prediction horizon [23], [24]. While, from a theoretical point of view, this approach produces suboptimal control laws (thus the term suboptimal LPV-MPC adopted in the literature), it is typically well suited for practical implementations. As a matter of fact, the resulting finite-horizon optimal control problem (FHOCP) to be solved at runtime is a standard quadratic program (QP), which can be efficiently solved even by relying on low computational power. Moreover, the suboptimality is very often well compensated by the receding horizon (RH) nature of the MPC algorithm [25], [26], [27].
In this work, the proposed control scheme addresses the problems of setpoint tracking and disturbance rejection and is complemented by a parameter-varying Kalman filter (PVKF) for state and disturbance estimation and by a parameter-varying steady-state auxiliary target calculator (SSATC), which computes, instant by instant, an admissible steady-state target for the SIMO system, given the current value of the scheduling variable and the current disturbance estimate [21], [28], [29].
The performance of the control scheme is assessed by performing simulated experiments with a detailed, pressuredriven [30] unsteady flow model of the WDN [31], which is well suited to replace the actual plant for a preliminary in silico evaluation of the approach [32]. In particular, a realistic case study is constructed, based on a real WDN topology and an accurate modeling of the users' demand pattern [33]. Additional tests are carried out to stress the reliability of the control scheme, by simulating realistic scenarios, such as the sudden opening of fire hydrants [34]. The results obtained in silico stress the effectiveness of the novel MN-RTC approach. Moreover, due to its low computational complexity, the proposed algorithm is suitable for in situ implementations, even on low-power computing units.
This article is organized as follows. Section II describes the realistic case study addressed in this work, Section III-A discusses the numerical model of the WDN, and Section III-B discusses the choice of the controlled nodes for MN-RTC. The definition of the WDN working point is given in Section III-C, while a discussion of the model identification procedure is given in Section III-D. The control scheme is thoroughly described in Section III-E. Simulated results are presented and analyzed in Section IV, while further discussion is given in Section V. Finally, Section VI summarizes the main findings of this work.

II. CASE STUDY
The case study adopted for this work is the skeletonized WDN of the town of Castelfranco Emilia, located in Northern Italy. The network consists of 27 nodes (26 demanding nodes and one source node) and 32 pipes. Moreover, two fire hydrants are present at nodes 3 and 13. The complete topology is shown in Fig. 1. Detailed features of network nodes, pipes, and hydrants are reported in [32] and [34]. A PCV with diameter of 250 mm is installed in pipe 26-20, linking the source to the rest of the network. This PCV models a T.I.S. GROUP "NUOVAL" plunger valve (details and datasheet can be found in [35]), which can be equipped with an electric actuator and controlled by a PLC for automatic control purposes. The actuation speed is limited to prevent fast valve operations, which may stress the WDN structure due to the resulting water hammer effect. Specifically, the valve moves from a completely open to a completely closed position in 300 s.
Residential demand profiles are generated according to the bottom-up procedure [33], [36]. The demand profiles follow patterns keeping into account human daily routine. Three main peaks can be identified in the daily profiles: one in the morning, one close to midday, and an evening one, while demand is typically lower and flatter during nighttime. In this work, two different residential demand patterns are considered, leading to two different trends of the total WDN demand [see Fig. 2(a)]: a flatter trend (profile A) and a more peaked trend (profile B). This allows the analysis of RTC robustness with respect to different nodal demand behaviors. In this work, an industrial user is assumed to be present around node 24. This user is modeled by adding to the residential demand a constant offset (0 : 007 m 3 /s plus random, white fluctuations), active between 11 A.M. and 5 P.M. The industrial demand contribution ramps from 0 to its maximum value ad vice versa in 180 s in both demand profiles [see Fig. 2(b)]. The minimum pressure head required for full demand satisfaction for all WDN users is set to 20 m. Finally, the source pressure head profile is reported in Fig. 2(c).
In the context of SN-RTC, node 1 would be chosen as the critical node of the network and, consequently, as the controlled node since it is the node reaching the minimum pressure head during the day [37]. With the aim of guaranteeing a minimum pressure head of 20 m throughout the whole WDN, a reasonable pressure setpoint value for the critical node would be 25 m [20], [37]. For the MN-RTC discussed in this work, both nodes 1 and 24 are chosen as controlled nodes, with pressure setpoints of about 25 m. Specific details about the choice of the two nodes, as well as of the corresponding setpoints, are given in the rest of this article. Here, we only anticipate that the daytime activation of the industrial user close to node 24 induces spatial changes in the pattern of pipe water discharges in the WDN. This causes the switch of the critical node from node 1 to node 24, motivating the need for the MN-RTC approach developed in this work.

III. MATERIALS AND METHODS
This section describes the hydraulic model adopted for simulations, as well as the main steps required for the design of the MN-RTC algorithm introduced in this work.

A. Hydraulic Model for Simulations
This section summarizes the approach adopted for the implementation of the WDN simulated environment; for further details, refer to [20] and [38]. In particular, a pressuredriven [30], unsteady flow modeling approach [31], [32] is adopted and allows for a proper analysis of the hydraulic transients resulting from rapid nodal demand and/or valve setting variations.
For a generic pipe of a WDN, the 1-D unsteady flow equations take the form where h p (m) and Q p (m 3 /s) are the pressure head and the flow discharge along the pipe, respectively, x p (m) is the position along the pipe, t (s) is time, A p (m 2 ) is the pipe crosssectional area, g (m/s 2 ) is the gravity acceleration constant, c (m/s) is the wave celerity, and J p is the friction slope.
To account for leakage from WDN pipes, the following outflow per unit of pipe length q l (m 2 /s) is considered: where α leak (m 2−γ /s) and γ (-) are the leakage coefficient and exponent, respectively. As for leakage evaluation, the exponent γ is set to 1, typical value for plastic pipes [39]. The coefficient α leak (-) is set to 9.4 × 10 −9 m/s to obtain a leakage percentage rate of 20%. Moreover, hydrant outflows q hd (m 3 /s) can be calculated as the outflow from a pressurized orifice based on the emitter equation [5] q hd = C hd h 0.5 p (4) where C hd (m 5/2 /s) is the emitter coefficient, which takes account of the outflow contracted area. If a hydrant is opened linearly in time, C hd increases linearly from 0 to its maximum value C hd,max . On the contrary, if it is closed linearly in time, C hd decreases linearly from C hd,max to 0.
Finally, each pipe friction slope is evaluated as where n [s/m (1/3) ] is the Gauckler-Manning coefficient. In addition, to improve the accuracy of the model and account for the unsteady flow effects, pipe friction slopes are increased using the correction proposed in [40].
The effect of the control valve is modeled by considering no link at the valve site and setting nodal outflow Q up from the upstream end node and nodal inflow Q down into the downstream end node at where A v (m 2 ) is the valve cross-sectional area, g (m/s 2 ) is the gravity acceleration constant, ξ (-)the valve cross-sectional is the valve local head loss coefficient, H v (m) is the head drop in the valve, and α (-) is the valve closure setting, ranging from 0 (fully open) to 1 (fully closed). The valve local head loss coefficient is a growing function of α, as shown in Fig. 3 for the specific valve model adopted in this work. This function is characterized and made available by the valve manufacturer [35]. From this point onward, let Q v denote the flow at the valve site.
In the model implementation, the water hammer partial differential equations are solved by relying on the method of the characteristics [31]. Network pipes are discretized with spatial steps x p , and the hydraulic variables of interest (pressure head and water flow) along the pipes are computed at each time integration step t, with x p and t such that Suitable boundary conditions are assigned in correspondence to source and demanding nodes, where fixed total pressure head and demands are prescribed, respectively. The discretized water hammer equations are coupled with the continuity equation, applied to each node of the WDN. Note that, with this approach, leakage is modeled along pipes, by leaning on the internal nodes of the method of characteristics.
The instantaneous, residential demand at each WDN node is first evaluated using the stochastic bottom-up approach proposed in [33], [36]. The industrial demand, generated as a constant offset with random (white) fluctuations, is then added to the residential demand at node 24 during the time in which this user is active. Finally, the nodal outflow is computed by multiplying the instantaneous demand by the correction factor proposed in [30], to account for the dependence of nodal outflow on service pressure head.
In addition, in order to obtain a more realistic framework, measurement noises are introduced in the WDN model, acting on the measured pressure heads h(t) and on the flow at the valve site Q v (t).

B. Controlled Nodes Selection
The design of any remote RTC algorithm starts with the choice of the controlled node. As briefly introduced earlier in this article, a typical criterion for the choice of the controlled node is the critical node criterion. Specifically, when dealing with SN control schemes, the controlled node is chosen as the node whose pressure head reaches the minimum value throughout the whole network, during the whole day [6], [37]. This worst case approach is usually effective in preventing the pressure head to fall below the minimum desirable value in any point of the network, provided that the pressure setpoint for the control scheme is chosen with some safety margin [38]. However, as discussed in Section I, several nodes, located far away in the WDN topology, may reach critical pressure head values during the day, due to spatial inhomogeneities in both elevation and users' demand. In these scenarios, SN control schemes may fail to maintain sufficient pressure head throughout the WDN, due to the inherent open-loop nature of the control scheme for the nonmonitored WDN nodes.
The MN-RTC framework proposed in this article can overcome this limitation, provided that a careful choice of the controlled nodes is carried out. For this purpose, the critical node criterion needs to be extended for the choice of a set of p controlled nodes. A naive approach would be sorting the WDN nodes according to the minimum pressure heads reached during the day and choosing the p nodes associated with the lowest values of pressure. However, due to the WDN physics, this approach is likely to select the critical node (in the strict sense), plus a set of p − 1 neighboring nodes. This is strongly undesirable since the need for MN-RTC is mainly motivated by a time-varying distribution of the peak demand values through the WDN, or by the presence of multiple, isolated nodes simultaneously reaching very low-pressure head values. Such naive choice of controlled nodes would provide very little advantage to the control scheme, which would continue to operate in an open-loop fashion for controlling the pressure head associated with distant-almost or temporarily critical-nodes.
Additional care is therefore required for the choice of the controlled nodes. In particular, the critical node criterion should be combined with a detailed inspection of the WDN topology. The procedure adopted in this work for the choice of p controlled nodes requires two main steps: 1) identification of p groups of neighboring nodes, whose minimum pressure reaches very low values during the day; 2) application of the critical node criterion to each group of neighboring nodes from step 1. As an example, Fig. 4 shows the minimum daily nodal pressure head in the no-control scenario for some potentially critical nodes in the case study. Based on the spatial distribution of the nodes in the WDN, two critical nodes are finally chosen, i.e., nodes 1 and 24.

C. Working Point Definition
This section discusses the definition of the working point of the system in the case of MN-RTC. Consider the multi-inputmultioutput (MIMO) system defined by the following input signals.

3) D(t) is the vector of water demands (m 3 /s).
The output signals are given as follows.
1) h(t) is the vector of the measured pressure heads at the controlled nodes of the WDN (m). 2) Q v (t) is the flow measured at the valve site (m 3 /s). Let ξ(α(t)) be the control variable, h(t) be the controlled variables, and H (t) and D(t) be stochastic disturbances acting on the process. Moreover, let D tot (t) be the overall demand of the WDN, i.e., where D i is the water demand at node i and N nodes is the number of demanding nodes in the WDN. The average values of typical H (t) and D i (t) profiles are considered as input to the system for the definition of the WP. Note that this information is usually available to the WDN manager. Fig. 2(a) and (c) shows the profiles of D tot (t) and H (t) adopted for simulations, respectively. Let the tuple WP = (ξ , H , D, h, Q) represent the working point for the MIMO system. In this work, the values of ξ and h are defined by simulated experiments performed on the WDN. In particular, the procedure consists in adjusting the valve closure α (i.e., the local head loss ξ ) until the minimum of h reaches the desired pressure head. The value of ξ associated with such valve closure is adopted as ξ . The corresponding pressure heads at the controlled nodes then define h. Note that, in a real scenario, this procedure can be directly performed on the real plant, by manual trial and error adjustments of α, or by means of SISO control schemes, as discussed in [41]. 1 Moreover, preliminary computations, based on a (possibly steady state) hydraulic model, can provide a convenient starting point for the search of ξ .
Remark: The choice of ξ as control variable, opposed to α, is motivated by the strong nonlinearity of the ξ(α) curve. As underlined in [42], when the system moves far away from its nominal WP, with α ≈ 1 in particular, the increase in the process gain may result in instability of the closedloop system. Since the ξ(α(t)) curve is typically available and invertible, the problem can be faced by exploiting ξ as a control variable and introducing a nonlinearity inversion block in the loop to compute at runtime the corresponding value of α [38].

D. Process Model Identification
This section summarizes the identification of a black-box LPV process model [22], to describe the WDN dynamic behavior. It should be remarked that the choice of explicitly modeling the transient response of the WDN is motivated by the detailed analysis performed in [38], [43], [44]. The analysis stresses the fact that the presence of unmodelled dynamics may result in wide pressure oscillations or even lead to closed-loop instability [20], [43], as the WDN drifts away from its nominal operating point, due to fluctuations in the users demand. This is extremely important, as the occurrence of the aforementioned events could not just impact on the service to end users but also stress or damage the WDN infrastructure. Note that this issue mainly occurs with low-level controllers in charge of pressure regulation, as the one proposed in this work. Moreover, as stressed in [42], standard robustness analysis may not provide accurate results if the WDN dynamic response is not modeled. On the other hand, even simple or partial information about the WDN transient behavior [45] could sensibly improve the design of more robust and efficient control algorithms [46]. The approach was originally introduced and discussed in detail in [41] and [45], which can be consulted for further details; further details on the linear identification techniques can be found in [47].
The procedure adopted in this work consists of three main steps.
1) Identification of Linear, Local Models of the WDN Dynamics Around WP: Consider a parametric, linear, local model, whose input is the variation signal δξ(t) = ξ(t) − ξ and output is the variation signal δh(t) = h(t)−h. Input-output data can be obtained by exploiting a nonlinear WDN simulator (hydraulic modeling discussed in Section III-A), if available, or by means of experiments performed directly on the real plant. While the experimental design can be focused on maximizing the information contained in the dataset in the former case, particular care is required in the latter case since the experimental session should not interrupt the service to users and should not stress the WDN structure [41]. Input-output identification techniques can then be applied (e.g., by means of MATLAB Identification Toolbox [48]) to compute the parameter values that provide the best fit between model prediction and identification and test data. 2) Discretization and Realization: If the linear model from step 1 is identified as a continuous-time one (e.g., as in [38] and [41]), it then needs to be discretized with a proper sampling time. It should be remarked at this point that the high-frequency behavior of the system plays a very important role and should be properly described by the linear process model to reduce the risk of instability [42], [46]. The choice of the sampling time should be therefore carried out to preserve the description of high-frequency oscillations in the model response. In addition, the control sampling time should also be consistent with the desired closed-loop response of the system [46], [49]. As a matter of fact, closedloop settling time plays an important role since sudden pressure drops occur regularly (e.g., due to the presence of the industrial user) but may also occur due to the opening of fire hydrants or sudden pipe breaks. In this scenario, pressure should be promptly restored to avoid problems to end users. In the case of sudden pressure increase due to a demand drop, pressure should be instead quickly lowered to reduce the stress on the WDN infrastructure [20], [34]. Also, note that, as discussed later in this article, the proposed control scheme requires the solution of two QPs at each sampling instant and is therefore computationally tractable for sampling times in the order of seconds. Finally, the discrete-time model is converted to a state-space form, and the transport delay terms that characterize the system response [45] are replaced by a proper number of poles at the origin of the complex plane. At this point, the local, linear model is completely defined as where δx ∈ R n is the state of the black-box local model, δξ ∈ R is its input, and δh ∈ R p is its output; A ∈ R n×n , B ∈ R n×1 , and C ∈ R p×n . 3) Extension to the LPV Structure: The valve equation (6) shows a quadratic dependence of the pressure loss induced by the valve on the flow through the valve itself Q v . In particular, previous works dealing with SN-RTC underline that a static description of this dependence can be very effective in improving the associated control design [17], [38], [45]. The static gain µ of (9) is given by µ = C(−A) −1 B. The input matrix B can be made parameter-dependent to accommodate the static nonlinearity as follows: where the measured flow at the valve site Q v (k) is measured online and assumes the role of scheduling parameter [22]. To conclude, the LPV system s modeling the process dynamics is given by (10).

E. Control Design Methodology
The control scheme proposed in this work is shown in Fig. 5 and consists of three main elements.
1) A PVKF for state and disturbance estimation.
2) An SSATC for the computation of admissible steady states for the process model, given the current disturbance estimate obtained via PVKF, and the current measurement of the flow at the valve site. 3) A model predictive controller for the computation of optimal control moves, given a suitable cost function for regulation to the current steady state and a set of physical and performance constraints to be fulfilled at each time instant. Moreover, two low-pass filters are introduced in the scheme, with the aim of smoothing the dynamics of the two scheduling variables (Q v and α). Each element of the scheme is in fact based on the LPV system s , which is updated at every time instant according to the current value of the scheduling variable Q v . In addition, the MPC controller features a further gain scheduling policy based on α, to improve the tradeoff between error and cost of control over a wider range of operating conditions [38]. The choice of removing high-frequency components from the scheduling signals is common in the context of LPV-based gain scheduling (see, e.g., [22], [38], [50], [51]) and typically results in improved control performances.
Each element of the control scheme is discussed in detail in the following.
1) State and Disturbance Estimation: This section describes the PVKF adopted for online estimation of the states of s . Moreover, the control scheme proposed in this work accounts for the presence of exogenous disturbances affecting the WDN (source pressure head and users' demand) as black-box, nonmeasurable process disturbances d ∈ R p acting on the output of s , i.e., The PVKF is therefore exploited for both state and disturbance estimations. To this end, the process model s is augmented to include the effect of disturbances, which are modeled as integrating disturbances, driven by the white Gaussian noise. Let KF be the corresponding augmented system included in the PVKF and where the system matrices are augmented to include the disturbance dynamics as follows: with 0 r ×c ∈ R r ×c the zero matrix and I r ×r ∈ R r ×r the identity matrix. Letδx KF (k|k − 1) be the one-step prediction of δx KF (k), given the knowledge of system inputs and outputs up to the time instant k − 1. Then, the PVKF estimator is described by the following equations: where the covariance of the estimation errorP(k|k − 1) is the solution of the Riccati equation reported in the following equation: initialized according toP Note that, due to the system augmentation introduced beforê whered(k|k − 1) represents the one-step prediction of d(k) given the knowledge of system inputs and outputs up to the time instant k − 1.
Remark: It should be stressed that, with the choice of the LPV structure described in Section III-D, both the estimator gain L and the covariance of the estimation errorP do not show any dependence on the scheduling parameter, which only appears in the system matrix B KF . The following result on the convergence of the proposed PVKF therefore holds.
Theorem 1 [21]: Let B q be a partition ofQ such that Q = B q B ′ q . If the pair (A KF , B q ) is reachable and that the pair (A KF , C KF ) is observable, then the proposed state estimator (21) is optimal with where P is the unique positive definite solution of the stationary Riccati equation Moreover, the estimator is asymptotically stable, that is, all the eigenvalues of (A KF − LC KF ) have modulus less than 1.
2) Steady-State Auxiliary Target Calculation: At each time instant k, given the pressure setpoint values δh sp = h sp − h, the current output disturbance estimated(k|k), and the current scheduling parameter Q v (k), the following constrained QP allows the computation of an admissible steady state SS(k) = (δξ ss (k),d(k|k), δx ss (k), δh ss (k)) for the system s,d : where ∥m∥ 2 M denotes the square norm of a vector m weighted by a matrix M, where R > 0, and Q ss ≥ 0 is a diagonal, positive definite matrix. Moreover, δξ min and δξ max are lower and upper bounds for the control variable δξ , respectively, and δh ss min and δh ss min are lower and upper bounds for the steady-state values of the controlled variable δh, respectively. Note that output bounds (32) can be implemented as soft constraints to preserve the feasibility of optimization problem, while input bounds (31) model physical limitations on the valve opening and closure and are therefore implemented as hard constraints. The matrix Q ss can be used to assign different weights to the distance of the steady-state outputs δh ss from their desired setpoints δh sp , for different controlled nodes. Finally, in the case of multiple equilibria yielding the same steady state/setpoint distance, the presence of R ss favors the choice of the minimum norm steady-state input and guarantees the uniqueness of the choice [21].
3) Model Predictive Control: Once a feasible steady-state target is available, a constrained MPC scheme is used to enforce optimal regulation. The MPC is based on an augmented version of the LPV process model s,d , where an integral action is introduced at the system input. This allows weighting the derivative of δξ in the cost function, thus providing a better measure of the energy required to perform regulation and of the wear of the PCV [6]. Let ξ denote the discrete-time derivative, and let MPC be the augmented system included in the MPC and where the system matrices are with T s the controller sampling time.
The following cost function is considered: where ξ [k,...,k+N −1] = [ ξ(k), . . . , ξ(k + N − 1)] is the future control sequence defined from instant k, δx MPC ss = [δx ′ ss δξ ss ] ′ is the steady-state reference for the state of the augmented system MPC , Q y = Q ′ y ≥ 0 is the weight for the regulation error, R(α(k)) > 0 is the (parametervarying) weight for the control variable, S(Q v (k), α(k)) = S ′ (Q v (k), α(k)) > 0 is the (parameter-varying) weight for the terminal state, and N > 0 is the prediction horizon. In this work, the weight for the control variable varies according to the following scheduling law: where R max > 0, R min > 0, c R > 0, and l R > 0 are design parameters. As an example, the scheduling law adopted in this work is shown in Fig. 6. As far as the terminal weight S(Q v (k), α(k)) is concerned, a gain-scheduled version of the infinite-horizon linear quadratic regulator (LQR) is considered, with associated parameter-varying, optimal cost P LQR (Q v (k), α(k)). Specifically, at each time instant k, P LQR (Q v (k), α(k)) can be computed by solving the stationary Riccati equation reported in (42), where the LQR weights, Q LQR and R(α(k)) LQR , are chosen to match instant by instant those of the MPC [21], i.e., To sum up, the following FHOCP is solved at each time instant k: subject to the dynamics of MPC (33) and subject to the following constraints: Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. with and where it is assumed The closed-loop control action, ξ ⋆ (k), is then obtained as the first element of the optimal control sequence, according to the RH principle.
Finally, the actual local head loss coefficient to be achieved at time k, ξ ⋆ (k), is computed as and the corresponding valve closure value, α ⋆ (k), is computed by inverting the valve curve ξ(α) and evaluating it at ξ ⋆ (k).
If necessary, the integration should be suitably saturated to account for the presence of bounds acting on δξ(k).
It should be remarked that the proposed MPC scheme belongs to the class of suboptimal LPV MPC schemes, due to the assumption of a frozen guess for the scheduling parameter along the horizon N . As already introduced in Section I, this allows formulating an FHOCP belonging to the class of QPs, thus easily solvable in real time, even with long horizons. This comes at the price of (possibly nonnegligible) mismatch between plant dynamics and model prediction, which is, however, well compensated by the RH nature of the control scheme in terms of control performances.
In addition, note that the state constraints acting on the local head loss coefficient (45) should be consistent with their steady-state counterparts (31), while the output constraint (46) may assume different values from their steady-state counterparts (32). Moreover, in order to ensure its feasibility of the FHOCP, output constraints (46) are implemented as soft constraints.
To mirror the presence of the valve speed limit, (44) is introduced to impose upper and lower bounds on the time derivative of the local head loss coefficient. It should be remarked that the straightforward implementation of the valve speed limit as bound on the local head loss coefficient would result in a nonlinear constraint, due to the nonlinear relation linking the two variables. To overcome this limitation and maintain a QP formulation of the FHOCP, the constraints are implemented as linear, parameter-varying bounds. Specifically, at each time instant k, the upper and lower bounds ( ξ max (k) and ξ min (k), respectively) are obtained by computing the approximate derivative (right and left, respectively) of the ξ(α) curve via finite difference. The finite difference discretization steps are chosen equal to the maximum (positive and negative) variation allowed for α, according to its speed limits (opening and closure).
Finally, the nonlinearity of the valve curve also motivates the choice of introducing the scheduling law for R, defined in (39). In particular, when the operating point of the WDN moves toward high flow values, α assumes values close to zero. In this region, the curve ξ(α) becomes almost flat, meaning that even a small control signal ξ requires a wide variation of α to be realized, thus resulting in strong oscillations in the valve closure. As demonstrated in [38], a proper scheduling policy can be effective in balancing the tradeoff between controller aggressiveness and cost of control across a wide range of operating conditions. 4) Low-Pass Filtering of Scheduling Parameters: Recall that the process static gain dependence on the flow at the valve site Q v (t) arises steady-state considerations about the process under control [17], [38]. However, in real-life situations, the behavior of the WDN is mainly driven by the demand at the different nodes, which show daily patterns as those shown in Fig. 2(a). For this reason, it is necessary to introduce a low-pass filter to remove the oscillations in Q v (t) arising from both high-frequency demand oscillations and from the water hammer effect in the WDN pipes. In addition, a part of such oscillations are transferred to α(t), and thus, the same considerations hold for the scheduling of the MPC weight R(α(t)). In both cases, the gain scheduling is therefore implemented by relying on Q v LP (t) and α LP (t), the low-pass filtered versions of, respectively, Q v (t) and α(t).

F. Offset-Free Regulation of Average Pressure
As discussed earlier in this article, the typical setup for MN-RTC consists of a single PCV and a set of p controlled nodes, with p > 1. The number of control variables is, therefore, lower than the number of controlled variables; thus, the system does not fulfill a necessary condition to guarantee an offset-free regulation [21], [52] of pressure at every controlled node.
However, it is sometimes desirable to favor an offsetfree regulation of the average of pressure, to ease higher level management of pressure distribution across the WDN. This goal can be achieved by defining the average pressure deviation from the WP one as follows: where 1 1× p ∈ R p is the one vector.
The following constraint is then included in the SSATC In this way, at each time instant k, the SSATC tries to compute an admissible steady state associated with null average pressure deviation. As discussed in Section III-E2, to preserve the feasibility of the associated QP, also (56) should be implemented as soft constraint. However, priority should be given to the relaxation of (56) with respect to the relaxation of output constraints (32).

IV. RESULTS
This section describes the application of the MN-RTC methodology proposed in this work to the case study of the Castelfranco Emilia WDN, as described in Section II.

A. Selection of Controlled Nodes
The first step of the methodology consists of the choice of controlled nodes, as discussed in Section III-B. To this end, a whole-day simulation is performed by considering realistic users' demand and source pressure head profiles, while the PCV is set to a constant value. The minimum daily pressure is recorded for each WDN demanding node and is shown in Fig. 4(a).
With the help of a topological representation of the WDN [ Fig. 4(b)], it is possible to identify two main areas of the network characterized by low-pressure values. Specifically, the first one is covered by nodes 1-3, while the second one is covered by nodes 24 and 25.
The application of the critical node criterion to the two areas leads to the choice of nodes 1 and 24 as controlled nodes. Note that the aforementioned areas could be extended to include nodes 14-17, which are anyway associated with slightly higher pressures and would not alter the choice of the controlled nodes.

B. Working Point Definition and Process Model Identification
Once the choice of controlled nodes is carried out, it is necessary to define a working point for the system under control. Note that the system outputs correspond to the measured pressures at the controlled nodes, arranged in increasing order of node labels. Following the directions given in Section III-C, the working point WP results: Different step response experiments are then performed around WP, to build up a dataset for the identification of a linear, local process model. At this stage, data are sampled at a sampling rate of 1 Hz, assumed as the maximum sampling rate available for the system. A continuous-time transfer matrix G(s), linking δξ to δh, is identified with the prediction error method from MATLAB Identification Toolbox [48]. The structure of the two transfer functions composing G(s) is characterized by a relative degree equal to 1 and by the presence of a transport delay term [41], [45].
In order to define the best model order for each of the two transfer functions, the original dataset is first split into identification and validation subsets, and a set of candidate models is identified using the former dataset and tested on the Fig. 7. Comparison of data (blue dashed line) and model predicted (red solid line) measured pressure δh(t) in case of a negative step variation of δξ(t). Top: pressure variation δ h 1 (t) (node 1). Bottom: pressure variation δ h 2 (t) (node 24). latter. The model resulting in the best performances on the validation dataset is then chosen as the definitive one. In the case of models delivering similar performances, the most parsimonious one is chosen. In particular, all the goodnessof-fit metrics adopted in this work (root-mean-square error (RMSE), coefficient of determination (COD), and FIT) [47] support the choice of order equal to 7 for both transfer functions.
The overall transfer matrix G(s) is reported as in (57), shown at the bottom of the page. Fig. 7 shows a comparison of data and model prediction in the case of negative step variation of δξ . The frequency response associated with each transfer function is shown in Fig. 8  This has to be expected since the static behavior of the plant is mainly due to the pressure loss at the valve site, which is in turn generated by the same PCV. A zero-order-hold (ZOH) discretization is applied to G(s) to obtain a discrete-time representation. The sampling time is relaxed to T s = 3 s. This choice still provides a sufficiently detailed characterization of the high-frequency behavior of the system while allowing for further time for the computation of the control action. The transport delay terms are then transformed into poles at z = 0, and a state-space representation for the discrete-time process model is computed. Finally, the static dependence on the flow at the valve site Q v is included by parameterizing the system matrix B s , as discussed in Section III-A. This defines the LPV process model s adopted for the design of the overall control algorithm. In particular, s has n = 21 states and p = 2 outputs.

C. Control Design
This section discusses the design of the main components of the MN-RTC scheme for the specific case study analyzed in this work. In particular, the tuning of the design parameters is carried out as follows. First, a rough estimate of disturbances acting on the measured pressure is obtained by performing a simulation of the WDN, with the valve closure set at its working point value. The disturbance estimate is computed at each time instant as the difference of the recorded pressure heads and the corresponding working point values. At this point, the closed loop based on the identified LPV model and the disturbance estimates can be simulated to evaluate the desired performance metrics. This allows for a trialand-error tuning of the main parameters of the algorithm. Note that this approach does not require any closed-loop simulation involving the complex hydraulic model discussed in Section III-A. The results of the tuning phase are described in the rest of this section.
1) State and Disturbance Estimation: As discussed in Section III-E1, a PVKF is designed by considering integrating disturbances acting at the system outputs and properly augmenting s to provide their estimate, together with the estimate of nonmeasurable states. The design parameters for the PVKF are set tõ Note that, due to the black-box nature of s , no specific tuning was carried out on the elements ofQ corresponding to the states of the system. On the contrary, the tuning process focused on the elements corresponding to the disturbances.
2) Steady-State Auxiliary Target Calculation: Next, the SSATC is setup according to the following design parameters: In particular, δξ min and δξ max correspond to a saturation of valve closure, which is restricted to the interval α ∈ [0; 0.9], while δh ss min and δh ss max are chosen as 0.8δh min and δh max , respectively. Finally, with the aim of favoring an offset-free regulation of the average pressure, the SSATC is setup according to the procedure discussed in Section III-F.
3) Model Predictive Control: Following the approach described in Section III-E3, the system s,d is extended to include an integral action on the control variable, and the FHOCP is formulated according to the following design parameters: (67) Fig. 6 shows the value of R as a function of α for this specific tuning. Specifically, the tuning parameters are chosen to smoothly raise the value of R as the valve closure approaches α = 0, starting from α ≈ 0.5. In this interval, the curve ξ(α) gets flatter (see Fig. 3), and a wide movement of the valve shutter is required to produce a small variation in the valve head loss coefficient. The proposed tuning is therefore chosen to reduce the controller aggressiveness and reduce the associated valve motion. For this case study, no specific tuning is carried out for each of the two controlled variables. The prediction horizon N is selected so as to cover the whole open-loop settling time of the plant, as shown in Fig. 7. As for constraints, δξ min and δξ max correspond to the valve saturation, while δh min and δh max correspond to the minimum and maximum pressure heads of h min = 20 m and h max = 40 m, respectively, for both controlled nodes. As already discussed in Section III-E3, ξ max (k) and ξ min (k) are computed online, at each time instant k, to fulfill the constraint on the valve speed Hz.
4) Low-Pass Filtering of Scheduling Parameters: Let LP(s) be the transfer function of the low-pass filters used to remove oscillations in Q(t) and α(t) Based on the analysis performed in [38], a filter time constant T LP = 180 s is used in this work. Both filters are implemented in discrete time, with a sampling time of 3 s.

D. Closed-Loop Tests
With the aim of evaluating the performances of the proposed MN-RTC scheme, whole-day, closed-loop simulations are performed with demand profiles A and B. The results of closed-loop simulations are evaluated according to the following metrics.
The regulation error evaluates the proximity of pressure head to the desired setpoint at each controlled node of the WDN.

2)
| α(k)|(−), where α(k) = α(k) − α(k − 1). The control cost impacts on the energy required to perform regulation and on the wear of actuators. Note that, for the sake of performance evaluation, all signals are sampled at the minimum sampling time assumed for the system (1 s). Moreover, to further evaluate the applicability of the proposed approach to a real WDN, both average and maximum computing times required for all closed-loop operations (i.e., filtering, state and disturbance estimation, auxiliary target, and closed-loop control action computation) are recorded for both simulations. In particular, simulations are carried out on an Intel 2 Core 3 i7-10700 CPU 2.90 GHz processor. Optimization problems are implemented with the help of YALMIP [53] and solved by means of MOSEK [54]. Fig. 9(a) shows the results of a whole-day, closed-loop simulation with demand profile A. For the sake of comparison, Fig. 10(a) shows the results of a simulation performed with the SN-RTC algorithm from [38]. In this case, node 1 is chosen as the only controlled node in the WDN. Both algorithms share similar performances, while the industrial demand is not active. However, in the presence of industrial demand, the SN-RTC algorithm does not properly compensate for the strong pressure loss at node 24. This can be expected, as node 24 is not monitored with the SN-RTC scheme, and the effect of the industrial demand on the controlled node is negligible for the control loop to adequately raise the corresponding pressure head. On the contrary, the MN-RTC immediately reacts to the activation of the industrial demand by allowing an increase in the pressure head at node 1 and preventing that at node 24 to fall down to undesirable values. Fig. 9(b) shows the results of a whole-day closed-loop simulation with demand profile B, which is designed to test the behavior of the control loop over a wider range of operating points. As a result, the valve spans from almost completely open to almost completely closed during the day. As in the previous scenario, the MN-RTC scheme quickly reacts to the activation of the industrial demand at node 24 and finds a balanced compromise for pressure at nodes 1 and 24. Again, the SN-RTC algorithm [see Fig. 10(b)] ignores the sudden pressure loss occurring at node 24, and keeps regulating the controlled pressure at the usual value. A similar situation occurs during the main demand peak around 9 a.m., with the MN-RTC algorithm allowing for a slightly higher reduction of pressure at node 1, in order to keep pressure at node 24 at a lower value, with respect to the SN-RTC algorithm.
For both demand profiles, the MN-RTC algorithm manages to deliver satisfactory performances across all operating points and quickly reacts to the activation of the industrial demand at node 24. It is also interesting to note that, in both cases, the choice of the steady-state auxiliary targets satisfactorily fulfills the offset-free requirement of the average pressure deviation (see Fig. 11).
For a quantitative assessment, Table I summarizes the performance metrics for whole-day closed-loop simulations.  In general, note that both regulation error at each node and control cost values are aligned with the results obtained with SN-RTC schemes, with the MN-RTC approach trading some regulatory performance at node 1 to improve it at node 24 while maintaining similar control costs. This is quite interesting since the proposed scheme controls two nodes with a single control variable. Moreover, for MN-RTC, it is also interesting to note that the regulation error is extremely similar at both nodes, as favored by the choices of Q ss and Q y as identity matrices. Finally, despite the long prediction horizon and the high number of states, the computing times are well below the sampling time of 3 s chosen for the loop elements. Note that no specific implementation was carried out to speed up the computations.
To further assess the ability of the proposed MN-RTC scheme to respond to sudden variations of the demand across the WDN, several tests are performed by simulating a sudden opening of a fire hydrant, combined with demand profile A. For all simulations, a fire hydrant is kept open for 1600 s. In particular, Fig. 9(c) shows the results of a closed-loop simulation where the fire hydrant at node 3 is activated during daytime. Note that this coincides with the time of the day when the industrial demand is also present, and thus, the overall flow reaches very high values. This in turn results in the PCV to operate close to saturation. Due to this physical constraint, the output constraint regarding pressure at node 1 is slightly violated. However, the pressure at node 24 is properly regulated to its setpoint. As far as transients are concerned, a quick settling time (about 180 s) is obtained both at hydrant opening and closure. Fig. 9(d) shows the results of a closed-loop simulation where the fire hydrant at node 3 is activated during nighttime. In this second case, both pressure heads are properly regulated at their setpoints, again with a settling time in the order of few minutes. Similar results are obtained with the activation of the fire hydrant at node 13. Results are reported in Fig. 9(e) and (f), which depict a daytime and a nighttime simulation, respectively.

V. DISCUSSION
The results presented in Section IV highlight the effectiveness of the proposed control scheme and suggest that it can be effectively applied to a real WDN. As a matter of fact, it should be stressed that the methodology proposed in this work is ready to be carried out directly on the plant, with no need for a detailed hydraulic model of the WDN nor for any estimate of users' demand signals. As already introduced in Sections III-C and III-D, both working point and process model identification phases can be in fact carried out in situ, due to the approach discussed in [41]. Moreover, the same tuning procedure adopted in this work and discussed in Section IV-C can be applied when directly working in situ, as it is only based on the LPV model and a data-driven estimate of disturbances. In addition, it is worth stressing that the proposed control methodology can be applied to an arbitrary number of controlled nodes and that the computational complexity of this methodology scales as that of LPV model, rather than that of the hydraulic one. This makes it suitable for large WDNs. Finally, the low computational complexity of the optimization problems involved in the computation of the control action allows for an online implementation of the control algorithm.
While the possibility of directly applying it to a real WDN is indeed a strong advantage, some of the steps of the proposed methodology leave room for further research. For example, the aforementioned tuning procedure could be optimized and automated, by formulating it as an optimization problem. Moreover, the estimate of disturbances, on top of which the tuning is based, could be further refined, e.g., by means of a further experimental session and the use of an observer (as the PVKF introduced in this work). A periodic [55] model of the disturbance could also be combined with the MPC algorithm for tracking periodic trajectories (e.g., [56], [57]), which are successfully applied in the context of high-level management of WDNs [58]. In addition, while this work focuses on underactuated pressure control system, future works can investigate possible issues arising from square or overactuated systems, for which LASSO-MPC could be a useful tool [59]. Another interesting research direction would be the design of a data-driven algorithm for the choice of controlled nodes. As a matter of fact, the guidelines proposed in this work require a domain-expert decision-making process, carried out by visual analysis of the WDN topology. An automated approach for controlled node selection would therefore ease the application of the proposed methodology to large WDNs. Finally, as the control algorithm investigated in this work explicitly targets some very specific needs of the application (weighting of system outputs, gain scheduling of input weight, and time-varying input rate constraints), a complete robust stability proof would require widely conservative arguments. Therefore, the possibility to guarantee robust stability by means of specific output feedback, LPV-MPC algorithms, could be investigated in future works, to discuss whether the associated performances could compare with those achieved with the present control scheme.

VI. CONCLUSION
This work focuses on MN-RTC of service pressure in WDNs. Standard RTC is in fact based on closed-loop control of pressure at a single node of the WDN. Specifically, the node characterized by the minimum daily pressure (critical node) is typically chosen as a controlled node. Service pressure at all other nodes is instead controlled in an open loop. This unavoidably leads to conservative setpoint values for the controlled pressure. Moreover, as discussed in this article, multiple nodes may simultaneously play the role of critical nodes or several critical nodes can alternate during the day, e.g., due to spatially unbalanced users' demand through the WDN. The proposed MN-RTC scheme allows for simultaneous control of pressure at several WDN nodes and overcomes these issues. The scheme is based on three main elements: a Kalman Filer for state and disturbance estimation, an SSATC to compute admissible references for the nonsquare system under control, and a constrained model predictive controller to regulate the system at the admissible setpoint. The proposed approach is tested on a detailed hydraulic model of an existing WDN, and the results obtained in simulations highlight the effectiveness of the approach. In particular, the novel MN-RTC scheme provides good disturbance rejection at all the considered nodes, in the presence of spatially unbalanced users' demand. Moreover, the proposed scheme is well-suited for online implementations due to its low computational complexity. Finally, the methodology can be directly applied to a real WDN, with no need for a hydraulic simulator of the plant, since both model identification and regulator tuning can be performed based on input-output data collected in situ.

CODE AVAILABILITY
The MATLAB code for the design of the MN-RTC control scheme discussed in this work is available at https://github. com/GiacomoGaluppini/Multi-Node-Real-Time-Control-of-Pressure-in-Water-Distribution-Networks-via-Model-Predictive-Control.