*Article* **Predictive Control in Water Distribution Systems for Leak Reduction and Pressure Management via a Pressure Reducing Valve**

**Jose-Roberto Bermúdez 1, Francisco-Ronay López-Estrada 1, Gildas Besançon 2, Guillermo Valencia-Palomo 3,***<sup>∗</sup>* **and Ildeberto Santos-Ruiz <sup>1</sup>**


**Abstract:** This work proposes a model predictive control (MPC) strategy for pressure management and leakage reduction in a water distribution system (WDS). Unlike most of the reported models that mainly consider EPANET-based models, the proposed method considers its dynamic representation given by ordinary differential equations. The proposed MPC uses a pressure-reducing valve (PRV) as a control element to regulate the pressure in the WDS to track the demand. The control scheme proposes a strategy to manage the high nonlinearity of the PRV and takes into account the demand profile throughout the day as well as the leaks that occur in the pipeline. The estimates of magnitude and location of the leak are provided by an Extended Kalman Filter from previous work and with the aid of a rule-based set point manager reduces the fluid loss in the event of a leak. Different scenarios are studied to illustrate the effectiveness of the proposed control system, achieving an approximate reduction of up to 5% of water losses, demonstrating robustness in the case of uncertainty in the leak location estimate.

**Keywords:** pipelines; control valves; leak reduction; water distribution system; pressure management

#### **1. Introduction**

Water distribution systems (WDS) are the most sustainable and efficient means of transporting fluids such as drinking water, natural gas, and oil [1]. These systems are composed of pipelines, pipe joints, connection nodes, and other components such as valves and pumps, which are prone to damage due to aging or unwanted events, such as earthquakes, floods, and lack of management, among others [2]. According to a study by the Organisation for Economic Co-operation and Development [3], these abnormal events cause water losses by leaks that reach almost 21%. The study examined water usage in 48 major cities across 17 countries, finding that in some cities of Mexico, the percentage was more than 40%. This highlights the need to propose technological developments to mitigate such losses.

From the control theory point of view, it is essential to investigate techniques to reduce the leaks without affecting the demand. In the literature, it has been found that the most effective approach to find a trade-off between maintaining the desired demand and reducing the water losses is by considering pressure water controllers on critical nodes of the water distribution network [4]. However, many challenges are associated with increasing demand and managing the pressure levels. In particular, pressure-reducing valves (PRV) are the recommended actuators to minimize these undesirable effects and operate the WDS effectively by following a pattern of demands [5].

**Citation:** Bermúdez, J.-R.;

López-Estrada, F.-R.; Besan¸son, G.; Valencia-Palomo, G.; Santos-Ruiz, I. Predictive Control in Water Distribution Systems for Leak Reduction and Pressure Management via a Pressure Reducing Valve. *Processes* **2022**, *10*, 1355. https:// doi.org/10.3390/pr10071355

Academic Editor: Alfredo Iranzo

Received: 1 July 2022 Accepted: 8 July 2022 Published: 12 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Modeling the WDS is vital for designing any control algorithm [6,7]. These models are analyzed taking into account the location and number of PRVs for a correct evaluation of the pressure [8]. For instance, Mazumder et al. [9] developed an optimization method based on genetic algorithms for pressure management in a WDS by adjusting a PRV of a hydraulic network designed in EPANET. Parra et al. [10] proposed a pressure management system in EPANET composed of a PRV and a pump as a turbine where the hydraulic model played a decisive role identifying critical nodes and predicting hydraulics properties in the network. García-Ávila et al. [11] designed a water leakage minimization system by optimizing pressure using a PRV. The hydraulic model of the WDN was developed using EPANET and WaterNetGen software. Dini and Asadi [12] designed a methodology based on particle swarm optimization (PSO) to identify the PRV that requires adjustment and obtain pressure management in the system; the network model was designed in EPANET. Hernández et al. [13] proposed a detrended fluctuation analysis (or DFA) to highlight some of the traits such as the head loss of high-viscosity gas–liquid flows. Navarro et al. [14] designed a leak diagnosis system for a pipeline and residue analysis using genetic algorithms (GA) to minimize location error. Jara-Arriagada and Stoianov [15] designed a sensitivity analysis system to evaluate the potential impact of pressure control in a WDS to reduce pipe breaks by applying a logistic regression technique. Mosetlhe et al. [16] proposed a review of techniques for locating PRVs and controlling pressure in a WDS, minimizing problems such as excessive pressure in the system. Mathye et al. [17] designed a pressure management system to reduce leaks through PRVs, taking into account consumption and leak flow. All those works developed algorithms for pressure management using hydraulic models based on EPANET and waterNetGen software, whose main limitation is that the models only work in steady-state behavior, while the effects of leaks and pressure changes are dynamic.

The present work proposes a dynamical approach to pressure management in a WDS. The transient effects due to pressure changes or leaks are modeled on the basis of waterhammer equations. Then, a constrained model predictive control (MPC) system is proposed to track the desired pressure profile driven by a set point manager that considers the water loss due to leaks and the demanded pressure profile. Moreover, a strategy to handle the high nonlinearity of the PRV control input is proposed. Finally, some simulations are proposed by considering the mathematical model of a real distribution system that can be configured as a single pipeline or a branched system. The results illustrate the effectiveness of the proposed method in the presence of physical constraints, noise, and transient behaviors due to leaks, saving up to 5% of water losses in the event of a leak and demonstrating robustness in the case of uncertainty in the leak location estimate. The rest of the document is organized as follows: Section 2 presents the considered case studies; Section 3 describes the problem formulation; Section 4 is devoted to the control strategy; Section 5 presents the simulation results; finally, Section 6 draws the conclusions.

#### **2. Case Studies**

Three different cases will be considered in this paper, corresponding to different hydraulic system structures with the pressure to be controlled at a specific position (called controlled node) and some leaks affecting the system at some other position.

#### *2.1. Case 1: Pipeline with a Leak before the Controlled Node*

In this case, the system is a pipeline under a leak like the one shown in Figure 1. It is composed of a reservoir that provides the fluid to the pipe divided into three sections for convenience. The first section is related to the distance between the inlet node and the leak (*z*-); the second one is related to the distance from the leak to the controlled node (*z*<sup>2</sup> = *z*<sup>1</sup> − *z*-); and finally, the third section is related to the distance between the controlled node and the outlet node (*z*3).

**Figure 1.** Case 1: Pipeline with a leak before the controlled node.

The mathematical model can be determined by considering water-hammer equations [18], which can be approximated on the basis of the three considered sections as follows (e.g., [19,20]):

$$\dot{Q}\_1 = \frac{\mathfrak{a}\_1}{z\_\ell} (H\_1 - H\_\ell) + \mu\_1 (Q\_1) Q\_1 |Q\_1| - \frac{\mathfrak{a}\_1}{z\_1} \Delta H\_{\mathfrak{v}\_\ell} \tag{1}$$

$$
\dot{H}\_{\ell} = \frac{\mathfrak{a}\_2}{z\_{\ell}} (\mathcal{Q}\_1 - \mathcal{Q}\_2 - \mathcal{Q}\_{\ell}),
\tag{2}
$$

$$\dot{Q}\_2 = \frac{\mu\_1}{z\_2} (H\_\ell - H\_3) + \mu\_2 (Q\_2) Q\_2 |Q\_2| \,\text{\textquotedbl{}Q-Dicke-V} \tag{3}$$

$$
\dot{H} \mathfrak{J} = \frac{\mathfrak{a}\_2}{z\_2} (Q\mathfrak{z} - Q\mathfrak{z}),
\tag{4}
$$

$$Q\_3 = \frac{\alpha\_1}{z\_3}(H\_3 - H\_4) + \mu\_3(Q\_3)Q\_3|Q\_3|\,\tag{5}$$

with

$$a\_1 = \lg A\_{r\prime} \quad a\_2 = \frac{c^2}{\lg A\_{r\prime}} \quad Q\_\ell = \lambda \sqrt{|H\_\ell|}$$

$$\mu\_i(Q\_i) = \frac{-f(Q\_i)}{2DA\_r} , \quad f(Q\_i) = \frac{1.325}{\left[\ln\left(\frac{c}{3.7d} + \frac{5.74}{(\frac{4Q\_i}{ndv})^{0.9}}\right)\right]^2} ^\prime$$

where *H*1, *H*-, *H*3, and *H*<sup>4</sup> are the piezometric heads (m) at the inlet, leak node, controlled node, and outlet, respectively; *Q*1, *Q*2, *Q*<sup>3</sup> are the volumetric flow rates in each section (m3/s); *g* is the gravitational constant (m/s2); *Ar* is the cross-sectional area of the pipe (m2); *c* is the wave speed (m/s); *d* is the pipeline diameter (m); *ν* represents the kinematic viscosity; the friction term is *f*(*Qi*); is the roughness of the pipe; and *λ* is the leak coefficient. Finally, Δ*Hv* describes the PRV effect.

A PRV is an actuator used to reduce downstream pressure. Internally, a PRV is made of a fixed orifice, a pilot valve, and a needle valve [21]. The shutter is the outer mechanism for adjusting the outlet pressure, either increasing or decreasing the pressure in a range from 0 to 100%. Figure 2 shows the schematic view of a PRV, where *Hv*<sup>1</sup> and *Hv*<sup>2</sup> are the piezometric head at valve ends, *r* ∈ (0, 1] is the valve adjustment (control input), and *Q*<sup>1</sup> is the flow through the valve.

The differential pressure in a PRV is described as [22]:

$$
\Delta H\_{\text{\textquotedblleft}1} = H\_{\text{\textquotedblleft}1} - H\_{\text{\textquotedblright}2} = \frac{Q\_1 |Q\_1|}{(rE)^2},\tag{6}
$$

where *E* = *CvAv*(2*g*)1/2 is the Torricelli expression, *Cv* is the discharge coefficient of the valve, and *Av* is the cross-sectional area of the valve (m2).

**Figure 2.** Schematic view of a PRV.

#### *2.2. Case 2: Pipeline with a Leak after the Controlled Node*

When the leak appears after the controlled node, as shown in Figure 3, then *z*<sup>2</sup> = *z*<sup>3</sup> − *z*-, and the model becomes

$$\dot{Q}\_1 = \frac{\mathfrak{a}\_1}{\mathfrak{z}\_1} (H\_1 - H\_2) + \mu\_1 (Q\_1) Q\_1 |Q\_1| - \frac{\mathfrak{a}\_1}{\mathfrak{z}\_1} \Delta H\_{\mathfrak{v}\prime} \tag{7}$$

$$
\dot{H}\_2 = \frac{\mu\_2}{z\_1} (Q\_1 - Q\_2),
\tag{8}
$$

$$
\dot{Q}\_2 = \frac{\mu\_1}{z\_\ell} (H\_2 - H\_3) + \mu\_2 (Q\_2) Q\_2 |Q\_2| \,\tag{9}
$$

$$
\dot{H}\_{\ell} = \frac{\mathfrak{a}\_2}{z\_{\ell}} (\mathcal{Q}\_2 - \mathcal{Q}\_3 - \lambda \sqrt{H\_{\ell}}),
\tag{10}
$$

$$
\dot{Q}\_3 = \frac{\alpha\_1}{z\_2} (H\_3 - H\_4) + \mu\_3 (Q\_3) Q\_3 |Q\_3| \,. \tag{11}
$$

**Figure 3.** Case 2: Pipeline with a leak after the controlled node.

#### *2.3. Case 3: Branched Water Distribution Network*

In this third case, a branched water distribution network is considered with the topology shown in Figure 4. With the same notations as before, the mathematical model describing branch flows *Qi* and node pressures *Hi* can be given by:

$$
\dot{Q}\_1(t) = \frac{\alpha\_1}{z\_1} (H\_{\rm B1} - H\_2) + \mu\_1 (Q\_1) Q\_1 |Q\_1| - \Delta H\_{\rm v} \tag{12}
$$

$$
\dot{H}\_2(t) = \frac{a\_2}{z\_1} (Q\_1 - Q\_2 - Q\_4)\_\prime \tag{13}
$$

$$Q\_2(t) = \frac{\alpha\_1}{z\_2} (H\_2 - H\_\ell) + \mu\_2 (Q\_2) Q\_2 |Q\_2| \tag{14}$$

