Real-Time Pressure Control in Water Distribution Networks: Stability Guarantees via Gain-Scheduled Internal Model Control

This article proposes a novel scheme for the real-time control (RTC) of service pressure in water distribution networks (WDNs), which is beneficial in terms of leakage reduction, energy recovery, pipe burst abatement, and extension of infrastructure lifetime. Compared with the other schemes previously proposed in the scientific literature, this novel scheme combines regulatory performance with proven guarantee of stability, which is obtained by framing gain scheduling in the context of internal model control (IMC) of linear parameter-varying (LPV) systems. Previous works relying on gain scheduling only prove stability for fixed scheduling parameter values, which is only a necessary condition for stability in case of possibly fast, time-varying parameters. The proposed RTC scheme guarantees instead stability of the closed loop for any admissible trajectory of the scheduling parameter. The novel control scheme is tested numerically against challenging operating conditions in a benchmark WDN, including two different demand patterns and four hydrant activation scenarios.


I. INTRODUCTION
I N THE context of water distribution networks (WDNs)   management, leakage reduction [1], [2], excess energy recovery [3], [4], [5], [6], pipe burst abatement [2], [7], [8], and overall infrastructure life extension represent absolute priorities.Real-time control (RTC) of service pressure [9] represents an effective approach to pursue these goals, as it allows reducing pressure excess across the WDN, while providing sufficient pressure to ensure satisfaction of users' demand.With the rise of the Water 4.0 approach [10], novel WDNs include sensors, actuators, and computing units, which are connected by wire (e.g., via optical fiber) or by high-end wireless networks [e.g., narrowband internet of things (NB-IoT)] to ensure fast and reliable communication.This technological advance enables the transition from local RTC, where the pressure is controlled right downstream of the control valve [11], [12], [13], to remote RTC, where the pressure can be controlled at any Giacomo Galuppini is with the Dipartimento di Ingegneria Industriale e dell'Informazione, University of Pavia, 27100 Pavia, Italy (e-mail: giacomo.galuppini@unipv.it).
Digital Object Identifier 10.1109/TCST.2023.3325546point of the WDN [14], [15], [16], [17].Specifically, the aim of remote RTC is regulation of service pressure at the socalled critical node, i.e., the WDN node characterized by the minimum daily pressure.As several works in the literature clearly demonstrated both in silico and in situ, a suitable choice of the desired pressure at the critical node, combined with an effective control scheme, can simultaneously reduce pressure excess and guarantee satisfaction of users' demand [13], [15], [16], [18].The majority of the literature dealing with RTC of pressure focuses on the performances of the control scheme, whereas the fundamental issues of closed-loop stability and robustness are often overlooked.Motivated by incidents occurred on real plants [19], a series of works [19], [20], [21], [22], [23] started analyzing stability and robustness of common RTC schemes and highlighted common pitfalls in the control design that can lead to instability.In particular, Galuppini et al. [22] stress the need for an accurate description of the high-frequency dynamics of the plant, in order to achieve a reliable controller design.The WDN dynamics is, in fact, characterized by a strong resonant behavior, which arises from the interaction of pressure waves traveling at finite speed through the pipes [24].A low-order model of the process, such as the first-or second-order linear models commonly adopted in the literature, may not adequately describe such resonant behavior and may, in turn, result in a poor evaluation of the stability margins of the control design [22].Moreover, Janus and Ulanicki [21] and Galuppini et al. [22] highlight the presence of strong static nonlinearities affecting the process and demonstrate how their combination with the high frequency resonance peaks can lead to closed-loop instability.Finally, Janus and Ulanicki [21] and Galuppini et al. [25] propose a combination of gain scheduling and nonlinearity inversion to robustify the control schemes and mitigate these issues.The scheduling policies aim at keeping the loop function design unaltered, as the process moves across a wide range of operating points.This guarantees that, for all fixed values of the parameters, the closed loop is (robustly) stable.The resulting gain-scheduled controllers are parametric in nature and belong to the class of linear parameter-varying (LPV) systems.It is well understood that, for LPV systems, stability for all fixed parameter values does not necessarily imply stability for time-varying parameters [26], [27].Whereas slow parameter variations can be tolerable, arbitrarily fast variations can be extremely dangerous and compromise the closed-loop stability.In the context of RTC, fast variations occur whenever the water demand increases or decreases abruptly (e.g., in case of fire hydrant operations [28] or presence of industrial users [29]).Therefore, all the aforementioned works fail to rigorously consider the issue of LPV stability and only focus on stability for fixed parameter values.The main aim of this work is to combine the findings from [21] and [25] with the rigorous gain scheduling approach discussed in [27], [30], and [31], to propose and test an RTC scheme for pressure control with LPV stability guarantees.The proposed approach is based on an LPV internal model control (IMC) scheme [27] and includes an LPV Smith predictor (SP) [32] to compensate for the presence of time delays.Low-pass filters are used to clean the scheduling signals from high-frequency oscillations [31].A filter reinitialization strategy is also proposed, to preserve LPV stability of the closed loop.
The novel control scheme is tested in a simulated environment, by relying on a detailed, pressure-driven [33] unsteady flow model of the WDN [34].This hydraulic modeling approach is the most suitable to replace the actual plant for a preliminary evaluation [35].The case study explored in this article is based on a real WDN topology and an accurate modeling of the users demand pattern [36].Simulations include the sudden opening of fire hydrants [28], to evaluate the behavior of the control system in case of fast scheduling parameter variations.All the results stress the effectiveness and the reliability of the novel RTC approach discussed in this article.
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.The definition of the nominal working point (WP) of the system is discussed in Section III-B, and the identification of the process model is presented in Section III-C.The overall control scheme is thoroughly described in Section III-D.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 analyzed in this article is the skeletonized WDN of the town of Castelfranco Emilia, Northern Italy.The WDN is composed of 27 nodes (26 demanding nodes and one source node) and 32 pipes.In addition, two fire hydrants are located at nodes 3 and 13, respectively.The complete topology is depicted in Fig. 1.For further details, including features of network nodes, pipes, and hydrants, refer to [28] and [35].
For this case study, the goal to be achieved is pressure regulation at h sp = 25 m at node 1, which is chosen as controlled node.A T.I.S. GROUP "NUOVAL" plunger valve with a diameter of 250 mm (other details and datasheet can be found in [37]) is installed in pipe 26-20, linking the source to the rest of the network.This pressure control valve (PCV) is well suited for RTC, as it can be equipped with an electric actuator and controlled by a programmable logic controller (PLC).For safety reasons, the actuation speed is limited to 0.0033 s −1 .
Residential demand patterns are generated by means of the bottom-up approach discussed in [36] and [38].The demand patterns keep into account human daily routine.With the aim of assessing the robustness of regulatory performances, two different demand patterns are considered, resulting in two different trends of the total WDN demand (see Fig. 2): a flatter trend (pattern A) and a more peaked trend (pattern B).Note that demand patterns A and B share the same average values for each single demand pattern.In addition, the presence of leakage is explicitly considered in this system.The source pressure pattern (see Fig. 2) follows the overall demand trend.

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 control algorithm introduced in this work.

