The automation adopted in the San Giovannello station operates with a very small control time step, i.e., one second, potentially resulting in a high number of setting corrections, and therefore in high wear of the actuators. Furthermore, it assumes that the maximum hydropower is always attained at the minimum downstream pressure. To overcome these limitations, an RTC based on the digital twin of the San Giovannello station was constructed. This digital twin includes the hydraulic model of the plant and, given the valve and the flap configuration, is able to accurately reproduce the behaviour of the San Giovannello station equipped with the PRS turbine, in terms of both hydraulic variables and hydropower generation.

The layout of the proposed digital twin is represented in Fig. 4, Table 1 and Table 2. Valves V1, V2, V3 and the turbine PRS are represented as fictitious pipes with fully turbulent flow and a specific resistance *r**k*, that is:

\(\left({H}_{i}-{H}_{j}\right)=\left({P}_{i}-{P}_{j}+{D}_{k}\right)={r}_{k}{Q}_{k}^{2}\) (1),

where *H* is the nodal total head, *P* is the pressure head, *D* is the elevation difference across the pipe, *Q* is the flow rate, and *k* is the index used for the fictitious pipe connecting upstream node *i* to downstream node *j*. The index *k* also applies to the valve or turbine present in the pipe. *O**1* is the opening setting in the by-pass valve V1; *O**2* and *O**3* are the opening settings of valves V2 and V3, respectively; finally, O4 corresponds to λ. For the valve-fitted pipes, the correlations observed in the data confirmed the resistance *r**k* to be a function of the actuator opening setting, i.e. *r**k*(*O**k*) for *k* = 1,…, 3. For the turbine-fitted pipe, the resistance *r*4 was assumed to be a function of both *O*4 and *Q*4, i.e., *r*4(*O**4*, *Q**4*), based on the correlations observed in the data. The relationships for evaluating *r**k* can be constructed starting from the historical data measured after a preceding time horizon (e.g. within the previous 10 days). The changing boundary conditions are the demand in node 5, corresponding to Qw1 plus Qw2 in Fig. 2, as well as the pressure head in node 2, which is assumed to be equal to the measured pressure head PU and not to be affected by the valve and turbine configuration changes. The same hypothesis is applied to node 6, where the pressure head acquired by PA is assumed as a/the boundary condition. This choice is motivated by the fact that the Raganzilli and Cava Ricevuto tanks, as well as their inlet mains, are not managed by Siciliacque SpA, and at least one regulation valve, with an unknown setting, is included inside.

Table 1

San Giovannello digital twin network nodes.

Node | 1 | 2 | 3 | 4 | 5 | 6 |

*H* m a.s.l. | 75 | Time function (measured) | - | - | - | Time function (measured) |

Elevation *Z* m a.s.l. | - | 88 | 88 | 88 | 37 | 37 |

Demand *Q* | QT2 | - | - | - | Time function Qw | QT1 |

Pressure Sensor | - | PU | - | - | - | PA |

Table 2

San Giovannello digital twin network pipes.

Pipe | 1 | 2 | 3 | 4 | 5 | 6 |

Upstream node *i* | 2 | 4 | 5 | 2 | 4 | 3 |

Downstream node *j* | 3 | 1 | 6 | 3 | 5 | 4 |

*r* | *r*1(*O*1) | *r*2(*O*2) | *r*3(*O*3) | *r*4(*O*4,*Q*4) | - | - |

The following optimization problem is embedded into the digital twin:

Maximize \(PW={Q}_{4}\left({P}_{2}-{P}_{3}+{D}_{4}\right)\eta\) (2)

subject to the following constraints:

*P* 3 ≥ *P**set*; *Q*3 = *Q**T*1; *Q*2 = *Q**T*2; (3)

0 ≤ *Ο**k* ≤ 100, *k* = 1,2,3; 30 ≤ *Ο*4 ≤ 100, (4)

*P* 2 = *PU, P*6 *= PA, P*1 *=* 0 (5)

where *Q* and *P* are the numerical model pipe flow rates and nodal pressure heads, *P**3*, *P**5*, *Q**1*, *Q**2*, *Q**3*, *Q**4* are functions, through the same model, of the assigned boundary conditions of the model and of the decision variables *O**1*, *O**2*, *O**3* and *O**4*.

The constraint on *O*4 is due to a turbine constructive requirement on the flap setting *λ*. The decision variables of the optimization problem are the four configuration-settings *O* of the three valves and of the turbine flap, which affect all the digital twin variables. The variable *η* is equal to the total (mechanic + hydraulic + electric) efficiency of the turbine, multiplied by the gravity acceleration. It is a function of the flap opening setting *O*4 and of the velocity ratio Vr, which is a function of flow rate *Q*4 and of λ. Therefore, in this work, it was assumed that *η = η*(*O*4,*Q*4).

Based on the assumption that the boundary conditions assigned to the model are not significantly affected by the valve and flap opening changes, the relationships *r*1(*O*1), *r*2(*O*2), *r*3(*O*3), *r*4(*O*4,*Q*4), as well as *η*(*O*4,*Q*4), were obtained by means of data-driven empirical relationships, to be calibrated based on recent measurements.

The optimization problem presented in equations (2–5) was solved by using a specifically developed algorithm, which includes a hydraulic solver (Creaco et al., 2022a) for enforcing preservation of energy balance on pipes and mass conservation at nodes. The digital twin was constructed in the Matlab2022b® environment and is linked to EPANET 2.2 code (Rossman et al., 2020) for visualizing the results.

The operation of the digital twin for the ARTC of the San Giovannello station can be explained as follows (see Fig. 5). The ARTC is applied at the end of each control time step Δ*t* by optimizing the actuator settings with the digital twin, given the measured boundary conditions and using the relationships *r*1(*O*1), *r*2(*O*2), *r*3(*O*3), *r*4(*O*4,*Q*4), and *η*(*O*4,*Q*4) calibrated with the measures collected in a previous *T**precal* time window. The time step Δ*t* must be small enough to follow the dynamics of the system and large enough to dampen the flow transient effects in the system and the accidental errors in the measurements. If the actuator setting variations proposed by the optimizer at the generic time instant *t* are smaller than 1%, the settings of the previous time steps are kept for the new *Δt*.

The calibration of the relationships *r*1(*O*1), *r*2(*O*2), *r*3(*O*3), *r*4(*O*4,*Q*4), and *η*(*O*4,*Q*4) is carried out at the end of each *T**CAL* time window using the measures collected in the previous *T**PRECAL* time window.