$$
\hat{H}\_{\ell}(\mathbf{t}) = \frac{\alpha\_2}{z\_2} (Q\_2 - Q\_5 - Q\_3 - Q\_\ell) \,\tag{15}
$$

$$
\hat{Q}\_3(t) = \frac{\mathfrak{a}\_1}{\mathfrak{z}\_3} (H\_3 - H\_{\mathfrak{B}2}) + \mu\_3(Q\_3) Q\_3 |Q\_3| \,\tag{16}
$$

$$
\dot{Q}\_4(t) = \frac{\kappa\_1}{z\_4} (H\_2 - H\_{B3}) + \mu\_4 (Q\_4) Q\_4 |Q\_4| \,\tag{17}
$$

$$
\dot{Q}\_5(t) = \frac{\varkappa\_1}{z\_5} (H\_3 - H\_{B4}) + \mu\_5(Q\_5) Q\_5 |Q\_5| . \tag{18}
$$

The mathematical models developed here are valid for any pipeline system with the configurations described in Figures 1–4. For validation tests, the physical parameters were taken from a real system located at the Hydroinformatics Laboratory of the Technological Institute of Tuxtla Gutiérrez, whose mathematical model was presented in [23].

**Figure 4.** Case 3: Branched water distribution network.

**Remark 1.** *Notice that λ and z are unknown parameters that must be estimated by considering a leak detection and localization method. For Cases 1 and 2, we consider our previous result published in [24] where a leak location and estimation method was designed using an extended Kalman filter (EKF). For Case 3, an EKF as the one in [25] can be used. Therefore, for the sake of simplicity, λ and zare assumed to be known.*

#### **3. Problem Formulation**

WDS are designed to meet the desired demands at their ends, even when affected by leaks. It is important to note that these demands have different profiles depending on the time and day. Typically, the demand is higher during the day and lower during the early morning [26]. In case of a leak event, the fluid loss rate (leak magnitude) can be reduced by reducing the pressure on the controlled node. However, this also reduces the flow at the node of demand, compromising the main objective of the WDS. In this regard, the control strategy for reducing leaks must consider a time-varying profile, with the primary objective to maintain a trade-off between reducing the leak magnitude and maintaining the desired demand.

To address this problem, the control scheme shown in Figure 5 is proposed. This scheme is made of three components: an extended Kalman filter (EKF), a pressure controller, and a set point manager. The EKF is used to estimate the flows and pressures along the WDS and to detect and estimate the leak position and its magnitude by using only pressure head and flow rate measurements at the pipeline ends. The EKF considered for Cases 1 and 2 is the one reported in our previous work in [24] and, for Case 3, the one in [27] is considered; however, any other leak location and estimation method can be used for the proposed scheme. The pressure controller is an MPC that takes into account physical constraints and the leak dynamics. For its operation, the MPC uses the estimated values obtained from the EKF (*z*ˆ-, *λ*ˆ , *H*ˆ -) and measured flows and pressures from the WDS. It is important to mention that the MPC also considers the PRV model whose behavior is highly nonlinear due to the inverse quadratic term of its control variable (*r*), which represents a challenge for the control system. Finally, the set point manager block provides the reference

pressure (*sk*) to the MPC; it handles the trade-off between fluid loss and fulfilment of a demanded pressure profile (DP) in the event of a leak.

**Figure 5.** Block diagram representation of MPC used in WDS.

#### **4. Pressure Control**

The control strategy adopted to regulate the pressure in the WDS is an MPC. Figure 6 shows the implemented scheme which will be detailed in this section. The basic idea in MPC is to calculate the control action at each sampling instant through the solution of an optimization problem, which is written in terms of a prediction model.

**Figure 6.** MPC scheme for pressure control.

#### *4.1. PRV Handling*

The PRV has the function of regulating the incoming water for a safer constant predetermined downstream level. The control signal *r* establishes the downstream pressure. One of the main challenges in using a PRV is its high nonlinearity as an inverse quadratic term (*r*) represents its control variable. This term affects the control input considerably because the more the PRV is closed, the greater the effects on the pressure. Then, to use linear models in the MPC formulation, the PRV control variable *r* is substituted in (6) by a virtual input *uv* to artificially hide the nonlinear behavior of the PRV for the prediction model. This is accomplished with:

$$
\Delta H\_{\text{\textquotedblleft}} = \frac{Q\_1 |Q\_1| \mu\_{\text{\textquotedblleft}}}{E^2}, \quad \text{with} \quad \mu\_{\text{\textquotedblleft}} = \frac{1}{r^{2}}.\tag{19}
$$

where *uv* ∈ [1, ∞), but for practical reasons *uv* ranges from 1 (for the PRV fully open, *r* = 1) to 100 (for the PRV 90% closed, *r* = 0.1). This new term is substituted in the three models discussed in Section 2. In this way, the MPC calculates *uv* instead of *r*. Next, *uv* is transformed back to *r* before it is applied to the WDS with *u* = *r* = <sup>√</sup><sup>1</sup> *uv* , as shown in Figure 6.

#### *4.2. Prediction Models*

The standard MPC strategy considers linear prediction models [28,29], so the models of the case studies listed before need to be linearized.

For a pipeline with a leak before the controlled node, the state and input vectors are defined as *x*¯ = [*x*1, *x*2, *x*3, *x*4, *x*5]=[*Q*1, *H*-, *Q*2, *H*3, *Q*3], and *u*¯ = [*H*1, *H*4, *λ*, *uv*] *T*. The system (1)–(5) is linearized at the operating point (*x*¯∗, *u*¯∗) taking the form *x*¯˙(*t*) = *Acx*¯(*t*) + *Bcu*¯(*t*) + *δc*, where *δ<sup>c</sup>* is the offset caused by linearization. The Jacobian matrices are given by:

$$\begin{aligned} A\_{\mathcal{L}} &= \begin{bmatrix} a\_{11} & -\frac{a\_{1}}{z\_{\ell}} & 0 & 0 & 0\\ \frac{a\_{2}}{z\_{\ell}} & 0 & -\frac{a\_{2}}{z\_{\ell}} & 0 & 0\\ 0 & \frac{a\_{1}}{z\_{\ell}} & a\_{33} & -\frac{a\_{1}}{z\_{2}} & 0\\ 0 & 0 & \frac{a\_{2}}{z\_{2}} & 0 & -\frac{a\_{2}}{z\_{2}}\\ 0 & 0 & 0 & \frac{a\_{1}}{z\_{3}} & a\_{55} \end{bmatrix},\\ B\_{\mathcal{L}} &= \begin{bmatrix} \frac{a\_{1}}{z\_{\ell}} & 0 & 0 & \frac{-gA\_{\ell}x\_{1}^{\prime}|x\_{1}^{\prime}|}{E\_{z}^{2}}\\ 0 & 0 & -\frac{-a\_{2}\sqrt{H\_{2}^{\prime}}}{z\_{2}} & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & \frac{-a\_{1}}{z\_{3}} & 0 & 0 \end{bmatrix},\end{aligned}$$

with

$$\begin{split} a\_{11} &= 2|x\_1^\*|\mu(x\_1^\*) + \frac{2(1.325)}{0.9\frac{5.74}{\left(\frac{4}{\pi d}\right)^{0.9}}} \frac{(x\_1^\*)^{1.9}}{\ln\left(\left(\frac{\epsilon}{3.74}\right) + \frac{5.74}{\left(\frac{4\pi}{\pi d}\right)^{0.9}}\right)^3} \\ &- \frac{2A\_{r\mathcal{G}}|x\_1^\*|\mu\_v^\*}{z\_1 E^2}, \\ a\_{33} &= 2|x\_3^\*|\mu(x\_3^\*) + \frac{2(1.325)}{0.9\frac{5.74}{\left(\frac{4}{\pi d}\right)^{0.9}}} \frac{(x\_3^\*)^{1.9}}{\ln\left(\left(\frac{\epsilon}{3.74}\right) + \frac{5.74}{\left(\frac{4\pi}{\pi d}\right)^{0.9}}\right)^3} \\ a\_{55} &= 2|x\_5^\*|\mu(x\_5^\*) + \frac{2(1.135)}{0.9\frac{5.74}{\left(\frac{4}{\pi d}\right)^{0.9}}} \frac{(x\_5^\*)^{1.9}}{\ln\left(\left(\frac{\epsilon}{3.74}\right) + \frac{5.74}{\left(\frac{4\pi}{\pi d}\right)^{0.9}}\right)^3} \\ \delta\_{\hat{c}} &= f\_1(\vec{x}^\*, \vec{\mu}^\*) - A\_{\hat{c}}\vec{x}^\* - B\_{\hat{c}}\vec{\mu}^\* + \Delta\_{1/2} \end{split}$$

where *f*1(*x*¯∗, *u*¯∗) is the nonlinear model of Case 1 evaluated at the linearization point, and Δ<sup>1</sup> gathers terms of order larger than 1.

For a pipeline with a leak after the controlled node, the state and input vectors are defined as *x*¯ = [*x*1, *x*2, *x*3, *x*4, *x*5]=[*Q*1, *H*2, *Q*<sup>2</sup> *H*-, *Q*3], *u*¯ = [*H*1, *H*4, *λ*, *uv*] *<sup>T</sup>*. Using *a*11, *a*33, *a*<sup>55</sup> from the previous case, the Jacobian matrices for system (7)–(11) are:

$$A\_{\mathcal{C}} = \begin{bmatrix} a\_{11} & -\frac{\underline{\alpha\_1}}{z\_1} & 0 & 0 & 0\\ \frac{\underline{\alpha\_2}}{z\_1} & 0 & -\frac{\underline{\alpha\_2}}{z\_1} & 0 & 0\\ 0 & \frac{\underline{\alpha\_1}}{z\_\ell} & a\_{33} & -\frac{\underline{\alpha\_1}}{z\_\ell} & 0\\ 0 & 0 & \frac{\underline{\alpha\_2}}{z\_\ell} & 0 & \frac{\underline{\alpha\_2}}{z\_\ell}\\ 0 & 0 & 0 & \frac{\underline{\alpha\_1}}{z\_2} & a\_{55} \end{bmatrix} \prime$$

$$B\_{\mathfrak{C}} = \begin{bmatrix} \frac{a\_1}{z\_1} & 0 & 0 & \frac{-\overline{\chi}A\_r x\_1^\* |x\_1^\*|}{E\_2^2 z\_1} \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & -\frac{a\_2 \sqrt{H\_4^\*}}{z\_\ell} & 0 \\ 0 & \frac{-a\_1}{z\_2} & 0 & 0 \end{bmatrix}.$$

,