A. Hydraulic Model for Simulations
The pressure-driven [33], unsteady flow modeling [34], [35] is the most suitable approach enabling an accurate 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 (2) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.] is the gravity acceleration constant, c [m/s] is the wave celerity, and J p is the friction slope.
An additional outflow per unit of pipe length q l [m 2 /s] models leakage from WDN pipes where α leak [m 2−γ /s] and γ leak [−] are the leakage coefficient and exponent, respectively.The exponent γ leak is set to 1, typical value for plastic pipes [39].Coefficient α leak [−] is set to 9.4 × 10 −9 m/s, resulting in a leakage percentage rate of 20%.Hydrant outflows q hd [m 3 /s] are obtained as the outflow from a pressurized orifice, by means of the emitter equation [40] 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 and vice versa.
Each pipe friction slope is evaluated as follows: where n [s/m ] is the Gauckler-Manning coefficient.Furthermore, pipe friction slopes are increased using the correction proposed in [41] to account for the unsteady flow effects.
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 ] is the valve cross-sectional area, ξ [−] 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 α.This function is typically available in the PCV datasheet [37] or can be characterized by laboratory experiments.Fig. 3 shows the function related to the T.I.S. GROUP "NUOVAL" plunger valve of the Castelfranco Emilia WDN.From this point onward, let Q 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 [34].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.Finally, measurement noises acting, respectively, on the measured pressure heads h(t) and on the flow at the valve site Q(t) are also included in the model.