$$\delta\_{\mathfrak{c}} = f\_2(\mathfrak{x}^\*, \mathfrak{u}^\*) - A\_{\mathfrak{c}} \mathfrak{x}^\* - B\_{\mathfrak{c}} \mathfrak{u}^\* + \Delta\_{\mathfrak{2}'}$$

where *f*2(*x*¯∗, *u*¯∗) and Δ<sup>2</sup> refer to the nonlinear model of Case 2 as before.

For the branched water distribution network, the states and the inputs are defined as *x*¯ = [*x*1, *x*2, *x*3, *x*4, *x*5, *x*6, *x*7]=[*Q*1, *H*2, *Q*2, *H*-, *Q*3, *Q*4, *Q*5], *u*¯ = [*HB*1, *HB*2, *HB*3, *HB*4, *λ*, *uv*] *T*. The Jacobian matrices for system (12)–(18) are:

$$A\_c = \begin{bmatrix} a\_{11} & -\frac{a\_1}{z\_1} & 0 & 0 & 0 & 0 & 0\\ \frac{a\_2}{z\_1} & 0 & -\frac{a\_2}{z\_2} & 0 & 0 & -\frac{a\_2}{z\_1} & 0\\ 0 & \frac{a\_1}{z\_2} & a\_{33} & -\frac{a\_1}{z\_2} & 0 & 0 & 0\\ 0 & 0 & \frac{a\_2}{z\_2} & 0 & -\frac{a\_2}{z\_2} & 0 & -\frac{a\_2}{z\_2}\\ 0 & 0 & 0 & \frac{a\_1}{z\_3} & a\_{55} & 0 & 0\\ 0 & \frac{a\_1}{z\_4} & 0 & 0 & 0 & a\_{66} & 0\\ 0 & 0 & 0 & \frac{a\_1}{z\_5} & 0 & 0 & a\_{77} \end{bmatrix} \prime$$

$$B\_{\mathcal{C}} = \begin{bmatrix} \frac{a\_1}{z\_1} & 0 & 0 & 0 & 0 & \frac{-\sqrt{A\_1 x\_1^\*} \left| \left| x\_1^\* \right|}{E\_2^\* z\_1} \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{-a\_2 \sqrt{H\_3^\*}}{z\_1} & 0 \\ 0 & \frac{-a\_1}{z\_3} & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{-a\_1}{z\_4} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{-a\_1}{z\_5} & 0 & 0 \end{bmatrix},$$

with

*a*<sup>11</sup> =2|*x*<sup>∗</sup> <sup>1</sup> |*μ*(*x*<sup>∗</sup> <sup>1</sup> ) + <sup>2</sup>(1.135) 0.9 5.74 ( <sup>4</sup> *<sup>π</sup>d<sup>ν</sup>* )0.9 (*x*∗ <sup>1</sup> )1.9 ln& ( 3.7*<sup>d</sup>* ) + 5.74 ( 4*x*∗ <sup>1</sup> *<sup>π</sup>d<sup>ν</sup>* )0.9 '<sup>3</sup> <sup>−</sup> <sup>2</sup>*Arg*|*x*<sup>∗</sup> <sup>1</sup> |*u*<sup>∗</sup> *v <sup>z</sup>*1*E*<sup>2</sup> , *a*<sup>33</sup> =2|*x*<sup>∗</sup> <sup>3</sup> |*μ*(*x*<sup>∗</sup> <sup>3</sup> ) + <sup>2</sup>(1.135) 0.9 5.74 ( <sup>4</sup> *<sup>π</sup>d<sup>ν</sup>* )0.9 (*x*∗ <sup>3</sup> )1.9 ln& ( 3.7*<sup>d</sup>* ) + 5.74 ( 4*x*∗ <sup>3</sup> *<sup>π</sup>d<sup>ν</sup>* )0.9 '<sup>3</sup> , *a*<sup>55</sup> =2|*x*<sup>∗</sup> <sup>5</sup> |*μ*(*x*<sup>∗</sup> <sup>5</sup> ) + <sup>2</sup>(1.135) 0.9 5.74 ( <sup>4</sup> *<sup>π</sup>d<sup>ν</sup>* )0.9 (*x*∗ <sup>5</sup> )1.9 ln& ( 3.7*<sup>d</sup>* ) + 5.74 ( 4*x*∗ <sup>5</sup> *<sup>π</sup>d<sup>ν</sup>* )0.9 '<sup>3</sup> , *a*<sup>66</sup> =2|*x*<sup>∗</sup> <sup>6</sup> |*μ*(*x*<sup>∗</sup> <sup>6</sup> ) + <sup>2</sup>(1.135) 0.9 5.74 ( <sup>4</sup> *<sup>π</sup>d<sup>ν</sup>* )0.9 (*x*∗ <sup>6</sup> )1.9 ln& ( 3.7*<sup>d</sup>* ) + 5.74 ( 4*x*∗ <sup>6</sup> *<sup>π</sup>d<sup>ν</sup>* )0.9 '<sup>3</sup> ,

$$\begin{split} a\_{77} &= 2|\boldsymbol{\chi\_{3}^{\*}}|\mu(\boldsymbol{\chi\_{7}^{\*}}) + \frac{2(1.135)}{0.9\frac{5.74}{(\frac{4}{\pi d \boldsymbol{\nu}})^{0.9}}} \frac{(\boldsymbol{\chi\_{7}^{\*}})^{1.9}}{\ln\left(\frac{\boldsymbol{\epsilon}}{(\frac{\boldsymbol{\epsilon}}{3.74})} + \frac{5.74}{(\frac{4\boldsymbol{\epsilon}\_{7}^{\*}}{\pi d \boldsymbol{\nu}})^{0.9}}\right)^{3/4}} \\ \boldsymbol{\delta\_{\boldsymbol{\varepsilon}} &= f\_{3}(\boldsymbol{\tilde{x}^{\*}}, \boldsymbol{\tilde{u}^{\*}}) - A\_{\boldsymbol{\varepsilon}}\boldsymbol{\tilde{x}^{\*}} - B\_{\boldsymbol{\varepsilon}}\boldsymbol{\tilde{u}^{\*}} + \Delta\_{3\prime} \end{split}$$

where *f*3(*x*¯∗, *u*¯∗) and Δ<sup>3</sup> refer to the nonlinear model of Case 3.

For all cases, the controlled output is *y*(*t*) = *Ccx*¯(*t*) = *x*2(*t*). It is also noted that the only manipulated input is *uv*. The rest of the elements in *u*¯ are non-manipulated in the three cases.

Now for a discrete-time implementation, the previous linear models are discretized, with a sample time *Ts*, leading to systems of the form

$$\mathbf{x}\_{k+1} = A\_d \mathbf{x}\_k + B\_d \boldsymbol{u}\_k + \boldsymbol{\delta}\_d; \quad \mathbf{y}\_k = \mathbb{C}\_d \mathbf{x}\_k. \tag{20}$$

The MPC calculates control input increments; then each model is augmented to express it in those terms:

$$\zeta\_{k+1} = \underbrace{\begin{bmatrix} A\_d & B\_d \\ 0 & I \end{bmatrix}}\_A \zeta\_k + \underbrace{\begin{bmatrix} B\_d \\ I \end{bmatrix}}\_B \Delta u\_k + \underbrace{\begin{bmatrix} \delta\_d \\ 0 \end{bmatrix}}\_{\delta} \tag{21}$$

$$y\_k = \underbrace{\begin{bmatrix} \mathbf{C}\_d & \mathbf{0} \\ \mathbf{C}\_d & \mathbf{0} \end{bmatrix}}\_{\mathbf{C}} \mathbf{y}\_k \tag{22}$$

where *ζ* = [*x<sup>T</sup> <sup>k</sup> <sup>u</sup><sup>T</sup> <sup>k</sup>*−1] *<sup>T</sup>* and <sup>Δ</sup>*uk* <sup>=</sup> *uk* <sup>−</sup> *uk*−1.

#### *4.3. Prediction Equations*

The predictions equations of the augmented-state, output and input over a prediction horizon *ny*, and control horizon *nu* are:

$$
\underline{\mathcal{G}}\_{k}{}^{k} = P\_{\tilde{\mathsf{V}}\tilde{\mathsf{V}}} \mathsf{\zeta}\_{k} + P\_{\tilde{\mathsf{V}}\Lambda\mathsf{u}} \underline{\Delta\mathsf{u}}\_{k}{}^{k-1} + P\_{\tilde{\mathsf{V}}\delta} \delta\_{\mathsf{v}} \tag{23}
$$

$$
\underbrace{y}\_{\rightarrow}k = P\_{y\zeta}\zeta\_k + P\_{y\Lambda u}\underline{\Lambda}\underline{u}\_{k-1} + P\_{y\delta}\delta\_\prime \tag{24}
$$

$$
\underline{u}\_{k}{}^{k-1} = P\_{\underline{u}\Delta\underline{u}} \underline{\Delta\underline{u}}\_{k-1} + P\_{\underline{u}\underline{\zeta}} \tilde{\zeta}\_{k\prime} \tag{25}
$$

where the arrow notation denotes prediction and is defined as *x* −→*<sup>k</sup>* = [*x<sup>T</sup> <sup>k</sup>*+<sup>1</sup> *<sup>x</sup><sup>T</sup> <sup>k</sup>*+<sup>2</sup> ...] *T*. Augmented state prediction matrices are

$$P\_{\widetilde{\zeta}\widetilde{\zeta}} = \begin{bmatrix} A \\ A^2 \\ \vdots \\ A^{n\_g} \end{bmatrix}, \ P\_{\widetilde{\zeta}\Delta u} = \begin{bmatrix} B & \cdots & 0 \\ AB & \cdots & 0 \\ \vdots & \ddots & \vdots \\ A^{n\_g-1}B & \cdots & A^{n\_g-n\_u}B \end{bmatrix}, \ P\_{\widetilde{\zeta}\delta} = \begin{bmatrix} I \\ I+A \\ \vdots \\ I + \sum\_{i=1}^{n\_g-1} A^i \end{bmatrix}.$$

Output prediction matrices are *Py<sup>ζ</sup>* =diag(*C*)*Pζζ* , *Py*<sup>Δ</sup>*<sup>u</sup>* =diag(*C*)*Pζ*Δ*u*, *Py<sup>δ</sup>* =diag(*C*)*Pζδ*. Input predictions matrices are given by:

$$
\underline{u}\_{k} = \underbrace{\begin{bmatrix} I & 0 & \cdots & 0 \\ I & I & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ I & I & \cdots & I \end{bmatrix}}\_{P\_{\text{adv}}} \underline{u}\_{k-1} + \underbrace{\text{col}(\begin{bmatrix} & 0 & I \end{bmatrix})}\_{P\_{\text{uf}}} \underline{\zeta}\_{k\prime} \tag{26}
$$

and for the states *x* −→*<sup>k</sup>* <sup>=</sup> *Px<sup>ζ</sup> <sup>ζ</sup>* −→*<sup>k</sup>* <sup>=</sup> diag *<sup>I</sup>* <sup>0</sup> ! *<sup>ζ</sup>* −→*k*.

To construct the prediction vectors, full state availability is required. Although *H*-(*t*) is not a measured state and *λ*, *z* are unknown parameters, all are considered known as they could be provided by the EKF estimator. The rest of the states are measured flows and pressures.

#### *4.4. Cost Function and Constraints*

The cost function to be optimized penalizes the future output error with respect to the desired output value *sk* and control input along the prediction and control horizons:

$$J = \sum\_{i=1}^{n\_y} (y\_{k+i} - s\_k)^T \mathcal{Q} (y\_{k+i} - s\_k) + \sum\_{i=0}^{n\_{u-1}} \Delta u\_{k+i}^T \mathcal{R} \Delta u\_{k+i} \tag{27}$$

with Q > 0 and R > 0.

Constraints of the form:

$$\text{diag}(A\_x)\underset{\rightleftarrows}{\underset{\rightleftarrows}{\rightleftarrow}} \leq \text{col}(b\_x); \text{diag}(A\_y)\underset{\rightleftarrows}{\underset{\rightleftarrows}{\rightleftarrow}} \leq \text{col}(b\_y);$$

$$\text{diag}(A\_u)\underset{\rightleftarrows}{\underset{\rightleftarrows}{\rightleftarrow}} \leq \text{col}(b\_u); \text{diag}(A\_{\Delta u})\underset{\rightleftarrows}{\underset{\rightleftarrows}{\left(}} \leq \text{col}(b\_{\Delta u})\tag{28}$$

are considered, where *Ax*, *Ay*, *Au*, *A*Δ*<sup>u</sup>* = [*I*, −*I*] *<sup>T</sup>* and *bx*, *by*, *bu*, *b*Δ*<sup>u</sup>* are vectors that contain the maximum and minimum values allowed in the form [*max min*] *T*.

#### *4.5. Control Law*

Cost function (27) and constraints (28) can be expressed in terms of the decision variable Δ*u* [30], and the following optimization problem is obtained:

$$\begin{aligned} \underline{\Delta\mu}\_{k-1}^{\*} &= \underset{\underline{\Delta\mu}\_{k-1}}{\text{arg min}} \left\{ \frac{1}{2} \underline{\Delta\mu}\_{k-1}^{T} \mathcal{H} \underline{\Delta\mu}\_{k-1} + \mathcal{F}^{T} \underline{\Delta\mu}\_{k-1} \right\}, \\ &\text{s.t. } \mathcal{M}\_{\mathfrak{c}} \underline{\zeta}\_{k} + \mathcal{N}\_{\mathfrak{c}} \underline{\Delta\mu}\_{k-1} \le f\_{\mathfrak{c}}. \end{aligned}$$

with

$$\begin{aligned} \mathcal{H} &= P\_{y\Delta u}^{\mathsf{T}} \operatorname{diag}(\mathcal{Q}) P\_{y\Delta u} + \operatorname{diag}(\mathcal{R}); \\ \mathcal{F} &= P\_{y\Delta u} \operatorname{diag}(\mathcal{Q}) (P\_{yx}\mathcal{I}\_{k} + P\_{y\delta}\delta - \underline{\hat{s}}\_{\delta}); \\ M\_{\mathsf{c}} &= \begin{bmatrix} \operatorname{diag}(A\_{x}) P\_{x\zeta} P\_{\xi\zeta} \\ \operatorname{diag}(A\_{y}) P\_{y\zeta} \\ \operatorname{diag}(A\_{u}) P\_{u\zeta} \\ 0 \end{bmatrix}; \\ N\_{\mathsf{c}} &= \begin{bmatrix} \operatorname{diag}(A\_{x}) P\_{x\zeta} P\_{\xi\Delta u} \\ \operatorname{diag}(A\_{y}) P\_{y\Delta u} \\ \operatorname{diag}(A\_{u}) P\_{u\Delta u} \\ \operatorname{diag}(A\_{\Delta u}) \end{bmatrix}; \\ f\_{\mathsf{c}} &= \begin{bmatrix} \operatorname{col}(b\_{x}) \\ \operatorname{col}(b\_{y}) \\ \operatorname{col}(b\_{u}) \\ \operatorname{col}(b\_{\Delta u}) \end{bmatrix} - \begin{bmatrix} \operatorname{diag}(A\_{x}) P\_{x\zeta} P\_{\zeta\delta} \\ \operatorname{diag}(A\_{y}) P\_{y\delta} \\ 0 \\ 0 \end{bmatrix} \delta. \end{aligned}$$

Finally, the control law is given by

$$
\mu\_k = \mu\_{k-1} + \Delta u\_{k'}^\* \tag{29}
$$

where Δ*u*∗ *<sup>k</sup>* is the first element of <sup>Δ</sup>*u*−→<sup>∗</sup> *<sup>k</sup>*−1.

#### *4.6. Offset-Free Tracking*

The MPC presented earlier does not contain an explicit mechanism to deal with disturbances and modeling errors that arise from set point changes that move the system away from the original operation point where it was linearized. Therefore, an integral action must be embedded into the control law to deal with these situations and guarantee zero steady-state error. This is performed including in the state space model an integrating state

$$
\begin{split}
\begin{bmatrix}
\mathfrak{F}(t) \\
\mathfrak{F}(t)
\end{bmatrix} &= \begin{bmatrix}
A\_c & 0 \\
\end{bmatrix} \begin{bmatrix}
\mathfrak{F}(t) \\
\mathfrak{F}(t)
\end{bmatrix} + \begin{bmatrix}
B\_c \\
0
\end{bmatrix} \vec{u}(t) + \begin{bmatrix}
\delta\_c \\
s(t)
\end{bmatrix}, \\
y(t) &= \begin{bmatrix}
\end{bmatrix} \mathsf{F}\_c \quad \text{0} \quad \Big|
\begin{bmatrix}
\mathfrak{F}(t) \\
\mathfrak{F}(t)
\end{bmatrix}.
\end{split}
$$

A stabilizing feedback gain can be computed to guarantee that in steady-state *y*(*t*) = *x*2(*t*) = *s*(*t*). Another approach to obtain integral action could be to estimate the disturbance with the use of a disturbance observer [31].

#### *4.7. Set Point Manager*

The set point manager (SPM) is a ruled-based algorithm that takes into account the periods of expected maximum and minimum demand through time conditions, that is, it identifies the hours (early morning) in which it is possible to reduce the pressure to reduce the leak magnitude without affecting the demand. Therefore, when the system is leak-free (i.e., *λ* = 0), the SPM computes the pressure reference (*sk*) by solving the hydraulic model to ensure that the demand is satisfied. Then, if there is no fluid loss, the priority is to supply to the WDS the pressure needed to fulfil the users' demand in both maximum and minimum demand hours.

On the other hand, the SPM receives the leak magnitude estimate through the leak coefficient *λ<sup>k</sup>* when a leak occurs. The SPM identifies if it is possible to reduce the pressure reference according to the day's time, prioritizing meeting the DP. If the DP is maximum, the pressure reference is set to the control pressure (*H*ctrl) calculated by simulation. Otherwise, the pressure reference is reduced to meet the minimum DP to minimize fluid loss without neglecting the required demand. This procedure is formally described in Algorithm 1.

#### **Algorithm 1:** Set point manager for pressure management

**Input**: One-day demand profile (DP), time-varying leak coefficient (*λk*), leak detection threshold (*λ*th). **Output**: Time-varying set point (*sk*). Initialization: **for each** index *i* **in** DP **do**

```
Compute the control pressure (Hctrl) to ensure
scheduled demand in DP assuming leak-free
conditions:
```

```
Hctrl[i] ← solveHydraulics(DP[i])
```
**end for**

Online operation: **for each** time-step *k* **do**

**if** *λ<sup>k</sup>* > *λ*th

```
Sk ← min(Hctrl)
  else
    tk ← getCurrentTime( )
    Sk ← interp(Hctrl, tk)
  end if
end for
```
The solveHydraulics() subroutine in Algorithm 1 solves the hydraulic model to compute the required pressure in the control node, ensuring that the demand is satisfied at each hour of the day. The interp() subroutine computes by interpolation the control pressure corresponding to each time step *tk*.

#### **5. Simulation Results**

This section presents simulations of the WDS built in the Hydroinformatics Laboratory at the Technological Institute of Tuxtla Gutierrez, whose mathematical model was validated in [23]. The system parameters are given in Table 1. The system can be configured as a single horizontal pipeline and a branched network as in Figures 1, 3, and 4. An EKF-based method [24] has been considered to estimate the leak position *z*-, its magnitude *λ*ˆ , and the leak pressure *H* for Cases 1 and 2. A leak occurring in an accessory (pipeline joint) located in an unknown position is assumed for Case 3; nevertheless, it can be estimated by an EKF, e.g., [25]. However, it is important to note that the proposed method is not attached to any leak estimation method, and it could be generalized to any WDS within the topology presented here.

The initial conditions for the flows *Qi*(0) and pressures *Hi*(0) for all cases are presented in Table 2. The sampling period is *Ts* = 0.01 s. Gaussian noise with a variance of 2.53 × <sup>10</sup>−<sup>10</sup> <sup>m</sup>6/s2 for the flow rate and 3.72 × <sup>10</sup>−<sup>4</sup> <sup>m</sup><sup>2</sup> for the pressure were added to the signals. These noise levels were characterized according to the response of the Yokogawa sensors installed in the physical system as described in [32]. Cases 1 and 2 take into account an uncertainty for *z*-, *λ*, and *H* of ±15%, that is, an error in the estimation of the leakage variables. For Case 3, the same value for *λ* as the previous cases is considered. For the MPC, the prediction horizon is *ny* = 15, the control horizon is *nu* = 3, and the weights are Q = 1 and R = 0.5. These values are based on well-known tuning methods [33]. The constraints were proposed according to the physical behavior and expected DP. For all cases, 0.1 ≤ *r* ≤ 1. For Case 1, the constraints on the states are

$$\begin{aligned} 1 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s} &\leq Q\_{1,2,3} \leq 1.3 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s},\\ 2 \,\mathrm{m} &\leq H\_3 \leq 2.87 \,\mathrm{m}. \end{aligned}$$

**Table 1.** System parameters.



**Table 2.** Initial conditions for the three simulation cases.

For Case 2,

$$\begin{aligned} 1.5 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s} &\leq Q\_{1,2,3} \leq 2 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s}, \\ 4.3 \,\mathrm{m} &\leq H\_2 \leq 5.4 \,\mathrm{m}. \end{aligned}$$

Finally, for Case 3,

$$\begin{aligned} 0.17 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s} &\leq Q\_1 \leq 3.2 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s}, \\ 0.7 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s} &\leq Q\_2 \leq 1.6 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s}, \\ 0.4 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s} &\leq Q\_3 \leq 0.7 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s}, \\ 0.95 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s} &\leq Q\_4 \leq 1.5 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s}, \\ 0.4 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s} &\leq Q\_5 \leq 0.7 \times 10^{-3} \,\mathrm{m}^3/\mathrm{s}, \\ 1.2 \,\mathrm{m} &\leq H\_2 \leq 2 \,\mathrm{m}. \end{aligned}$$

These constraints were chosen to satisfy a minimum demand, even in the case of a leak.

The simulation covers a period of 24 h whose maximum and minimal magnitudes are concerning a typical demand profile. For Case 1, the results are displayed in Figure 7. The top plot of Figure 7 shows the pressure in the demand node. The dotted line is the DP, the dashed is the reference driven by the SPM that adapts *Sk* during the leak period to reduce the water losses, and the solid line represents the controlled pressure (*H*3) in the node. The leak occurs at *t* = 3 h with the conditions given in Table 1. The SPM identifies a time of minimum demand and reduces the pressure. This pressure reduction remains until *t* = 7 h where the maximum demand period starts. During the period of maximum demand, the SPM does not adjust the reference because the priority is to satisfy the demand. After *t* = 17 h, the SPM again reduces *Sk* to reduce the leak magnitude. The middle plot of Figure 7 shows the effect of leak reduction due to the SPM. The solid line is the leak magnitude without the SPM and the dashed-line with the SPM. The bottom plot of Figure 7 shows the PRV opening to track the desired reference. It is important to note that all the constraints are satisfied for the MPC. We can conclude that the MPC tracks the set point with good performance even during the leak period, which demonstrates the effectiveness of the proposed method.

**Figure 7.** Case 1: Pressure control and management in a pipeline with a leak before the controlled node.

For Case 2, the results are displayed in Figure 8. As can be seen, the MPC tracks the pressure reference provided by the SPM even in a leakage scenario. Similar to Case 1, the SPM automatically adjusts the controlled node pressure at the maximum and minimum demand periods to reduce water losses, and all the constraints are satisfied for the MPC. To test the robustness of the controller, ±15% of uncertainty in the location of the leak is added, demonstrating the performance of the controller even when there is a leak location error.

To contrast the results with the traditional method of control and management of pressure in a pipe, Figure 9 presents a comparison between a traditional PID controller and the proposed predictive control scheme. The PID gains were computed with the Matlab- PIDtuner-, obtaining the optimal values for proportional *KP* = 1, integral *KI* = 0.5, and derivative *KD* = 0.5 gains. It can be seen that the PID tracks the demand profile with good performance. However, in the event of a leak, the set point remains the same, and the leak is seen by the controller as a disturbance. So the PID adjusts (increases) the pressure to track the original demand profile. This way of operating guarantees the pressure in the demand node but at the cost of fluid loss. In contrast, the MPC is aware of the leak event and with the aid of the set point manager adjusts the pressure operating point to a lower value reducing the fluid loss and delivering the minimum allowed pressure to the demand node.

**Figure 8.** Case 2: Pressure control and management in a pipeline with a leak after the controlled node. Including ±15% of uncertainty in the leak location.

**Figure 9.** Case 2: Comparison between a traditional PID controller and the proposed control scheme.

Results for Case 3 are shown in Figure 10, with a leak being simulated at a connection. The MPC still achieves a good tracking of *Sk* that the SPM sets, which proposes a pressure reduction at the time of the leak. This pressure reduction remains constant in the first period until the maximum demand. The middle plot of Figure 10 shows the leak magnitude with and without the SPM. Note that Case 3 considers three demand nodes (due to its branches) controlled by adjusting the pressure on the controlled node. Therefore, it is essential to analyze these flows to evaluate the MPC effectiveness. The magnitudes with and without the SPM are represented by a dashed and a solid line, respectively. As can be seen, the priority is to satisfy each node demand during the maximum demand periods. However, during the minimum demand periods, the pressure is reduced, which reduces water losses.

Finally, Figure 11 displays the flow rates at the demand nodes for Case 3: *Q*3, *Q*4, and *Q*5. The solid lines are the flows without the effect of the SPM, and the dotted lines

represent the reduction due to the SPM. The reduction of flows is within the minimum permissible limits given by constraints on the states taken into account by the MPC and only during the hours of minimum demand.

**Figure 10.** Case 3: Pressure control and management in a branched water distribution network.

**Figure 11.** Case 3: Flows in the main water intakes for the WDS.

For all the cases, the trade-off between maintaining a minimum demand and reducing the leak magnitude was achieved due to the good performance of the MPC and the SPM algorithm. As a result, a reduction of water losses of ≈5% was accomplished, which is important due to the fact that water distribution systems must operate without interruptions throughout the year.

#### **6. Conclusions**