B. WP Definition
Consider the multi-input multioutput (MIMO) system characterized by input signals as follows: 1) ξ(α(t)), the local loss coefficient, function of the valve closure ([−]); 2) H (t), the source pressure ([m]); 3) D(t), the vector of water demands ([m 3 /s]).Consider the multi-input multioutput (MIMO) system characterized by output signals as follows: 1) h(t), the measured pressure ([m]) at the critical node of the WDN; 2) Q(t) the flow at the valve site ([m 3 /s]).The local head loss coefficient ξ(α(t)) is chosen as control variable, whereas pressure h(t) is chosen as controlled variable.Source pressure H (t) and water demands D(t) are stochastic disturbances acting on the process.
For convenience, define the overall demand of the WDN D tot (t) as follows: where D i (t) is the water demand at node i and N nodes is the number of demanding nodes in the WDN.
Let the tuple WP = (ξ , H , D, h, Q) represent the WP for the MIMO system.The average values of typical H (t) and D i (t) patterns are usually available to the WDN manager (e.g., via billed consumption) and are adopted for the definition of the WP.The value of ξ must be chosen, so that the corresponding steady-state pressure at the controlled node, h, is sufficiently close to the desired one.In a simulated environment, this can be achieved by adjusting the valve closure α until the pressure h reaches the desired pressure h sp .
The value of ξ associated with such valve closure is adopted as ξ .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 [42].
Remark: The analysis carried out in previous works [19], [22], [25] underlines that the choice of ξ as control variable is fundamental to face the strong static nonlinearity affecting the system and improving the robustness of the control scheme.If α is chosen as control variable, due to the ξ(α) nonlinearity, the increase in the process gain occurring at high values of α represents the main source of instability for closed-loop pressure control systems [17], [19].However, as the ξ(α) function is typically available and invertible, this issue can be conveniently faced by selecting ξ as control variable and introducing a nonlinearity inversion block in the loop to compute at runtime the corresponding value of α [25].

C. Process Model Identification
This section summarizes the two-step procedure originally introduced in [24] and [42] for the identification of an LPV process model [26].The methodology combines the black-box identification of a local, linear model, and physical knowledge of the valve static behavior, to obtain an LPV description of the plant dynamics.
Start by considering a continuous-time, strictly proper transfer function P(s), relating local head loss variations δξ If available, a hydraulic WDN simulator (as discussed in Section III-A) can be exploited as source of input-output data.If, instead, a hydraulic model of the WDN is not available, or it is not sufficiently accurate, input-output data can be directly collected by means of experiments performed on the real plant.Note that, when working in silico, the experimental design can focus on maximizing the information contained in the dataset, but the collected data can be affected by hydraulic model versus plant mismatch.On the other hand, working in situ avoids this issue, but requires a more careful design of experiments, which should not interrupt the service to users and should not stress the WDN structure [42].Once a dataset is available, input-output identification techniques can then be applied to compute suitable parameter estimates for P(s) (e.g., via prediction error method (PEM) [43]).At this point, P(s) should provide an accurate description of the WDN dynamics in a neighborhood of WP, and a testing dataset could be effective in assessing the predictive performances of P(s).Due to the dissipative nature of the system, P(s) must be asymptotically stable.Moreover, due to the interaction of pressure waves in the WDN, P(s) is typically characterized by several high-frequency resonance peaks, possibly associated with low damping factors [22], [24].In addition, in the context of remote RTC, P(s) is also characterized by a time delay arising from the time required by pressure waves to travel the WDN pipes from the PCV to the pressure sensor, along the quickest path [24].This process model typically takes the form where µ is the static gain, T z,i and T p,i are time constants associated with real zeros and poles, ζ i and ξ i are the damping factors of complex pairs of zeros and poles, α ni and ω ni are the natural modes of complex pairs of zeros and poles, and τ is the time delay.At this point, physical knowledge of the valve can be incorporated in P(s), to extend its range of validity.In particular, note that, in the valve equation [see (6)], the pressure loss induced by the valve, H v , quadratically depends on the flow through the valve, Q.Furthermore, several works in the literature suggest that a static description of this nonlinearity can be very effective in improving the associated control design [14], [24], [25].As Q(t) can be measured online, it is possible to adopt it as scheduling parameter w(t) and move to an LPV description of the plant.Let µ be the static gain of P(s), corresponding to the flow at the valve site at WP, Q. Define P(s, w) as the LPV extension of P(s), and let its parameter-dependent static gain µ(w(t)) be Moreover, define and consider a state-space realization of P(s, w) of the form (15) (16) with u(t) = δξ(t), y(t) = δh(t), A P and C P matrices of suitable dimensions, and B P (w) a continuous function of the parameter w, which is supposed to belong to a compact set W. 1 Similarly, denote by ′ (w(t)) = {A P , B P (w(t)), C P } the state-space realization of P ′ (s, w).With the proposed class of possible state-space realizations, the system ′ (w(t)) is both asymptotically stable and LPV stable (see Appendix, Definition 1) if and only if the matrix A P is Hurwitz.