This work proposed a model-based predictive controller for managing and controlling water distribution systems. The proposed method seeks a trade-off between maintaining water supply and reducing water losses due to leaks. To achieve this goal, the MPC is based on the water-hammer equations of the hydraulic system together with physical constraints and an adaptive demand profile managed by a set point manager algorithm. The control input was calculated by minimizing the output errors with respect to the demand profile driven by the set point manager. The control objective was formulated with state, input, and output constraints on the cost function. Moreover, the control scheme includes a strategy to handle the strong nonlinear behavior of the PRV. The controller was tested with numerical simulations in a model characterized by a real pipeline and pipe network located at the Hydroinformatics Laboratory of the Technological Institute of Tuxtla Gutiérrez. Therefore, this work presents a mathematical model based on a real system and realistic operating conditions. The simulation results illustrate the performance and robustness of the MPC for pressure management in the system and the reduction of leaks due to water losses, with an average of 5% in the presence of noise, disturbances, and uncertainty in the leak location estimate. Future work will extend this work to an integrated methodology of a multi-leak tolerant control algorithm.

**Author Contributions:** Conceptualization, J.-R.B., F.-R.L.-E. and G.B.; methodology, J.-R.B., F.-R.L.-E. and G.B.; software, J.-R.B. and G.V.-P.; validation, J.-R.B., G.V.-P. and I.S.-R.; formal analysis, J.-R.B., F.-R.L.-E. and G.B.; data curation, G.V.-P. and I.S.-R.; writing—original draft preparation, J.-R.B., F.-R.L.-E. and G.B.; writing—review and editing, G.V.-P. and I.S.-R.; visualization, I.S.-R.; supervision, F.-R.L.-E. and G.B.; project administration, F.-R.L.-E. and G.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has been supported by the Consejo Nacional de Ciencia y Tecnología (CONA-CyT) and by Tecnológico Nacional de México under the program *Proyectos de Investigación Científica y Desarrollo Tecnológico e Innovación 2022*.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Marko Miloševi´c 1,\*, Milan Radi´c 1, Milica Raši´c-Amon 1, Dragan Litriˇcin <sup>2</sup> and Zoran Staji´c <sup>1</sup>**


**Abstract:** This paper indicates the importance and advantages of the application of hybrid models in the control of water supply systems. A range of possibilities provided by this scientific approach is presented in the practical examples related to the fault diagnostics and fault-tolerant control in the pumping station (PS) control systems. It is presented that continuous monitoring and recording of the data of the pumping stations operation processes (electrical parameters such as electrical power, pressure or flow in the pipelines, water levels in the tanks, changes in various discrete states, etc.) could represent a significant resource that can be used to develop various hybrid models using the appropriate "data-driven" techniques. During this process, data are transformed into information, and thereafter, information into knowledge. Based on this knowledge, the control of PS operation can be significantly improved and a significant increase in the user's satisfaction can be achieved while the maintenance and operation costs can be reduced.

**Keywords:** water supply systems control; pumping stations operation; hybrid models; "data-driven" models; fault diagnostics; fault-tolerant control

#### **1. Introduction**

A hybrid model is a combination of two or more different modeling techniques, which aims to improve the modeling adequacy of complex systems [1–9]. In the processes of system design or optimization, the principles of mathematical modeling are usually used. These modeling principles include the representation of physical models using appropriate systems of equations (differential, algebraic) that describe the relations between the system's variables and the system's parameters. These models are known as parametric models [10–17]. Due to known relationships between the model parameters and variables, they are usually presented as "white boxes" ("white box models"). Despite the complexity of these models, they are not able to present the ideal behavior of the real system in all operation conditions during its exploitation. This is especially difficult for systems that are characterized by certain unexpected occurrences, irregular conditions, or failures that are difficult to predict in the systems' (re)designing phase. In situations where certain data recorded during such systems' operation are available, it is possible to explore and discover various phenomena related to the systems' processes that have not been previously noticed or sufficiently researched. It is also possible to recognize some additional relationships between the process parameters and process variables. There are also examples where it is not possible to mathematically describe these relationships. These hidden relations in data are mostly extracted through some machine learning systems [9,18–22], which are used in gaining a better understanding of the process or to help in the decision making [9,15,23,24]. These models are known as non-parametric models, presented as "black boxes" ("black box models"), and they are classified in the "data-driven" models [2,8,9,25–28]. The combination and integration of classical parametric models with one or more non-parametric models, connected into the serial, parallel, or serial–parallel structures is known as the

**Citation:** Miloševi´c, M.; Radi´c, M.; Raši´c-Amon, M.; Litriˇcin, D.; Staji´c, Z. Diagnostics and Control of Pumping Stations in Water Supply Systems: Hybrid Model for Fault Operating Modes. *Processes* **2022**, *10*, 1475. https://doi.org/10.3390/pr10081475

Academic Editors: Francisco Ronay López-Estrada and Guillermo Valencia-Palomo

Received: 18 April 2022 Accepted: 24 May 2022 Published: 27 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

hybrid model in the literature. An example of the serial–parallel hybrid model structure is presented in Figure 1.

**Figure 1.** Serial-parallel hybrid model structure.

This modeling method allows that all of the available knowledge can be integrated into a single model without compromising the adequacy of the basic parametric model [2,8,9]. On the other hand, the quality of the model increases and processes from the reality are better modeled [6,9,13]. In addition, the information and knowledge that can be obtained from the non-parametric models or simulation results from hybrid models can be used to improve the quality of the system control. Thereafter, these models represent significantly better solutions for the modeling of the complex processes and systems in comparison to the models that use only one source of knowledge (only one modeling method) [1–9,29].

Some authors in this area recognize that a "data-driven" approach can also be applied in situations where there is a lack of information about the functioning of the real system, where hybrid models can be expanded by adding data from different sources, and that hybrid models are especially suitable for a range of applied targeted research [2,9,30]. Other authors have concluded that hybrid modeling techniques are not unknown and that they are mainly applied in scientific (academic) research, but they are insufficiently applied in practice [9].

Therefore, this paper presents the applied research in the area of pumping station control, and the main goal was to demonstrate the advantage of using hybrid models in the field of diagnostics and the tolerant control of the irregular operation modes as well as the fault-tolerant control in the water supply systems. The authors chose the concept with a number of simpler non-parametric models, because it allows for "step-bystep" expansion and an increase in the complexity of the hybrid model and facilitates its application in practice.

#### **2. Hybrid Models of the Pump Stations Automatic Control Systems**

In the process of creating a hybrid model of a certain PS, it is necessary to start from its physical model and form an appropriate parametric model based on it.

As an example of this procedure, in Figure 2, the physical model of PS with two pumping units (PU) is presented. One of the units is in operation and the other is spare, which is a common situation for a smaller PS in water supply systems. For easier understanding, in Figure 2, the water flows are presented in blue, the power flows that feed the drive motors are shown in red, the control commands are in orange, and the reference value according to which the control is performed (control signal) is in green.

**Figure 2.** The physical model of a PS with two PUs pumping water into the discharge–distribution piping system with one counter-tank.

It is assumed that it was a booster PS (i.e., that PUs pump water from the common suction line that branches to the suctions of both pumps). Pump discharge lines are also connected within the PS into the single discharge–distribution pipeline that represents the situation regularly presented in practice. This pipeline supplies the connected customers and excess water is usually transferred to a counter-tank located at the appropriate location above the PS to ensure the stability of pressure in the distribution system and reduce its variations during PU switching on and off and/or major changes in the water consumption of the system.

Depending on the water consumption and water flows in the system, it is possible at certain periods that the total water flow to consumers (∑<sup>n</sup> <sup>i</sup>=<sup>1</sup> Qconsi ) is higher than the flow provided by the PU in operation (Qout), which leads to the decrease in water levels in the tank. This is presented in Equation (1):

$$\mathbf{Q\_{res}} = \mathbf{Q\_{out}} - \sum\_{i=1}^{n} \mathbf{Q\_{cons\_i}} < \text{ } \text{ } \tag{1}$$

However, such systems are usually designed so that in regimes when one PU operates, the total water flow to consumers is less than the flow provided by PU in operation for much longer periods, so that excess water fills the tank and increases its water levels (Equation (2)):

$$\mathbf{Q\_{res}} = \mathbf{Q\_{out}} - \sum\_{i=1}^{n} \mathbf{Q\_{cons\_i}} > \begin{array}{c} \text{0.} \tag{2}$$

All examples presented in this paper refer to the PS Knez Selo, whose operation has been monitored since 8 June 2009. During that period, there were several PU replacements and reconstructions of parts of the electrical installations belonging to the PS:


For the analyses presented in this paper, it is important to point out that in all of the presented examples, there was no regulation of water flow in the system. All flow changes were exclusively a consequence of the process of self-regulation due to natural changes in the system pressure in the suction and discharge pipelines of the PU (changes in the PU head). Changes in the discharge pipeline pressures come as a consequence of changes in the water consumption, or changes in the counter-tank water level. On the other hand, pressure in the suction pipeline changes due to the switching of the PUs in the PS supplying water to consumers in lower altitude zones (PS Matejevac 2), and variations in the water consumption in that part of the system. Even in the presented examples in which the PUs were started via frequency converters, their output frequencies were constant and there was no regulation of the speed of the pump impellers. Any small changes in the speed of rotation of the pump impellers were due to changes in the load on the motor shaft, which was also a consequence of the natural change in the pump head. In all illustrative cases, the valves located on the suction and discharge branches of each PU were fully open and no flow control was performed with them. For this reason, their positions are not shown in the physical model presented in Figure 2.

#### *2.1. Parametric Models of PS*

There are a large number of scientific papers that have described the different parametric models of PS, which can be used to simulate their operation [10–17,26,31]. Bearing in mind the topic of this paper, the focus of research was not placed on the well-known parametric models of the PS. For the investigation covered in this paper, the parametric model of PS presented by the authors in [12] can be used as the basis.

#### 2.1.1. A Simple Parametric Model of PU Convenient for the Analysis of PS Operation

This model indicates that relations between electrical (phase and line voltages and currents, active, reactive, and apparent powers, power factors, frequency) and non-electrical (flow, pressure, and pumP s efficiency) parameters can be established. This can be achieved by providing appropriate H(Q) and P (Q) characteristics [12,15–17,32]. These characteristics can be represented by the corresponding polynomials of the fourth- and third-order, respectively:

$$\text{H(Q)} = \text{H}\_0 + \text{H}\_1 \ast \text{Q} + \text{H}\_2 \ast \text{Q}^2 + \text{H}\_3 \ast \text{Q}^3 + \text{H}\_4 \ast \text{Q}^4;\tag{3}$$

$$\mathbf{P}'(\mathbf{Q}) = \mathbf{P}\_0 + \mathbf{P}\_1 \ast \mathbf{Q} + \mathbf{P}\_2 \ast \mathbf{Q}^2 + \mathbf{P}\_3 \ast \mathbf{Q}^3. \tag{4}$$

In certain situations, it is assumed that the characteristic P (Q) is linear, meaning that P2 = 0 and P3 = 0, which will be adopted in the following examples.

The PU efficiency can be calculated based on the following formula:

$$
\eta(\mathbf{Q})\,\left[\%\right] = \frac{\mathbf{P}(\mathbf{Q})}{\mathbf{P}'(\mathbf{Q})} \ast 100 = \frac{0.981 \ast \mathbf{Q} \ast \mathbf{H}(\mathbf{Q})}{\mathbf{P}'(\mathbf{Q})}.\tag{5}
$$

In Equations (3)–(5), P is given in [kW], H in [m], and Q in [l/s].

#### 2.1.2. The Parametric Model Application—Obtaining the Simulation Results

Knowing the parametric models of PU, historical data of the measured electrical active power P can be used to indirectly determine the flow and pressure of the pump, and the PU efficiency in the same historic period. Furthermore, the results of modeling obtained in this way can be used for the estimation of the total volume of water pumped into the system in a certain period, or for energy efficiency analyses. Additionally, the basic idea of using the historical data could be to find out what really happened during the PS operation in the past. Fixed data sequences used as data input in the parametric model enable the detection of different types of irregularities, faults, or failures that have appeared during the PS operation as well as to explore the dependencies between the different system variables and corresponding system parameters. This is an effective way in which to transform the data into valuable information that can help in diagnosing different technical problems and in finding appropriate solutions for them.