D. Control Design Methodology
This section discussed the control design methodology proposed in this article.The overall control scheme (Fig. 4) has, as fundamental building block, the IMC scheme for LPV-stable time delay systems originally developed in [27] and [30], which is adapted to meet the specific needs of this application.Each element of the control scheme is discussed in detail in the reminder of this section.
1) IMC Scheme for LPV-Stable Time Delay Systems: The time delay present in the LPV system arises from the propagation of pressure waves across the WDN, and its time scale can pose sensible limitations to the performances of the control loop [24].For LPV-stable systems as (w(t)), the effect of time delay can be effectively compensated by means    of the LPV-SP formulation proposed in [30].This approach frames the SP in the context of IMC and rearranges the control scheme to obtain a Youla-Kucera (YK) parameterization of the overall controller.Fig. 5 depicts the standard linear timeinvariant (LTI) control system based on the SP, with C(s) the transfer function of the primary controller, located in the forward path of the overall controller, R(s) (delimited by a red dashed line).Fig. 6 depicts the scheme rearranged as an IMC scheme.The negative feedback loop inside the blue dotted rectangle, whose transfer function is It must be stressed that (w) may be unstable for time varying w, even if A Q (w) is Hurwitz for all constant values of w [26], [27].To overcome this issue, it is necessary to compute an LPV-stable realization of Q(s, w), every time the scheduling parameter w varies in time.As proved in [27], every parametric transfer function Q(s, w) that is stable for all admissible constant values of w admits an LPV-stable realization (w), i.e., a realization that is stable for all functions w(t) taking values in W. Procedure 1 [27] allows the computation of such LPV-stable realization.
Procedure 1: Given w ∈ W, the following hold.1) Compute a generic realization with ÂQ (w) Hurwitz.2) Compute the positive-definite solution ϒ(w) of the following Lyapunov equation: 3) Factorize ϒ(w) as follows: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
with (w) an upper triangular matrix (Cholesky's decomposition).4) Realize Q(s, w) as follows: where Theorem 1 ensures LPV stability of the closed-loop control scheme in Fig. 7.
Theorem 1 [30]: If the system (w) is LPV stable, the control system in Fig. 7 is LPV stable (see Appendix, Definition 2) if the state-space realization of (w) is LPV stable.
The considered IMC scheme for LPV-stable time delay systems entails the design of the parametric transfer function of the controller C(s, w).In turn, this consists of two steps, namely, the design of a nominal controller transfer function C(s) based on P ′ (s) and the definition of a scheduling policy based on the online knowledge of w(t).
2) Nominal Controller Design: Several techniques are available in the literature for the design of C(s) (e.g., pole placement and loopshaping [44], [45]).
In this article, a loopshaping procedure is carried out, adopting the controller structure and the design rationale discussed in [22].As the control scheme proposed in this article includes an SP for compensating the effects of time delay, the design of the nominal primary controller transfer function C(s) is based on the delay-free transfer function P ′ (s).The controller structure consists of a proportional-integral (PI) term and an additional filter term.Specifically, let C pi (s) be the PI transfer function with K p the proportional gain, K i the integral gain, and T i = K p /K i the integral time constant.The gain of the transfer function µ C coincides with K i .Then, let C f (s) be the filter transfer function Finally, the primary controller transfer function C(s) is Note that C f (s) has unitary static gain, so that µ C coincides with the gain of the overall regulator.The integral term responds to static performances requirements (i.e., robust regulation to the constant set points, and disturbance rejection in presence of constant and slowly varying process disturbances).The filter term improves the deamplification of the high-frequency resonance peaks of P ′ (s).The two additional zeros in the controller can be used to compensate for the effect of low-frequency poles in P ′ (s) and to recover phase margin.In this article, after sufficient deamplification of resonance peaks is provided, the closed-loop bandwidth is set to the largest possible value to assess the performances of the control scheme in a demanding situation.
Let ω r p be the angular frequency associated with the resonance peak of P ′ (s) located at the lowest frequency.Then, the controller design can be carried out by setting Let T p be the time constant associated with the lowest frequency pole of P ′ (s) with null imaginary part.One can then set where ω c is the desired closed-loop bandwidth expressed in rad/s.With this procedure, ω c is the only free design parameter, which must fulfill ω c ≪ ω r p in order to provide robustness margins against gain and phase uncertainties.Note that, if 1  T f < ω c , the actual closed-loop bandwidth is smaller than the expected one, due to the effect of the pole associated with time constant T f .To conclude the nominal design phase, nominal closed-loop stability can be proved by means of Bode criterion [44].
Remark: Note that Theorem 1 ensures that the closed-loop system remains LPV stable as long as the YK parameter is stable for every fixed parameter values and is properly realized.While this work takes advantage of the control design considerations from the literature (e.g., [22], [24], [25]), the proposed control methodology would also allow for an online, trial-and-error tuning of the YK parameter transfer function [27].
3) Gain Scheduling for Static Nonlinearity Compensation: As discussed in Section III-C, if the local head loss coefficient ξ is used as control variable, the process model is characterized by a static nonlinearity depending on the flow at the valve site Q [see (10)].In order to keep loop functions constant over a wide range of operating points, the primary controller transfer function C(s) can be made parameter-dependent [25] by acting on its static gain.The static gain µ GS−Q C (w) of the parameter-dependent controller transfer function C(s, w) can be chosen as follows: 2 (28) so that the dependence on the scheduling parameter does not appear in the loop functions.
4) Gain Scheduling for Regulation Error Versus Cost of Control Trade-Off: The control approach discussed above is able to compensate the static nonlinearities of the process because of the choice of the local head loss coefficient ξ as control variable and the LPV design based on online measurements of the flow at the valve site Q.In this way, the loop functions can be kept as constant as possible over a wide range of operating points of the WDN.However, this choice may not always be the most efficient one.Consider an operating point characterized by high flow values and associated valve closure α assuming values close to zero.Under these conditions, the function ξ(α) is almost flat, and small variations of the control variable ξ require a wide variation of α to be obtained.As a result, the cost of control associated with high flow conditions is much higher than that of low flow conditions.This issue can be explicitly accounted for by leveraging the gain-scheduled nature of the control scheme, to progressively reduce the closed-loop bandwidth as the flow at the valve site Q increases.To this end, include the valve closure α in the scheduling parameter w, and note that α ∈ [0; 1], so that w = [Q α] ⊤ still belongs to a compact set.The static gain µ GS−α C (w(t)) of the new parameter-dependent primary controller transfer function, C(s, w), can be chosen according to the following scheduling law, inspired by Galuppini et al. [25]: where µ ratio (α(t)) ∈ (0, 1] is a multiplicative factor.Assuming the nominal design of C(s) to result in the maximum desired closed-loop bandwidth, µ ratio (α(t)) can be chosen as follows: with α * ∈ [0; 1], µ * ratio ∈ (0; 1], and p > 0 are design parameters. In particular, µ * ratio represents the value of µ ratio (α) when α = 0, p the power of increase of µ ratio (α) as α → α * , and α * the upper limit of the gain scheduling policy, which affects the design only for α ∈ [0; α * ].An example is reported in Fig. 8.The three parameters are selected, so that the gain of the controller is progressively reduced, as soon as the valve closure α reaches the flat region of curve ξ(α).The value of α * can, therefore, be selected by inspecting the curve ξ(α) (see Fig. 3).The values of µ * ratio and p are then selected to obtain a sufficiently fast reduction of the controller gain.Based on [25], convenient starting values are µ * ratio = 0.1 and p = 4, which should be adjusted to fit the specific case study.
5) Filtering of the Scheduling Parameter: The scheduling policies discussed above are based on steady-state considerations about the WDN and involve online measurement of flow at the valve site Q and valve closure α.Since the WDN is mainly driven by the users' time-varying demand, pressure and flow oscillations arise from both high-frequency demand oscillations and from the water hammer effect in the WDN pipes.Part of these oscillations is transferred to α(t), through the closed-loop system.Low-pass filtering of the scheduling parameter can be introduced in the control scheme to mitigate this issue.Then, let LP(s) be the transfer matrix of the low-pass filters for the scheduling parameter w(t) with T lp,Q and T lp,α the filter time constants.Denote as w lp (t) = [Q lp (t) α lp (t)] ⊤ the low-pass filtering of the scheduling parameter w(t).Note that T lp,Q can be selected offline, by analyzing the spectrum of daily records of Q(t), which can be directly collected from the plant.A straightforward choice for T lp,α could then be setting T lp,α = T lp,Q [25].
Simulations show that this low-pass filtering improves the performances of the gain-scheduled control approach.Nevertheless, this introduces a mismatch between the measured and the filtered scheduling parameter.In principle, both gain scheduling and SP design require a perfect knowledge of the parameters.In order to retain the benefits of the low-pass filtering and avoid stability issues, the following solution based on the small gain theorem [31], [45] is proposed.
1) Implement a state-space realization of the low-pass filter, such that the filter state x lp (t) coincides with the filter output w lp (t).
2) (Re)initialize the filter state and output so that every time the following condition on the loop gain is violated: where || • || is a generic operator norm, γ > 0 is a robustness threshold (accounting for measurement errors for the scheduling parameter), T (s) is the time delay operator [see (16)], and (s, w lp (t), w(t)) is given by (s, w lp (t), w(t)) = P(s, w lp (t)) − P(s, w(t)).(35) Similarly, it is possible to define filter reinitialization policies based on thresholds on the absolute difference between measured and filtered scheduling parameters.Whereas the reinitialization policy discussed earlier ensures closed-loop stability, these can improve the control performances in case of sudden scheduling parameter variations (due to, e.g., fire hydrants opening [28] and temporary activation of industrial user demands [29]).6) Input Saturation and Rate Limiter: Note that the valve closure α is upper and lower bounded.Moreover, under normal WDN operations, the valve must never be completely closed, i.e., the upper bound for α is always lower than 1.As a consequence, the head loss coefficient ξ is also upper and lower bounded by finite values.Moreover, as discussed in Section II, the valve speed is limited for safety reasons.This, in turn, imposes maximum and minimum rates of variations to ξ .In order to account for these limitations, it is possible to include models of input saturation, sat(•), and rate limiter, rate lim(•), in the control scheme, as shown in Fig. 9, and prove bounded-input bounded-output (BIBO) stability of the closed loop.Due to the saturation, u(t) is bounded.As the system (w) is quadratically stable, it is also BIBO stable, i.e., y(t) is bounded whenever u(t) is bounded.Note that the same holds for the system model inside the controller; thus, ŷ(t) is bounded whenever u(t) is bounded.The difference y(t) − ŷ(t) is also bounded.It follows that the output of the YK parameter (w) is bounded for bounded set-point signals y sp (t).In summary, y(t) is bounded for bounded set-point signals y sp (t).