Values of all of the electrical quantities measured by power loggers PL1 and PL2 (Figure 2) were recorded with a sampling period of 12 s. Some of them are shown in Figures 3–6 as the daily diagrams of changes in the recorded electrical quantities on 30 December 2009, but the principle of using the parametric model is identical for any other period.

**Figure 3.** The diagrams of the phase and line voltage changes U(t) at PS Knez Selo on 30 December 2009.

**Figure 4.** The diagrams of the PU2 load current changes I(t) on 30 December 2009.

**Figure 5.** The diagrams of the PU2 active, reactive, and apparent power changes on 30 December 2009.

**Figure 6.** The diagrams of the PU2 power factor changes on 30 December 2009.

It should be pointed out that the only electrical quantity needed for the presented parametric model application is active power P . In order to obtain the corresponding simulation results, only this diagram is actually needed (turquoise line in Figure 5). All of the other diagrams presented in Figures 3–6 can be used to explore the dependencies between the different system variables and the corresponding system parameters.

The following values of the coefficients in Equations (3) and (4) were used when calculating the non-electric quantities: H0 = 144, H1 = −4.7273, H2 = 0.3413, H3 = −0.07267, H4 = 0.00102, P0 = 9.95, P1 = 0.620448. The same coefficients were used for both PUs, because the experimental tests performed in the period presented in the figures showed that they had approximately the same operating characteristics.

The calculation procedure was as follows: Based on the information on the electrical active power P that the pump drive motor draws from the supplying electrical system at a certain time moment, the corresponding flow value Q can be determined using Equation (4). For the presented linear function, it is not complicated at all. By substituting this value of the flow in Equation (3), the value of the head H is determined, and by substituting the same value in Equation (5), the corresponding value of the PU efficiency η at that time moment is obtained. This calculation procedure is repeated, so that for each subsequent recorded value of the active power P , the values of the corresponding quantities Q, H, and η are obtained. Presenting these obtained results in the time moments in which there are available data on the active power P (every 12 s), the diagrams of the Q, H, and η changes in time are obtained (Figures 7–9).

**Figure 7.** The PU2 flow change diagram Q(t) on 30 December 2009.

**Figure 8.** The PU2 head change diagram H(t) on 30 December 2009.

**Figure 9.** The PU2 efficiency change diagram η(t) on 30 December 2009.

Bearing in mind that the values of all of these quantities are known at each of the observed moments, their interdependencies can also be shown. Figures 10–12 show the functional dependencies H(Q), P (Q), and η(Q). Unlike the diagrams presented in Figures 3–9, which show the daily changes in the different process variables on 30 December 2009, the images in Figures 10–12 refer to the entire month (December 2009).

**Figure 10.** The operating point positions on the PU1 and PU2 H(Q) characteristics in December 2009.

**Figure 11.** The operating point positions on the PU1 and PU2 P (Q) characteristics in December 2009.

**Figure 12.** The operating point positions on the PU1 and PU2 η(Q) characteristics in December 2009.

2.1.3. Practical Considerations and Simulation Analysis

When observing the time diagrams or functional dependencies in the previously presented figures, the following considerations must be taken into account:


in the suction branches of the PU at PS Knez Selo, which resulted in an increased pump head (Figure 8) and decreased water flow (Figure 7). Changes in the active power diagram shown in Figure 5 are also direct consequences of these phenomena. Precisely, such recognized relations make it possible to identify certain phenomena occurring in the larger part of the hydromechanical system, based on the recorded electrical quantities changes.


with another topic, this model will not be further developed, and the original recorded data will be used in further analyses.


#### 2.1.4. Control Algorithms and Automatic Control System Functioning

In order to achieve a better understanding of the phenomena being observed in some PSs and to be able to form a parametric model of PS operation, it is necessary to understand the logic (i.e., the control algorithm applied during the implementation of the automatic control system (ACS) in the PS, the reference value used for the control, and the applied equipment).

In pump stations of the same type as the pump station presented in Figure 2, the "on–off" model of control is usually provided. The control is based on the information on the actual water level in the tank (h). The main goal of the regulation is the requirement that the water level (h) always remains between the previously defined minimum (hmin) and maximum (hmax).

The Knez Selo water tank has a total volume of 100 m3, and the changes in the water level was monitored using five conductive probes placed in the tank at different heights. The system monitors and records changes that occur when any of the probes becomes submerged (which is signaled with the blue dot) or when it remains dry (which is signaled with the red x sign). Figure 13 shows the reactions of the upper- and the lower-limit probe in the Knez Selo water tank on 30 December 2009. Based on the time schedule of these signals, the change in the water level in the tank can be modeled. Comparing Figure 13 with the diagrams of changes in the quantities shown in Figures 3–9, it is easy to note that the limit probes controlled the operation of PU2 in the PS. PU1 was not put into operation on 30 December 2009.

**Figure 13.** The diagram of the water level change h(t) in the Knez Selo water tank on 30 December 2009.

Observing Figure 13, one can see that only three of the existing five probes detected a change in the state during 30 December 2009. Two of them were used as limit probes (the upper one was set at hmax = 3.3 m and the lower one at hmin = 1.1 m height from the bottom of the tank). The third one was placed in the middle between them (i.e., at a height of 2.2 m). There were two more probes in the tank. One of them signaled the appearance of an overflow in the tank and was placed at a 3.5 m height, between the upper limit probe and the overflow pipe, which was physically located in the tank at a 3.6 m height (maximal water level in the tank). The other one was placed at a height of 0.2 m from the bottom of the tank, and when it remains dry, it means that the water tank is almost empty. During 30 December 2009, there was no change in the condition of these two probes.

With this arrangement of probes, it was chosen to turn off a PU when the volume of water in the tank reached a value of 91.6 m3, and to turn it on when the volume of water in the tank dropped to about 30 m3. Setting the lower limit probe to a lower level is not recommended because a certain volume of water in the tank must serve as a reserve to supply consumers from the moment of any problem in the PS until the arrival of the operator tasked with troubleshooting. A water volume of 30 m<sup>3</sup> allows for a water supply autonomy of at least 1 h. In most situations, it should be sufficient for the operator to be notified and to attend to the PS in order to fix the problem.

Considering the position of the upper limit probe, it might make sense to raise it by 0.1 m, which could increase the total available volume of water in the tank to about 94.4 m3, but this could result in more frequent overflows in the tank.

In some cases, such as in the examples shown in the previous figures, the choice of PU to be in operation was selected manually. By switching the appropriate switch, the operator determines which PU will be in operation during the next period, and until the next switching, the local automation turns on/off the selected PU only.

#### 2.1.5. The PU Parametric Model Application in More Complex Cases

After the last PS Knez Selo reengineering, performed in October 2017, two power loggers that were previously used to measure and record the electrical parameters of each PU were dismounted. Instead, a single power logger was installed in order to measure and record the parameters related to the electrical energy consumption of the complete PS. On that occasion, the pump units and the complete ACS were replaced. Two completely new PUs run by the drive motors with a rated power of 9 kW were installed. During this period, the PUs were started via appropriate frequency converters, but there was no frequency control (their output frequency was set to 50 Hz).

Relating to the parametric model application, it is a significantly more complex case. Namely, the main advantage of the previously applied concept represents the possibility of direct measurements on each PU (direct model application is possible), but there was not any information on the consumption of other electronic devices or appliances installed in the PS. The second concept, applied after October 2017, provides information about the total PS's electrical energy consumption and all the relevant parameters (Figure 14). However, the parameters related to the consumption of one or another PU cannot be obtained directly. This information has to be extracted from the originally recorded data using some kind of machine learning approach. Due to the single-phase consumers installed in the PS, it is hard to perform without the simultaneous monitoring of the parameters in all three phases.

During the implementation of the new ACS, the requirement of the cyclic operation was set—at every following achievement of the hmin, the command for the start of the PU that in the previous cycle was off is given. A part of a simplified algorithm that meets the control requirements is presented in Figure 15. The practical result of its application in the normal operation of the PS is presented in Figure 14. Actually, in Figure 14, the daily diagrams of the load current changes at PS "Knez Selo" are presented. These diagrams were recorded on 3 January 2022 with the sampling period of 10 s. The order of PUs starting is indicated by the appropriate ordinal number in Figure 14.

According to the goals of the research, in this paper, the focus will be on the "datadriven" non-parametric models, algorithms that can be used to implement certain hybrid models, and/or knowledge that can be used to improve the ACS or its control algorithms [7–9,15,33–35].

**Figure 14.** The diagrams of the PS load current changes recorded on 3 January 2022. 1: the PU1 is in operation, 2: the PU2 is in operation.

**Figure 15.** A part of the control algorithm allowing cyclic PU operation.

#### *2.2. What Is the Purpose of Non-Parametric Models in the Simulations of PS Operation?*

All predictable and expected situations and phenomena that occur in the PS operation can be modeled by some of many parametric models available in the literature. In all these cases, based on the simulation results, it is possible to (re)design the complete PS including appropriate automatic control and protection systems.

However, there are some specific situations and phenomena that appear in the PS operation that are not easily predictable in the system-design phase such as many of the ones belonging to the field of fault diagnostics and fault-tolerant control. These phenomena could have many negative consequences on the PS operation such as installed equipment damage, reducing the operating life, increase in the operating and maintenance costs, etc. In the hybrid modeling theory, the combined use of physics-based and data-driven models is of great importance. Some authors call this approach hybrid analysis and modeling (HAM), and consider it the third monumental step in the development of science [24].

As one of the examples showing that this is indeed an approach with huge potential, it can be stated that the considerations given in Section 2.1.3 have shown that the leaps in active power observed in the corresponding diagram in Figure 5 (in periods when some of the PU works) can be converted into useful information about the on/off of the PUs installed at PS Matejevac 2, supplying the lower altitude zone as well as PS Knez Selo. For example, this can be used to signal possible disturbances in the operation of PS Matejevac 2 (e.g., longer periods in which the PUs in this PS are not switched off, which could be due to overflows in some tanks, burst pipelines, or leaching in parts of the water supply network located in lower altitude zones). Generally speaking, the possibility of obtaining information on the operation of other PSs in the system or other parts of the system, based on data on the operation of subjected PUs, provides enormous opportunities for various innovations and improvements in the quality of system operation in regular or irregular conditions, malfunctions, or various emergency situations.

In this sense, for all of the data-driven models, the data recorded during the real system exploitation are incomparably more important for this area than any of other data (e.g., data from different simulation models). For this purpose, data recorded during the real system's operation were selected as the subject of the research in this paper. Based on these data, the appropriate "data-driven" parametric model was created in Section 2.1.

The task of the parametric models that were the subject of this paper was to process the available historical data on the system's operation and to identify the different unexpected phenomena, irregular modes, or faults that occurred in the PS. Many of them could not be described by some of the new parametric models, and would be treated as "black boxes". However, the aim should be to provide enough information about them for the appropriate analyses. These analyses can be used to obtain other useful information and knowledge based on the changes in the values or frequencies of the changes in certain process data and their interrelationships. The obtained information and knowledge could be used to create appropriate non-parametric models that will, in combination with the existing parametric model, enable better modeling and simulation of that phenomena. On the other hand, it can be expected that the obtained information and knowledge can be used to improve the existing automatic systems' control and protection or applied control algorithms. Practical reasons for these improvements are usually to raise the reliability and quality of the customer service as well as the level of the protection and safety of the facilities and equipment, thus increasing the energy efficiency of the system and reducing the maintenance and operating costs, etc.

It was previously stated that the hybrid model consists of a parametric and a number of non-parametric models whose integration should be performed in a certain way. When hybrid models are applied for the purpose of modeling certain processes, then the result of such models should be appropriate simulation results or appropriate predictions. In order to facilitate the model integration, the authors also deliberately chose to use the experimentally recorded data of the electrical quantities related to the loads of the pump drive motors as input data for the parametric model.

#### **3. Example of the Creation of Non-Parametric and Corresponding Hybrid Models for Specific Phenomena Detected in the PS Operation**

#### *3.1. The Outdoor Lightning Detection Based on PS Load Currents*

The diagrams of the PS load current changes presented in Figure 14 can be used as a simple example that can practically verify the previously described approach to nonparametric model creation. Namely, with a more precise look at the load current diagrams, it can be noticed that at the beginning and end of the day, the load current in phase L1 (red line) was significantly higher than the other two currents, while in the middle part of the

day (approximately between 7:30 and 16:40), it was very close to the load currents in phases L2 and L3 (its value was between two other currents). Knowing that this phenomenon is repeated every day, it can be concluded that the cause of its appearance was the operation of the outdoor lighting of the PS, which was turned off in the morning and turned on in the evening. This diagram shows that at the time when the external lighting was off (the electrical circuit powering the external lighting was interrupted), the currents in all phases had approximately the same values. During these periods, the electricity in the PS is mostly consumed by the PUs' three-phase induction motors (strictly speaking, part of the energy is also spent on the no-load consumption of frequency converters and elements of control and protective devices, but it can be neglected). In periods when the outdoor lighting system is on, the energy required for its operation is taken over through phase L1, to which this system is connected. The characteristic time moment of this system switching on, when the peak in the load current appears, can be clearly seen in this diagram. It is a direct consequence of the characteristic transient processes taking place in light sources during their heating from the ambient to operating temperature.

Furthermore, based on the load current change diagrams, the moments of switching off and on the external lighting could be accurately detected every day. Knowing the expected periodicity during the year and voltage changes in the power distribution network, this algorithm can be extended to detect any failure in this system (e.g., burnout of light sources or the absence of command to turn off/on at certain times of the day). Such failures, with a certain time delay, can be signaled to the competent operators and their managers. For this purpose, for example, GSM/GPRS modems can be used as a communication channel through which the system can send appropriate SMS messages of predetermined content. However, this paper will deal with more complicated situations that occur in the exploitation of PS, so this simple algorithm will not be further developed.

#### *3.2. The Detection and Analysis of the Irregular Operation of PU*

For the analysis of the irregular modes of operation of the PS, the most important period was when the monitoring of its operation began. The recording of data from the power logger started on 8 June 2009. At that time, PUs with 15 kW drive motors were installed in the PS. The control concept was still based on the "on–off" regulation of the water level in the tank, and the choice of operation of one or the other PU was performed manually, using a changeover switch with two positions. Characteristic examples of regular ACS operation that was applied in this PS at that time are shown in Figures 16 and 17. The daily diagrams of the load current changes of PU1 (27 June 2009) and PU2 (23 June 2009) are presented in these figures, respectively (data sampling period was 12 s).

**Figure 16.** The diagrams of the load current changes in PU1 recorded on 27 June 2009.

**Figure 17.** The diagrams of the load current changes of PU2 recorded on 23 June 2009.

In contrast, Figures 18 and 19 show that there were modes of operation that can be considered irregular and that could hardly have been predicted by any system modeling.

**Figure 18.** The diagrams of the load current changes of PU2 recorded on 21 June 2009.

**Figure 19.** The diagrams of the load current changes of PU2 and PU1 recorded on 25 June 2009.

The irregularity of such a system's operation is reflected in periods in which a large number of consecutive cycles of the PU switching on/off occurred. These periods are easily noticeable in the diagrams, because, due to the large number of consecutive PU2 on/off, the diagrams are not clearly visible lines, but look like colored surfaces. Both cases shown in Figures 18 and 19 relate to the operation of PU2, but such modes of PS operation were common in both PUs.

The analysis of these phenomena revealed that in such situations, the existing ACS tries many times to turn on the PU and put it into operation, while a certain protection gives the command to turn it off, not allowing it to start. The diagrams in Figure 20, which actually represent a zoomed in part of the diagram in Figure 18, show that the automatic reconnection (AR) of the PU was set to about 30 s (two starts per minute).

**Figure 20.** The diagrams of the load current changes of PU2 recorded on 21 June 2009. (Zoomed in period 16:10–16:20).

More detailed analyses have shown that the reason for this behavior of the automatic control system and protection was the dry-running protection system. At this period, this protection system consisted of a pressure switch, which under a certain pressure on the suction pipeline did not allow for the operation of the PU, giving the command to turn it off. Frequency of this phenomenon was defined by the AR system, which was set to attempt the PU reconnection about 30 s after the protection system is triggered.

Taking into account the fact that the PS Knez Selo is a booster pump station pumping water from the second to the third altitude zone, it should not be surprising that the occurrence of pressure drop on the suction pipeline that causes the pressure switch to react is relatively common.

From the diagrams shown in Figures 16–19, it can also be noticed that these phenomena are often transient (reactions of this protection system are visible, after which the local automation succeeds to turn on the PU and it continues to work for a certain period of time). However, when the pressure drop is not transient, a large number of unnecessary shutdown cycles and automatic reconnection (AR) occur. They expose the PU and additional equipment to high electro-mechanical stresses and significantly reduce their operational life. Based on the available data from June 2009, when the first data on the operation of this PU were recorded, it could be estimated that the total number of unnecessary attempts to start the PU during the month could be in the hundreds.

#### A Significance of Using Detailed Databases

It is very easy to show that the many phenomena described in the paper could not be detected or recognized using lower detailed databases. For example, if we try to use data from the existing SCADA (Supervisory, Control, and Data Acquisition) system (with sampling period of 15 min for all analog values), instead of the diagrams presented in Figures 18 and 19, we will obtain completely different diagram shapes (Figures 21 and 22, respectively).

**Figure 21.** The diagrams of the load current changes of the PU2 recorded on 21 June 2009. (Sampling period of 15 min).

**Figure 22.** The diagrams of the load current changes of PU2 and PU1 recorded on 25 June 2009. (Sampling period of 15 min).

It is obvious that based on the diagrams with a sampling period of 15 min, the described irregular PU operation could not be detected at all. Knowing that the usually used control systems in the PS do not perform data acquisition (this option is eventually performed by corresponding SCADA systems), one should not be surprised by the fact that many of the important phenomena occurring in the PS remain undiscovered.

For this reason, the presented data-driven model approach, using the previously collected detailed databases, enables the detection and analyses of various phenomena, and can be considered as one of the paper's scientific contributions. The more different the data in the databases, the more qualitative analyses can be performed.

#### *3.3. Non-Parametric Model of the Described Phenomenon (Non-Parametric Model 2)*

The previously presented description and detailed analysis of the observed phenomena are sufficient to start the design of a non-parametric model, using which the described phenomena can be modeled.

The basic idea in a non-parametric model design should be to prevent unnecessary cyclical switching on and off of the PU because they can cause failure and damage to the equipment. However, in the model with the fixed input data sequences, it is not possible to simulate the mutual effects between the input data and ACS operation. For this reason, the main purpose of designing a non-parametric model should be to enable a simulation analysis, which will show how many unwanted cycles could be reduced by using the certain control logic that will be implemented in it. A corresponding software algorithm that will apply the selected control logic and that will be able to process the input data and

give corresponding simulation results, can be actually treated as a non-parametric ("black box") model.

3.3.1. Practical Considerations about the Control Logic That Should Be Implemented in Non-Parametric Model Designing

The non-parametric model 2 can be formed based on the loading and processing of the original recorded data that will be used as the input data. The control logic that should be implemented should be based on the following practical considerations:


be assigned a value of 0 until the end of the cycle. If the total cycle time is greater than time t1, after the expiration of time t1 from the input data should be left as the originally recorded values for the first-subsequent cycle of the AR and re-shutdown, while all of the other measured values (currents, power, power factors), except for the input voltages and frequencies in the electrical network, the value of 0 should be assigned again by the end of the cycle.


#### 3.3.2. A Simulation Analysis of the Proposed Non-Parametric Model

Before the practical application of the knowledge used for the creation of the presented non-parametric model for the PS operation control, an appropriate simulation analysis should be performed.

The parametric model described in Section 2.1 was chosen to use the values of the electrical quantities related to the powers (currents) that PU drive motors take from the grid at certain moments as input data. All of these data are available and can be used as the input parameter list in the proposed non-parametric model (software designed based on the proposed control algorithm). If this model (software) uses these input data, bearing in mind that it has information on the motor load currents every 12 s, it is clear that it can detect any of the previously described cycles (protection system triggering and AR system attempt to reconnect the PU about 30 s after). The input data in the non-parametric model will be the sequences of all of the electrical parameters recorded in a certain period of time. The output data will be all of these values, changed in accordance to the previously described control logic, but only in the periods of irregular PU operation (i.e., in intervals in which all of the other measured parameters (currents, powers, power factors), except the input voltages and frequencies in the electrical network, is assigned the value 0).

The simulation results obtained by applying this model for the situations presented in Figures 18 and 19 (for the values of the parameters n = 2 and t1 = 15) are shown in Figures 23 and 24, respectively. A simple comparison showed that now the number of unnecessary PU on/off cycles had been significantly reduced.

An overview of the exact number of switching on and off in real conditions as well as after the application of the simulation model is shown in Table 1. The first column is filled using the parametric model of the system (original data), and the second using the previously described non-parametric simulation model (simulation data).

**Figure 23.** The simulation of the PU2 load current changes on 21 June 2009. (n = 2, t1 = 15 min).

**Figure 24.** The simulation of the PU2 and PU1 load current changes on 25 June 2009. (n = 2, t1 = 15 min).

From Table 1, it can be noticed that for a little more than 22 days in June 2009, the number of on/off cycles for PU1 was 209, and for PU2 425 (the system started recording the electrical energy consumption on 8 June 2009 around 18:00 h). The previously described simulation analysis showed that in this period, there should have been 138 on/off cycles for PU1 and 113 cycles for PU2 (i.e., that the unnecessary number of cycles on/off PU1 was 71, and PU2 was 312). The number of unnecessary cycles in this case was 383 and was higher than the expected number of cycles that would occur during the normal operation of PS in this period (251) by about 52.6%.

The drastic reduction in the unnecessary number of on/off cycles in both PUs, which is evident in Table 1, since July 2009, is a direct consequence of the applied knowledge obtained from previous analyses. Namely, the reference value on the pressure switch had been reduced before the end of June 2009 from 0.5 bar to a value of about 0.1 bar.

This indicates the fact that the dry-running protection system implemented by measuring the suction pressure of the PUs was not adequate and was causing shutdowns, even in situations when the pumps' impeller had been submerged.

The previously described example showed that irregular operating modes occurred in the PS operation, and that they were often unnoticed or insufficiently researched. This is especially characteristic for smaller power PSs, which are usually equipped with a simpler ACS, designed according to the expected operating modes, without connection to remote monitoring and control systems. Operators usually visit them only after certain failures appear, and leave them soon after their elimination, and this practice is common for almost all countries, regardless of their development level.


**Table 1.** The overview of the total number of switching off cycles and AR.

The specific example of the irregular operating mode of the PS described above is not an isolated example. The author's experience shows that identical or very similar situations occur in a large number of PS (especially in booster type PS), regardless of the dry-running protection type (pressure switch, pressure transmitter, conductive probes). Practice has shown that similar phenomena occur even in situations where the ACS is designed using modern equipment (motor-starters (MS) or frequency converters (FC) controlled by PLCs), which is also indicated next to the symbol for the corresponding switching elements on the physical model presented in Figure 2. This is evidenced by the large number of examples of PSs and the large amount of experimentally collected data that the authors possess.

#### *3.4. Hybrid Model for More Efficient and Effective Control in PS Operation*

When hybrid models are applied for the purpose of modeling certain processes, then the result of such models should be appropriate simulation results or appropriate predictions. The hybrid model consists of a parametric and a number of non-parametric models, whose integration should be performed in a certain way. In the example shown, it was very easy to perform this using software, because the parametric and non-parametric models used the experimentally recorded data on the electrical quantities related to the loads of the pump drive motors as input data (Figure 25).

The model integration box in this situation could actually have the same program lines as the parametric model. When the simulation results from a non-parametric model should be obtained, its output data (changed in accordance to the previously described control logic) are passed through the same procedure as in the parametric model. In this process, the corresponding nonelectrical values are obtained (pump flows, PU heads, and efficiencies). It should be pointed out that this procedure will not take place in parallel with the procedure that takes place in the parametric model. A certain time delay between these two procedures, needed for the non-parametric model procedure to perform, will have to be implemented.

**Figure 25.** The hybrid model for the tolerant control of the described irregular PS operating mode.