IV. RESULTS
This section applies the RTC methodology proposed in this work to the Castelfranco Emilia WDN, described in Section II, and presents the results of simulations under several operating conditions.

A. Nominal Controller Design and Gain Scheduling Policies
For the Castelfranco Emilia WDN, a reasonable set point for pressure at the critical node (node 1) is h sp = 25 m.Based on the average values of typical users' demand and source pressure, the corresponding WP is Remark: As the average value of each demand D i (t) is the same for demand patterns A and B, a single WP is sufficient to cover both scenarios.
Input-output data for the identification of the nominal process transfer function, P(s), are obtained from a step response simulation around WP.The parameters of P(s) are obtained by means of PEM identification [46], by relying on MATLAB identification toolbox [43].Fig. 10 compares identification data and model prediction.The complete definition of P(s) is reported in (36), as shown at the bottom of the next page, and the Bode diagram of P(s) is reported in Fig. 11.From the Bode diagram of the modulus of P(s), it holds that ω r p = 0.1 rad/s.
The design of the nominal primary controller C(s) is based on the delay-free transfer function, P ′ (s) and is carried out following the loopshaping approach discussed in Section III-D2.
The parameters of the filter term of the controller, C f (s), are set to The required closed-loop bandwidth is ω c = 0.0306 rad/s.As the nominal gain of P ′ (s) results µ = −0.0631m, the nominal static gain of the controller is set to As P ′ (s) does not contain any real pole, T p is not defined, and the integral time constant T i can be treated as a free design parameter.In this case, T i is set to to recover some phase margin, which results φ m = 70 • .The gain margin results K m = 18.9.As for resonance peak deamplification, the analysis of Bode diagram of the loop gain function L(s) = C(s)P ′ (s) (Bode diagram in Fig. 12) shows that the highest resonance peak, located around ω r p = 0.1 rad/s, reaches −15 dB of magnitude. 2he parameters of the gain scheduling policy GS − α are The low-pass filters time constants are set to The threshold for filter reinitialization based on the small gain theorem is set to γ = 0.8.Moreover, the filter processing the flow at the valve site Q(t) is also reinitialized at time t if and the filter processing the valve closure α(t) is also reinitialised at time t if B. Closed-Loop Simulations Fig. 13(a) and (b) shows the results of whole-day, closedloop simulations for demand patterns A and B, respectively (details of results at the time scale of seconds are depicted in Fig. 13(c) and (d), respectively).In both scenarios, the control scheme provides accurate regulation of pressure to the     [25] saturation limit is reached around t = 8 h.However, because of the saturation model discussed in Section III-D6, the control loop quickly responds to the demand decrease and does not introduce any noticeable windup effect.Also note that, in both scenarios, the magnitude of high frequency oscillations of the valve closure α remains almost constant throughout the entire simulation and across very different working conditions.This highlights the benefits of the additional gain scheduling policy discussed in Section III-D4.For whole-day simulations, a quantitative assessment of the performances of the control scheme is also given in this article, by means of the following metrics (defined for discrete-time signals, with k the discrete-time instant).1) . The regulation error, accounting for the proximity of the measured pressure to the set point.I are quite similar for both demand scenarios and are aligned with the state-of-the art results in the literature (see [25], [47]), which, however, do not provide rigorous guarantees in terms of LPV stability.This further stresses the effectiveness and reliability of the approach proposed in this article.
As already discussed in this article, one of the main advantages of the proposed control scheme is its ability to face arbitrarily fast changes in the scheduling parameter.Specifically, as the flow Q appears in the scheduling parameter, this issue can occur in case of pipe bursts or fire hydrant operations.Four simulations involving fire hydrant operations are then carried out in order to assess the reliability of the control scheme under these circumstances.In particular, Fig. 13(e) and (g) shows the results of fire hydrant operations at nodes 3 and 13, respectively, for demand pattern A. Fig. 13(f Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.remain open for about 1600 s, from t = 0.83 h to t = 1.27 h.In all four cases, pressure is promptly regulated to the set point within few minutes from the fire hydrant opening or closure, with negligible overshoot in the response (recall that the usual sources of disturbance are always active).The dynamic response is consistent with the expected one, based on the loop function design discussed earlier in this section.Interestingly, the sudden variation in both flow at the valve site and valve closure results in low-pass filter reinitialization [see Fig. 13(e), (f), and (h)] due to (43) and (44).On the contrary, the stability-preserving reinitialization was never triggered in this simulations.

C. Comparison With the State of the Art
As discussed earlier in this article, several gain-scheduled control schemes are available in the literature of pressure RTC in WDNs (e.g., [21], [25], [35]).However, none of them rigorously considers the issue of closed-loop stability for any admissible time history of the scheduling parameters.This is a fundamental step toward the implementation of pressure RTC schemes in real WDNs.Loss of stability may trigger wide pressure oscillations that can stress and possibly damage the infrastructure and disrupt the service to end users [19], [22].The control scheme proposed in this work provides stability guarantees, while retaining the advantages of [25], i.e., a straightforward design of the nominal regulator, and simple gain-scheduling policies.A comparison with the results achieved by the filtered proportional-integral with smith predictor and gain scheduling (FPI-SP-gs) control scheme from [25] highlights that the novel control scheme can result in similar regulatory performances, while guaranteeing closedloop stability.To ensure a meaningful comparison, simulations are carried out by requiring the same nominal loop functions and using the same gain-scheduling policies.The results of the FPI-SP-gs control scheme are reported in Table I.A comparison with the results achieved by the LPV IMC scheme stresses that the performances are almost identical for demand profile A, as expected from the same design of the loop functions and similar scheduling policies.With demand profile B, the performances of the two control schemes are still similar, but not identical.Further investigations highlights that the main differences occur while the control action is saturating.This is consistent with a different handling of saturation in the two schemes.

V. DISCUSSION
The closed-loop simulations presented in the previous section highlight the effectiveness of the RTC scheme proposed in this work and support its implementation on a real plant.Note that the definition of the nominal WP, as well as the identification of the LPV process model, can be both performed directly on the real WDN with no need for an accurate hydraulic model, by following the guidelines discussed in [42].In this case, the knowledge of average users' demand would not even be required.
Moreover, as previously introduced, the proposed control scheme enables online tuning of the YK parameter Q(s), with no risk for the stability of the control loop.This would enable a safe application of reinforcement learning methodologies for the design of the controller.In fact, it is usually hard to directly relate relevant features of the loop functions and in situ performance metrics, as the dynamics of the main disturbances affecting the plant (i.e., users' demands) are difficult to model, especially at a fine timescale.For these reasons, Galuppini et al. [47] proposed an alternative approach for the design of the controller, based on biobjective optimization.In the proposed formulation, the cost functions to be minimized/maximized encode the conflicting requirements of control performances (e.g., closed-loop settling time and process disturbance rejection) and moderation of the control action/rejection of measurement noise/robustness to model uncertainties.An integral action is still included in the controlled to ensure robust static performances.This approach provides a Pareto front of tunings to be tested on the real plant and evaluated according to the desired performance metrics.The possibility of changing the controller parameters online, without compromising the stability of the closed loop, represents an interesting feature of this work, which well complements the biobjective optimization approach, as it enables an automated choice of the "best" tuning over the available candidates forming the Pareto front.Therefore, the development of a fully automated controller design procedure represents an interesting research direction.Another possible extension of this work could focus on the identification of a full LPV model from experimental data.In fact, while this article only considers a static parametrization of the process model, the introduction of parameter-dependent dynamics and time delay may improve the predictive performances of the model and, possibly, of the overall control scheme.Note that the control scheme proposed in this article would not require any modification, as long as the new process model is LPV stable.Finally, run-to-run [48] or Bayesian optimization [49] approaches could be leveraged to automate and optimize the design parameters of the gain-scheduling policies as well as the time constants of the low-pass filters.

VI. CONCLUSION
This article deals with RTC of service pressure in WDNs and proposes a gain scheduling approach guaranteeing stability of the closed-loop system for arbitrary variations of the scheduling parameter.The proposed approach combines an IMC scheme with gain scheduling, in order to handle an LPV description of the WDN dynamics.The control scheme also includes an SP, to compensate for time delay effects, as well as models of the static nonlinearities affecting the control action, to avoid windup issues.Finally, to improve the performances, low-pass filtering of the scheduling parameter is also implemented, and filter reinitialization policies are included in the scheme to prevent loss of stability due to measured/filtered scheduling parameter mismatch.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 even in case of very fast variations of the scheduling parameter.

Fig. 3 .
Fig. 3. Local head loss coefficient ξ as a function of the valve closure α.

Fig. 4 .
Fig. 4. Block scheme for the proposed RTC control algorithm, including parameter-varying internal model controller, saturation sat(•), and rate limiter rate lim(•) of control action, parameter-varying SP, gain scheduling policies, and low-pass filters.

Fig. 5 .
Fig. 5. LTI control system for a time delay plant based on the SP.
) represents the YK parameter of a controller whose positive feedback path (outside the blue dotted rectangle) contains a model of the time delay system.In case of LPV systems, transfer functions should be replaced by state-space realizations, as depicted in Fig. 7, to correctly address the issue of LPV stability.Let (w) = {A Q (w), B Q (w), C Q (w), D Q (w)} be the state-space realization of the proper LPV YK parameter Q(s, w), where A Q (w), B Q (w), C Q (w), and D Q (w) are continuous functions of w ∈ W.

Fig. 8 .
Fig. 8. Example of gain scheduling policy for regulation performances versus cost of control trade-off.

Fig. 11 .Fig. 12 .
Fig. 11.Bode diagram of the nominal process transfer function chosen set point, by rejecting the disturbances generated by time-varying users' demand and source pressure.For demand pattern B [Fig.13(b)], due to its first, high peak, the valve |α(k) − α(k − 1)| [−].The control effort, accounting for energy consumption and wear of actuators.All signals are sampled with a 1-s sampling time; for wholeday simulations, then K tot = 86400.The values reported in Table

Fig. 13 .
Fig. 13.Closed-loop simulations.Top: pressure head at node 1 h 1 (t) (blue solid line) and pressure set point h sp (red dashed line).Middle: flow at the valve site Q(t) (red solid line) and low-pass-filtered flow at the valve site Q lp (t) (green dashed line).Bottom: valve closure α(t) (black solid line) and low-pass-filtered valve closure α lp (t) (cyan dashed line).(a) Whole day, demand pattern A. (b) Whole day, demand pattern B. (c) Detail of whole-day simulation, demand pattern A. (d) Detail of whole-day simulation, demand pattern B. (e) Hydrant opening at node 3, demand pattern A. (f) Hydrant opening at node 3, demand pattern B. (g) Hydrant opening at node 13, demand pattern A. (h) Hydrant opening at node 13, demand pattern B.

TABLE I RESULTS
OF CLOSED-LOOP SIMULATIONS WITH DEMAND PROFILES A AND B FOR THE LPV IMC CONTROL SCHEME PROPOSED IN THIS ARTICLE AND FOR THE FPI-SP-GS CONTROL SCHEME FROM