In Section 2.1, it was clearly shown how the input data, with the help of the selected parametric model related to the measured electrical quantities, were converted into significant information about a number of non-electrical quantities in the system. On the other hand, Sections 3.2 and 3.3 described how, based on the same set of input parameters, the occurrence of irregular PS operation (resulting in numerous negative effects for electrical and hydraulic equipment in PS and for consumers) can be detected. The performed analysis explained which logic should be used to eliminate them, and the results of the simulation analysis obtained by applying the formed non-parametric model are shown. The results present what these phenomena might look like (i.e., how the appropriate electrical quantities diagrams would look like after applying the proposed algorithms). It is clear that treating the output quantities from the non-parametric model as input parameters and passing them through the parametric model would allow us to obtain the simulation values of all of the non-electrical quantities and to determine what effects the application of such control algorithms would have on the system (by how much the number of unwanted on/off cycles would be reduced, and by how much less the amount of electricity consumed would be, etc.). This could also be easily conducted in the model integration module (Figure 25).

However, the following issues that have practical significance are much more important:


It is obvious that in the presented process of the hybrid analysis and modeling, the input data (data on electrical parameters of PUs in different operating modes) are converted into knowledge (Figure 25), which allows the described situations to be reliably recognized and to create a control algorithm that allows for more efficient and effective control of the described occurrence in the irregular operation of a PS.

Part of this process has already been presented in the previous section, and the knowledge obtained from the presented information was enough to design an algorithm for the more efficient and effective management of PS, without the need of the existing ACS to be replaced (Figure 26).

**Figure 26.** The algorithm for the tolerant control in the irregular PS operating modes.

From the control algorithm shown in Figure 26, it can be noticed that the basic idea for upgrading and improving the control function consists of the following:



#### **4. An Example of Hybrid Model Extension: Modeling New Phenomena in the PS Operation**

#### *4.1. Consequences of Long-Term Irregular PU Operation*

A large number of the dry-running pump protection system activations led to a new, more serious problem in the PS operation, which also required an adequate solution. Namely, due to a large number of protection system tripping during the month as a result of material fatigue, the pressure switch relay failures occurred, after which this protection system was out of order and the PU entered into the "dry-running" operation. This state often lasted for several hours (Figure 27). Due to worsening cooling conditions, the pump bearings are often damaged in such situations, or other mechanical damage occurs due to the excessive thermal stresses, after which the pumps must be repaired. Figure 28 clearly shows that PU2 was replaced on 8 March 2010. It is obvious that before the replacement of the pump unit PU2, which was performed on 8 March 2010, PU1 and PU2 had very similar operating characteristics, which can be seen from the fact that in that period, in normal operating modes, both PUs drew approximately equal currents (the largest number of points on the red (PU1) and blue (PU2) curves are grouped around the values of the currents in the range between 24 and 25 A). On the other hand, after the replacement of PU2, the newly installed PU drew slightly higher load currents from the network, which in operating modes was around 28 A, which indicates that it had different operating characteristics from the previous one. Additional information confirmed that the old PU was replaced by a PU

of the same type that had been repaired in the meantime. After the replacement of the PU, its performance characteristics should be changed in the parametric model, according to the procedure explained in [12]. On the other hand, the diagram of the change in the current's average value in Figure 27 shows that only 3 days after the replacement (on 11 March 2010), the repaired PU2 entered into a "dry-running" operation that had lasted for more than 3 h. This is the operating mode of pump units that should be avoided by applying an adequate protection system.

**Figure 27.** An example of the dry-running pump operation, visible through the PU2 load current change diagram.

**Figure 28.** An example of the PU2 replacement performed on 8 March 2010, visible through the load current change diagram.

The diagram in Figure 27 shows that this mode is not problematic for the pump drive motor, because the PU load currents in such modes are significantly lower than the load currents in normal operating modes. As the PU drive motor operates at a relatively low load, such operating modes cannot in any way adversely affect its operation. The main problem that prohibits the operation of pumps in such modes is that, due to a lack of water in the impeller, the cooling of the pump was much worse compared to normal operating modes. Prolonged operation in such modes can lead to increased heating of the pump shaft bearings, which in situations where they are not in good condition can lead to thermal stresses and possible damage to their housings.

Based on the results of the conducted experiments, the authors concluded that when running dry, the PUs enter into a mode where their motors draw less current from the grid, even less than the PUs' no-load currents. Measurements performed on the PU led to the conclusion that the average value of the current that the drive motors of both the PUs draw from the grid at no-load (with the submerged impeller and maximally closed valves on the thrust) was greater than 23 A.

#### *4.2. An Extended Hybrid Model for More Efficient and More Effective Control of PS Operation*

This idea can be used to form another non-parametric model that can be used to simulate the operating modes of the PS in which the failure of the PU dry-running protection system fails, which are obviously situations that can realistically occur in practice. However, it is much more practical based on this information and knowledge to design a backup protection system that will, even in the event of the failure of basic protection systems, prevent the further operation of PUs in such regimes. Through the formation of a new hybrid model and simulation analysis, it is possible to simulate its operation.

The basic idea for such a backup protection system can be described as follows: As soon as the system detects that the average load current of the drive motor is less than 23 A at several consecutive points (e.g., 5 points, which corresponds to a period of 60 s), it gives a signal to turn off the PU after this period. This command in the system has the treatment of an identical command issued by a pressure switch under normal circumstances. This is enough to form the non-parametric model 3, an extended hybrid model corresponding to this situation (Figure 29).

**Figure 29.** The extended hybrid model of the serial-parallel structure.

As in the previously described case, this non-parametric model would use the same input parameters and its integration into a more advanced hybrid model would not be complicated again.

The simulation result of the load currents' average values of the PU2 drive motor, obtained after the application of the extended hybrid model, is shown in Figure 30.

**Figure 30.** The simulation of the PU2 load current average values obtained by the application of the extended hybrid model.

The same idea was used to develop an algorithm for the tolerant control after pressure switch failure (Figure 31), and to implement it in a backup dry-running pump protection system.

**Figure 31.** The algorithm for the tolerant control after a pressure switch failure.

The presented model is just an obvious example of a tolerant system control after a fault occurs. In this case, the application of non-parametric (i.e., hybrid models) enabled the realization of a more advanced ACS that has the possibility of tolerant management in the case of the failure of one protective component (in this case, the pressure switch, which served to protect the PU from dry-running).

In a similar way, non-parametric and hybrid models can be realized that will enable the tolerant control of PS in the case of the failure of other protective devices (voltage protection, current protection, overload protection, frequency protection, etc.). The presented model can be applied in all PS booster types because the described phenomena occurred in all similar PS.

The authors designed and suggested the implementation of this backup protection system in March 2010. The system was put into operation in September 2013 at PS Knez Selo after obtaining their permission for its practical implementation. Since then, the backup protection system has been tested in the operation of PS Knez Selo for many years, and its efficiency has been confirmed in the conducted experiments as well as in practice.

#### **5. Conclusions**

In the hybrid modeling theory, the combined use of physics-based and data-driven models is of great importance. In this paper, some of the numerous advantages of datadriven models and the hybrid analysis and modeling approach have been demonstrated with the example of one booster PS (the PS Knez Selo).

The authors used and processed the available historical data on the system's operation, in order to discover their mutual interrelationships and to identify different unexpected phenomena, irregular modes, or faults that occurred in the PS, especially those that previously could not be detected. The presented results pointed to a significance in the detailed database usage and showed that many phenomena described in the paper could not be either recognized or detected at all when using lower detailed databases. This was practically demonstrated by the corresponding comparison between the results taken from the existing SCADA database (with a sampling period of 15 min for all of the analog values) and the database previously collected by the authors (with the applied sampling period of 12 s).

The presented comparisons demonstrated that using diagrams with a sampling period of 15 min, the described irregular PU operation could not be detected at all. Knowing that the usually used control systems in PS do not perform data acquisition (this option is eventually performed by corresponding SCADA systems), one should not be surprised by the fact that many of the important phenomena occurring at the PS remain undiscovered. For this reason, the presented data-driven modeling approach, using previously collected detailed databases, enables the detection and analyses of various phenomena that could not be detected before, and that can be considered as one of the paper's scientific contributions. In several examples presented in the paper, it has also been shown that the more different the data in databases that are recorded, the more qualitative the analyses that can be performed.

Relating to the phenomena that could not be detected before, two of them were investigated in the paper.

The first case was the irregular PU operation that was reflected in many unnecessary on/off cycles. The analysis of these phenomena revealed that in such situations, the existing ACS tries many times to turn on the PU and put it into operation, while the dry-running protection system gives the command to turn it off, not allowing it to start. The frequency of this phenomenon is defined by the AR system, which is the part of the ACS and that was set to attempt PU reconnection about 30 s after the protection system was triggered.

The simulation results obtained by the corresponding hybrid model showed that for a little more than 22 days in June 2009, the total number of unnecessary and unwanted on/off cycles for both PUs was 383, which was higher than the expected number of cycles that would occur during the normal operation of PS in this period (251) by about 52.6%. The simulation results also showed a drastic reduction in the unnecessary number of on/off cycles in both PUs since July 2009. This was a direct consequence of the applied knowledge obtained from the presented analyses. Namely, before the end of June 2009, the authors suggested that the reference value on the pressure switch should be reduced from 0.5 bar to a value of about 0.1 bar, which was implemented.

This analysis directly pointed out the fact that the dry-running protection system based on measuring the suction pressure of PUs was not adequate and was causing shutdowns even in situations when the pumps' impeller had been submerged.

Without the information and knowledge acquired during the presented research, the authors supposed that this specific phenomenon could not be discovered. The best evidence that this assumption is real represents the fact that even today, many PS designers apply a dry-running protection system based on pressure monitoring. The results of the research presented in this paper undoubtedly confirm that this is not an adequate solution at all.

The second specific case analyzed in the paper was the failure of the dry-running protection system (caused by the actuator relay failure due to the material fatigue) occurred in March 2010. After this failure, the protection system was completely out of order, and the PU entered into the dry-running operation modes.

Based on the analyses presented in the paper, the authors designed and suggested the implementation of a backup protection system in March 2010. This protection system was based on monitoring the drive motor load currents, which became lower than some of the pre-set reference values only in the PU dry-running operation modes. This backup protection system was put into operation in September 2013 at PS Knez Selo after obtaining their permission for their practical implementation. Since then, it has been tested in the operation of PS Knez Selo for many years, and its efficiency has been confirmed in the conducted experiments as well as in practice.

According to the authors' opinions, it should be pointed out that the examples presented in the paper undoubtedly show that a data-driven approach to hybrid analysis and modeling is indeed an approach that has huge potential. Generally speaking, it is shown that it is possible to obtain information on the operation of other PSs in the system or other parts of the system based on the data on the operation of the subjected PUs. In the related example is shown that the leaps in the active power observed in the corresponding diagram in periods when some of the PUs in the PS Knez Selo operate can be converted into useful information about the ons/offs of the PUs installed at PS Matejevac 2, which supplies the lower altitude zone as well as PS Knez Selo. These possibilities provide enormous opportunities for various innovations and improvements in the quality of the system operation in regular or irregular conditions, malfunctions, or various emergency situations.

Although the complete analysis presented in this paper refers to one PS, the conclusions are of general importance because similar phenomena occur in all booster PSs.

According to the authors, one of the most important conclusions, indirectly indicated by this paper, is the high importance of continuous data recording on the operation of controlled facilities (in this case, the PU in the PS). The possibilities of data recording with a sampling period of a few seconds are especially important because the data collected in this manner enables the detection of numerous specific phenomena during the operation of the system.

Hybrid modeling techniques enable the application of various innovations that lead to the reduction in system exploitation and maintenance costs.

These types of applications could also have importance for the future development of "smart grid" and "smart city" projects.

**Author Contributions:** Conceptualization, M.M., M.R., M.R.-A., D.L. and Z.S.; Methodology, Z.S.; Software, M.M.; Validation, M.R., M.R.-A. and Z.S.; Formal analysis, D.L.; Investigation, M.R.-A.; Resources, Z.S.; Data curation, M.M.; Writing—original draft preparation, M.R.-A. and Z.S.; Writing—review and editing, M.M.; Visualization, M.M.; Supervision, D.L. and Z.S.; Project administration, M.R. and Z.S.; Funding acquisition, Z.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This paper was supported by the Ministry of Education, Science, and Technological Development of the Republic of Serbia through the following projects and/or contracts (Projects: NP EE413-118A, NPV-19.B, 401-00-00144/2008-01-IP, III 44006; Contract No. 451-03-68/2022-14/200102, 4 February 2022).

**Data Availability Statement:** The data presented in this study are available upon request from the submitting author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

