**Reduced-Scale Models of Variable Speed Hydro-Electric Plants for Power Hardware-in-the-Loop Real-Time Simulations**

#### **Baoling Guo 1,2,\*,†, Amgad Mohamed 3, Seddik Bacha 2, Mazen Alamir 3, Cédric Boudinet <sup>2</sup> and Julien Pouget <sup>1</sup>**


Received: 30 September 2020; Accepted: 29 October 2020; Published: 3 November 2020

**Abstract:** Variable Speed Hydro-Electric Plant (VS-HEP) equipped with power electronics has been increasingly introduced into the hydraulic context. This paper is targeting a VS-HEP Power Hardware-In-the-Loop (PHIL) real-time simulation system, which is dedicated to different hydraulic operation schemes tests and control laws validation. Then, a proper hydraulic model will be the key factor for building an efficient PHIL real-time simulation system. This work introduces a practical and generalised modelling hydraulic modelling approach, which is based on 'Hill Charts' measurements provided by industrial manufacturers. The hydraulic static model is analytically obtained by using mathematical optimization routines. In addition, the nonlinear dynamic model of the guide vane actuator is introduced in order to evaluate the effects of the induced dynamics on the electric control performances. Moreover, the reduced-scale models adapted to different laboratory conditions can be established by applying scaling laws. The suggested modelling approach enables the features of decent accuracy, light computational complexity, high flexibility and wide applications for their implementations on PHIL real-time simulations. Finally, a grid-connected energy conversion chain of bulb hydraulic turbine associated with a permanent magnet synchronous generator is chosen as an example for PHIL design and performance assessment.

**Keywords:** hydro-electric plant; variable speed operation; 'Hill Charts'; reduced-scale model; power hardware-in-the-loop; real-time simulation; testing and validation

#### **1. Introduction**

Hydropower, as a type of cost-competitive renewable source, has played a significant role in electricity mix [1,2]. In recent years, hydropower contributes to balancing the intermittent electricity resources [3,4]. However, the conventional hydro-electric plants do not perfectly meet these new demands because of its slow reaction dynamics. The power adjustment dynamics can be enhanced by introducing Power Electronics (PE) converters [1]. Moreover, the efficiency of hydraulic energy conversion can be improved through the variable speed operation when the hydraulic conditions change frequently [5,6]. In addition, the previous research indicates that variable speed operations can alleviate water hammer effects and optimize various hydraulic transients [2]. Therefore, a Variable Speed Hydro-Electric Plant (VS-HEP) equipped with PE has been increasingly introduced into the hydraulic context [1,2,6–8]. Limited to the PE capacity, the variable speed technique is mainly implemented into small hydro-electric plants [9–13]. The Run-of-River (RoR) scheme is economic and environmentally friendly as no dam is required, which have been suggested for small hydropower applications [9]. In an RoR scheme, a class of hydraulic turbines such as semi-Kaplan or Kaplan [5,10,13], propeller [11], and bulb turbine [12] has become the favourable prime mover because of their high efficiency under low water heads and fast flow-rates. In the work, to highlight, we will therefore focus on a class of small-scale hydraulic turbines.

VS-HEPs present many advantages in hydraulic applications but bring problems of the controls, disturbances rejections, and induced constraints as well [7,8]. Advanced control techniques continue to be needed to achieve higher robustness and enable faster control dynamics for VS-HEPs [14–17]. Reduced-scale PHIL benchmarks become necessary to test and validate control laws in hydraulic schemes because repetitive deterministic experiments cannot be conducted under a real disturbed environment [18–20]. The key question of designing a successful VS-HEP PHIL real-time simulation system is how to build a proper mathematical model for the concerned hydraulic turbine. However, unlike well-developed wind turbines [19], there is no general efficiency expression for hydraulic turbines. The hydraulic efficiency depends on both the turbine geometrical design and its practical operating conditions [21]. Various models have been proposed in previous works. The Kaplan turbine is described by a linear decreasing line of the mechanical torque versus the turbine speed in [10], while both the pitch adjustment and the guide vane opening are not considered. The hydraulic losses model of a bulb turbine is approximated by a 2D expression of the flow rate and the turbine speed in [12]. In another work [11], a 2D efficiency model of a propeller turbine is achieved by numerical approximations as well. The 2D models generally fix the guide vane angle or neglect its effects on the obtained hydraulic models. The real-time efficiency model is considered in [13]; however, precise measurements of water flow rate are still a critical challenge. A group of static data, measured from a real-life full-scale hydraulic electric plant, are then discretely stored in the look-up tables in [22,23]. If a control method works on the condition that the system must be continuous, such as the optimal control (Casadi, IPopt, etc.), a continuous model would be required [24]. Then, a class of look-up table based models is not able to be employed in such control schemes. A class of nonlinear or linear dynamic hydraulic models has been discussed in [25,26]. They are commonly used for hydraulic dynamics control test under offline simulations environments. However, the computational burden would increase by introducing large hydraulic dynamics. If the digital processor in use cannot satisfy the computational rapidity for real-time simulations, the modelling complexities of hydraulic turbine have to be sacrificed. Otherwise, the behaviours of implemented controllers could be deteriorated [19]. Hydraulic features are regressed to be singly linked to the rotational speed in [27], which is particularly suitable for a propeller turbine with a fixed guide vane angle. In a recent work of [28], the water flow rate is approximated by a third-degree polynomial, and the efficiency model is determined by the Multi-Layer Perception (MLP) of a neural network. However, the model determination process based MLP is complicated; meanwhile, a large scale of measurements are required to ensure the efficiency accuracy [28]. The Euler turbine equations are considered in [17,29]. The input is the turbine's main geometry at the best efficiency point of operation. However, the performance of hydraulic turbine can be affected by the hydraulic conditions of its localisation; moreover, the system efficiency feature can also change due to the ageing [28]. The discussions above indicate that a practical and generalised modelling approach regarding various hydraulic turbines is still required for achieving efficient PHIL real-time simulations.

In this paper, a generalized hydraulic modelling approach is presented for achieving efficient PHIL real-time simulations. The static models are based on the commonly used 'Hill Charts' measurements provided by industrial manufacturers. The regression models of water flow rate and hydraulic efficiency are computed by using optimization routines. The optimization routine is defined with the objective of minimizing the difference between the measured data points and the determined regression models. In addition, we can update the models with recent 'Hill Charts' measurements to avoid the model deviations due to the ageing or localisations. Furthermore, this work introduces the nonlinear dynamic model of the guide vane actuator. The dynamic model helps to more accurately evaluate how the hardware parts react to the induced dynamics and what would be the impacts on the whole closed-loop control performances. The scaling rules of hydraulic turbine and the mechanical drive train are described as well. Then, reduced-scale models can be flexibly obtained according to various laboratory conditions. Some features are enabled for its implementations on PHIL real-time simulations: decent accuracy, high computation efficiency, significant flexibility, and wide applications. Furthermore, a bulb hydraulic turbine associated with a Permanent Magnet Synchronous Generator (PMSG) is chosen as an example for the performance assessment. Based on the flexible hybrid PHIL benchmark in the laboratory [19], a VS-HEP experimental benchmark is built, which is adapted to the obtained reduced-scale turbine model.

This rest of this paper is organized as follows: A 'Hill Charts' based modelling approach is presented in Section 2, where a bulb turbine with adjustable guide vane actuator is chosen as an example to assess the performance. In addition, scaling rules of hydraulic turbine and the mechanical drive train have been described in this section as well. In Section 3, a VS-HEP PHIL real-time simulation system is developed that is adapted to the obtained reduced-scale turbine model. Experiments have been conducted to assess the performance of the proposed reduced-scale models under the built PHIL simulation benchmark in Section 4. Remarks are concluded in Section 5.

#### **2. 'Hill Charts' Based Hydraulic Modelling**

In this section, the static features are firstly modelled based on the 'Hill Charts' data. The dynamic model of guide vane actuator is then considered and described. Moreover, the scaling rules regarding both the hydraulic turbine and the mechanical drive train are discussed as well. Finally, the modellings discussed in Section 1 and the proposed hydraulic model in this work for their implementations on PHIL simulations are generally compared and analysed. The following assumptions are made before the modelling:


#### *2.1. Static Features Modelling*

The general procedure of obtaining the reduced-scale hydraulic model is schematically described in Figure 1, which will be detailed in the following parts. This method differs from the scenario discussed in [34], as an actual full-scale 'Hill Charts' data of a bulb turbine is used in computing the mathematical model, without the need of removing the unstable region (S-region) [35]. In addition, being different from the modelling approach proposed in [28], the models of water flow rate and efficiency are regressed through solving optimization problems.

**Figure 1.** Procedure of a 'Hill Charts' based modelling.

#### 2.1.1. Static Model Generation

Hydraulic turbines conventionally use the 'Hill Charts' (see Figure 2) to describe the steady-state relationships of the turbine speed *N* (r/min), the guide vane opening *γ* (ratio) if it exists, the water head *H* (m), the water flow rate *Q* (m3/s), and the mechanical torque *T* (Nm). The unitary variables *N*11, *Q*11, and *T*<sup>11</sup> are generally used to represent the 'Hill Charts' in the hydraulic context. The 'Hill Charts' measurements are usually provided by the industrial manufacturers.

**Figure 2.** Diagram of 'Hill Charts' and the variable speed operation, *η* represents the turbine efficiency.

The variables *N*11, *Q*11, and *T*<sup>11</sup> with turbine's diameter *D*(*m*) for each guide vane opening *γ* can be derived as follows [35]:

$$N\_{11} = \frac{ND}{\sqrt{H}} ; \ Q\_{11} = \frac{Q}{D^2 \sqrt{H}} ; \ T\_{11} = \frac{T}{D^3 H} \tag{1}$$

*Energies* **2020**, *13*, 5764

Based on the data points (*Q*11, *T*11, *N*11, *γ*) provided by the industrial manufacturer, the regression models *Q*˜ <sup>11</sup> and *T*˜ <sup>11</sup> are defined by

$$\mathcal{Q}\_{11}(B\_{\mathbb{Q}'}\gamma, N\_{11}) = A\_{\mathbb{Q}}(\gamma)B\_{\mathbb{Q}}\mathbb{C}\_{\mathbb{Q}}(N\_{11})^T \tag{2}$$

$$\mathcal{T}\_{11}(B\_T, \gamma, N\_{11}) = A\_T(\gamma) B\_T C\_T (N\_{11})^T \tag{3}$$

where:

$$A\_Q(\gamma) = \begin{bmatrix} 1 & \gamma \dots \gamma^{m\_Q} \end{bmatrix} \in \mathbb{R}^{1 \times (m\_Q + 1)}, \ A\_T(\gamma) = \begin{bmatrix} 1 & \gamma \dots \gamma^{m\_T} \end{bmatrix} \in \mathbb{R}^{m\_T + 1}$$

$$\mathbb{C}\_Q(N\_{11}) = \begin{bmatrix} 1 & N\_{11} \dots N\_{11}^{n\_Q} \end{bmatrix} \in \mathbb{R}^{1 \times (n\_Q + 1)}, \ C\_T(N\_{11}) = \begin{bmatrix} 1 & N\_{11} \dots N\_{11}^{n\_T} \end{bmatrix} \in \mathbb{R}^{n\_T + 1}$$

where the parameters *mQ*, *nQ*, *mT*, and *nT* need to be chosen. The optimization routine *lsqlin* defined in Matlab is chosen to solve the related optimization problems (4) and (5) to get the regression parameters *BQ* <sup>∈</sup> <sup>R</sup>(*mQ*+1)×(*nQ*+1) and *BT* <sup>∈</sup> <sup>R</sup>(*mT*+1)×(*nT*+1). By solving the least square optimization problems (4) and (5), we minimize the mismatch between the chosen regression model and the provided data points, while making sure that, whenever *γ* = 0, we get *Q* = 0 regardless of the rotational speed, in addition to ensuring that *T* = 0 when *γ* = 0 and *N* = 0:

$$\begin{aligned} \min\_{B\_Q} & \quad \frac{1}{2} \left[ Q\_{11} - \tilde{Q}\_{11} (B\_{Q'}, \gamma\_\prime N\_{11}) \right]^2 \\ \text{s.t.} & \quad \tilde{Q}\_{11} (\gamma = 0) = 0 \end{aligned} \tag{4}$$

$$\begin{array}{ll}\underset{B\_{T}}{\min} & \frac{1}{2} \left[ T\_{11} - \tilde{T}\_{11} (B\_{T\_{\prime}} \gamma\_{\prime} N\_{11}) \right]^{2} \\ \text{s.t.} & \mathcal{T}\_{11} (\gamma = 0, N\_{11} = 0) = 0 \end{array} \tag{5}$$

The choice of the regression parameters *mQ*, *nQ*, *mT*, and *nT* is made such that the relative regression errors can be minimized; meanwhile, the high sensitivity due to over-fitting must be avoided. If higher order polynomials are considered, then we may get undesirable oscillatory behaviours, since it does not represent the nature of the measured data used for regressions and causes higher regression errors. The procedure is as follows: firstly, the parameters *mQ*, *nQ*, *mT*, and *nT* can be gradually increased; the regression errors are then computed; the optimal values of *mQ*, *nQ*, *mT*, and *nT* are chosen just before the data points start to be over-fitted, which in turn can cause oscillatory models. Figures 3 and 4 show the comparative curves the original unitary industrial data points and the generated regression models with the defined parameters *mQ* = 3, *nQ* = 2, *mT* = 2, and *nT* = 5.

**Figure 3.** *Q*˜ <sup>11</sup> vs. *N*<sup>11</sup> for different values of *γ*.

**Figure 4.** *T*˜ <sup>11</sup> vs. *N*<sup>11</sup> for different values of *γ*.

Figures 5 and 6 show the relative regression errors of *Q*<sup>11</sup> and *T*11. The relative regression errors are always less than 5% in the range of interest. These optimization models can provide satisfied regression results due to the fact that the main patterns in the 'Hill Charts' have been captured by the chosen functions.

**Figure 5.** Relative error of *Q*<sup>11</sup> regression (%).

**Figure 6.** Relative error of *T*<sup>11</sup> regression (%).

#### 2.1.2. Efficiency and Mechanical Power

The mechanical torque and the flow rate can be derived by combining Equations (2) and (3) and the definitions of the unitary variables (1) to get:

$$Q(\gamma, \omega, H, D) = \tilde{Q}\_{11} D^2 \sqrt{H} \tag{6}$$

$$T(\gamma, \omega, H, D) = \tilde{T}\_{11} D^3 H \tag{7}$$

Once the water head *H* and the turbine's diameter *D* are specified, Equations (2)–(7) are then used to compute *Q* and *T* for each given *N* and *γ*. The efficiency *η* and the mechanical power *Pm* (W) are then computed using Equations (8) and (9), for *H* = 1 m and *D* = 0.25 m, which are shown in Figures 7 and 8, respectively:

$$\eta(\gamma,\omega,H,D) = \frac{T(\gamma,\omega,H,D)N}{\rho g H Q(\gamma,\omega,H,D)}\tag{8}$$

$$P\_m(\gamma, \omega, H, D) = \eta(\gamma, \omega, H, D)\rho \boxtimes \mathbb{Q}(\gamma, \omega, H, D) \tag{9}$$

where *ρ* (kg/m3) is the volume density of water, *g* (m/s2) is the gravitational acceleration, and *ω* = <sup>2</sup>*π<sup>N</sup>* 60 (rad/s) is the rotational speed.

The water head is supposed to be fixed with respect to the derived static hydraulic features *η*(*γ*, *ω*, *H*, *D*) and *Q*(*γ*, *ω*, *H*, *D*). The reason and its rationality have been analysed in the beginning of Section 2. However, the water head variations are essential to be considered for the hydraulic power computation regarding the variable speed control. The expression (9) is then modified to (10), in which the water head variations Δ*H* are taken into account:

$$P\_{\mathfrak{m}}(\gamma,\omega,H,D) = \eta(\gamma,\omega,H,D)\rho\mathfrak{g}(H+\Delta H)\mathbb{Q}(\gamma,\omega,H,D) \tag{10}$$

**Figure 8.** *Pm*(*γ*, *ω*, *H*, *D*) where negative values are set to zero.

This section provides a flexible philosophy on obtaining mathematical models for a class of hydraulic turbines with 'Hill charts' efficiency characteristics. Different definitions of *Q*˜ <sup>11</sup> and *T*˜ 11 can be adapted depending on the provided 'Hill Charts' data. The application scope of the proposed modelling approach is thus expanded.

#### *2.2. Guide Vane Actuator Dynamic Model*

The work in [34] assumes that the guide valve can be opened and closed almost instantaneously. However, it makes sense to consider the guide valve acting process in order to replicate the hydraulic dynamics. The introduced dynamics can help more precisely evaluate how the physical hardware reacts to the dynamics and what would be the impact on the electric control performances [19].

The full nonlinear dynamic model of the guide vane actuator used at GE is schematically presented in Figure 9, which is introduced in this work. The saturation blocks are used to impose limits on the guide vane angle *γ* and the change rate of the guide vane *γ*˙ , respectively. This model can correctly replicate the nonlinear dynamic behaviours of the guide vane actuator used for a bulb turbine.

**Figure 9.** Actuator's schematic, *Ts* and *Td* are time constants of the chosen actuator, *γ<sup>c</sup>* represents the reference of guide vane opening.

#### *2.3. Scaling Operations*

#### 2.3.1. Scaling of Hydraulic Turbine

The full size 'Hill charts' data provided by GE have been used to generate the models; these data are related to a 59 MVA bulb turbine rotating at around 8.5 rad/s. The 'Hill charts' can be generated under different working conditions by using the homothetic scaling laws [36,37]. In order to adapt to the laboratory experimental conditions, the diameter of the hydraulic turbine and the water head both can be adjusted, and the resulting regression models can be rescaled as follows:

$$
\mathcal{Q}\_{11}^s(\gamma\_\prime \mathcal{J} \times \mathcal{N}\_{11}) = a \times \mathcal{Q}\_{11} \tag{11}
$$

$$
\hat{T}\_{11}^s(\gamma\_\prime \zeta \times \mathcal{N}\_{11}) = \beta \times \tilde{T}\_{11} \tag{12}
$$

where the scaling coefficients *α*, *β*, *ζ* > 0.

#### 2.3.2. Scaling of Drive Train

The hydraulic turbine is coupled with the PMSG through the drive train. The dynamic equation between the hydraulic driven torque *T*(*γ*, *ω*, *H*, *D*) and the electromagnet torque *Te* (Nm) of PMSG can be formulated by

$$J\frac{d\omega}{dt} = T\_t - T(\gamma\_\prime \omega, H, D) - B\omega \tag{13}$$

where *<sup>J</sup>* (kg·m2) is the total inertia, and B (Nm·s) is the friction factor [14].

A Direct Current Motor (DCM) is used to emulate the hydraulic torque in response to the turbine speed (detailed in Section 3). Some similitude concerns must be confirmed to precisely emulate the turbine dynamic model [38]. The original inertia *J* in Equation (13) needs to be corrected for the emulated system. The scaling rules for a category of turning systems can be found in [38]. The launching time *Ts* (s) of the turbine emulator is required to be identical to that of the desired turbine as given by

$$T\_s = \frac{\text{J}\omega\_n^2}{P\_n} = \frac{\text{J}\_\epsilon \omega\_{n\epsilon}^2}{P\_{n\epsilon}} \tag{14}$$

where *Je* (kg·m2) is the emulator's inertia, *Pn* (W), and *Pne* (W) are the rated power, *<sup>ω</sup><sup>n</sup>* (rad/s), *<sup>ω</sup>ne* (rad/s) are the rated rotational speed, and index '*n*' and '*ne*' indicate rated parameters of the original turbine and its emulator, respectively.

The inertia *Je* employed to compute the emulated inertial torque can be corrected as the following:

$$J\_e = J \cdot m^2 / n \tag{15}$$

where the scaling factors are respectively given by *m* = *ωn*/*ωne* and *n* = *Pn*/*Pne* [38].

#### *2.4. General Comparison with Previous Modellings*

The hydraulic models discussed in Section 1 and the proposed model regarding their implementations on PHIL simulations are generally compared in Table 1, which are discussed from the following aspects: accuracy, computational rapidity, and scalability.

**Table 1.** Comparisons of different modellings for their implementations on PHIL benchmarks.


#### 2.4.1. Accuracy Analysis

Compared to the linear torque model [10] and 2D hydraulic models [11,12], more effectors including water flow rate *Q*, guide vane ratio *γ*, water head *H*, and turbine speed *ω* have been taken into consideration. In addition, we introduce the dynamic model of guide vane actuator, which can help more definitely evaluate the induced effects by the dynamics on the whole control system. Note that the efficiency feature of hydraulic turbine can deteriorate due to the ageing [28]. Then, we can update the models with recent 'Hill Charts' measurements.

#### 2.4.2. Computational Rapidity Analysis

Transient dynamics are usually considered for various offline simulation models for hydraulic dynamics study in [25,26]; however, the increased computational burden would be a challenge for real-time simulations. If the execution time devoted to the real-time simulation model exceeds the sampling period, the control system behaviour can be deteriorated. Consequently, the complexity of hydraulic model has to be sacrificed in order to achieve a rapid computational speed. In fact, the hydraulic transient dynamics are much slower compared to that of the electrical control; therefore, there are no great effects induced when the transient dynamics are neglected.

#### 2.4.3. Scalability Issues

The scalability refers to the analyticity, the flexible scaling, and the application field. In this work, an analytical modelling approach is proposed by taking use of optimization routines. By following the similarity laws, reduced-scale models are flexibly established for various laboratory operation conditions, the modelling flexibility is thus improved. In addition, this approach is based on 'Hill Charts' data, which can be adapted to various types of hydraulic turbines but not be limited to a particular turbine type. Furthermore, either a continuous model or a discrete model can be chosen, which depends on its application fields.

It can be remarked that the proposed modelling approach enables positive features for their implementations on VS-HEP PHIL real-time simulations systems: decent accuracy, light computational complexity, and high scalability.

#### **3. PHIL Real-Time Simulation System**

Based on the flexible hybrid PHIL benchmark in the laboratory [19], a VS-HEP test rig is built as shown in Figure 10, which is adapted to the obtained reduced-scale hydraulic model. A systematic PHIL design procedure implemented to wind power systems is described in [19], whose main design steps can be followed. The global design schematic is illustrated in Figure 10a. The experiential benchmark under test is shown in Figure 10b. The primary elements are:


**Figure 10.** PHIL-based hybrid real-time VS-HEP simulation system, (**a**) global design diagram, with *Thdy* the emulated hydraulic torque, *Km* the torque constant, *iDCM* the direct current of DCM, (**b**) PHIL experiment benchmark

Furthermore, the detailed descriptions of the main subsystems illustrated in Figure 10 will be respectively provided in the following subsections.

#### *3.1. dSPACE Modular, dSPACE Controller, and dSPACE Environment*

The dSPACE as a rapid-prototyping system allows the implementation of numerical models and any control structures for energetic conversion systems [39,40]. In this work, the model of hydraulic turbine and the control diagram are both achieved with the dSPACE real-time digital simulator based on the MATLAB/Simulink environment. The processor of the dSPACE card is dedicated to compute the model and synchronised on a timer; this allows for running the model in real time [41]. Moreover, some extra hardware interfaces are provided to facilitate the design of hardware prototype, such as Analog to Digital (ADC), Digital to Analog (DAC), Pulse Width Modulation (PWM), Encoder, Resolver, etc. In addition, dSPACE provides add-in toolboxes to Matlab/Simulink software which allows for easy setting of interface cards and easy development of control system in a Simulink environment [39,40].

The dSPACE system under use is composed of a processor card and a number of input/output cards. From the dSPACE host PC, the system can convert the Simulink models into real-time executable codes via the Real-Time Workshop. Then, the codes are executed on one or several processors [41]. The RTI dSPACE 1005 system is currently used in the laboratory, which consists of a DS1005 PowerPC processor card (as the motherboard) that operates at 480 MHz, a DS2003 measurement acquisition card with 32 analog inputs, a DS2101 display card with five analog outputs, a DS3002 speed card with six high resolution inputs for incremental encoders, a DS4003 card with 96 logic inputs/outputs, and finally a DS5101 card with 16 PWM outputs of a resolution 25 ns [41]. The parameters of the dSPACE real-time digital simulator in our work have been configured as follows: the discretization time step Ts = 0.1 ms, the analogue outputs of ±10 V, and the computation method ode1 (Euler).

More details of the controls development and implementation with the dSPACE software under Matlab/Simulink environment can be found in [39–42].

#### *3.2. Real-Time Physical Emulator*

Real-Time Physical Emulator (RTPE) is the core module of the PHIL benchmark, which includes two key elements: the RTSS and the actuator (see Figure 10). The RTSS can digitally emulate the hydraulic behaviours by using the obtained mathematical model. In addition, the computed torque and the measured rotational speed information exchange at the interface. The RTSS computes the torque reference and sends it to the DCM control system, and the coupling turbine speed value is then sent back to the RTSS. The hydraulic model is generated in the Matlab/Simulink environment with the dSPACE modular (DS1005). The torque control is implemented in the DSP TMS320F240. The DC/DC power electronics converter is composed of a four-quadrant Pulse Width Module (PWM) IGBT chopper connected with a diode rectifier [19].

The RTPE could correctly emulate the hydraulic physical dynamics and the control algorithms could be executed in real time, and the following conditions must be met:


#### *3.3. Control Laws under Test*

The controllers under test are running with the Matlab/Simulink based on the dSPACE DS1005. The control input is then converted to switch commands for the power electronics converters in real time [19]. The digital variables and the analog measurements are exchanged via the I/O interface.

The control diagram of the grid-connected VS-HEP is presented in Figure 11. In the figure, *ig*(*a*, *b*, *c*) represents the grid currents, *im*(*a*, *b*, *c*) are the machine currents, *vg*(*a*, *b*, *c*) are the grid voltages, *i* ∗ *gd* and *i* ∗ *gq* are the *d*-axis and the *q*-axis grid currents references, respectively, *i* ∗ *md* and *i* ∗ *mq* are the *d*-axis and the *q*-axis machine currents references, respectively, *Pg* is the grid-injected active power, *ω* is the rotational speed, *ω*<sup>∗</sup> is the rotational speed reference, *Lf* and *Rf* are the filter inductance and the filter resistance, respectively, *Vdc* is the DC-link voltage, and *C* is the DC-link capacitance.

**Figure 11.** Control design diagram of grid-connected VS-HEP.

The VS-HEPs have two mostly common operating configurations: the grid following mode and the grid forming mode [43]. This paper is targeting on the grid following case. Then, the control design is mainly involved in two aspects as follows:


In fact, the PI controller can ensure a unit gain and zero phase lag for DC components or very low-varying signals [45]. Otherwise, if the reference signal has AC components, it is almost impossible to ensure the zero steady-state error in terms of its amplitude and phase. Moreover, the error would increase with the increasing frequency. The steady-state error can be reduced by increasing the controller gain. However, some limitations have been proved in previous works [45,46]. By introducing a high control gain, the closed-loop bandwidth would be increased. Finally, the system stability is possible to be deteriorated [46].

A PR controller (16) has a proportional term and a resonant integrating term with the frequencies to be corrected. The resonant term introduces theoretically infinite gains at the resonance frequencies, and the static error is hence eliminated:

$$H\_{PR}(\mathbf{s}) = k\_{pr} + \frac{2k\_r \mathbf{s}}{s^2 + \omega\_0^2} \tag{16}$$

where *ω*<sup>0</sup> is the resonant frequency, *kr* is the integration gain, and *kpr* is the proportional term. To reduce the noise sensitivity [46], the PR controller applied in this work is modified by integrating the parameter *ω<sup>c</sup>* as follows:

$$H\_{PR}(s) = k\_{p\tau} + \frac{2k\_r s}{s^2 + 2\omega\omega\_c s + \omega\_0^2} \tag{17}$$

The controller achieves infinite gains at the resonance frequencies for expression (16). The design corresponding to the (17) can still achieve a relatively high gain at the defined resonant frequency. The steady-state output phase and magnitude error in the closed-loop system would be approximately zero [46]. Furthermore, the multi-resonant controller (18) can be used to correctly control the currents with harmonics injected [47]:

$$H\_{MR} = k\_{pr} + \sum\_{i=1}^{h} \frac{2 \cdot k\_{ir} \cdot s}{s^2 + (\omega\_i^2)} \tag{18}$$

with *ω<sup>i</sup>* being the angle frequency of *i*-order harmonic, and *kir* the corresponding integration gain.

Note that, due to the variable speed operation, the frequency of machine currents are not constant but varying with time. Consequently, it is difficult to control machine side currents by using a PR controller with a fixed resonant frequency *ω*<sup>0</sup> [45]. This is why a PI controller upon the *dq* frame is more appropriate for the machine side converter control.

Lastly, it is worth mentioning that this work is dedicated to the modelling of hydraulic turbines for its implementation on PHIL benchmark. We have not provided many results regarding the comparisons between PI and PR controllers. More comparative experiments between them therefore will be conducted based on the PHIL real-time benchmark discussed in this work.

#### *3.4. Physical System Configurations*

The emulated hydraulic shaft by DCM is coupled to a PMSG connecting with power grids through back-to-back power electronics converters as shown in Figure 10b. Some ancillary devices such as voltage, current, and position sensors are included as well [34]. The key parameters of the VS-HEP PHIL simulation system are given in Table 2.


**Table 2.** Parameters of the PHIL experiment benchmark.

#### **4. PHIL Performance Verification and Analysis**

Experiments are conducted to assess performance of the obtained reduced-scale model for its implementation on the built PHIL benchmark. The static hydraulic features are firstly verified to be correctly replicated by the actuator (the DCM). Then, the key parameters regarding dynamics and the dynamic response of the real-time physical emulator are verified as well. The control laws employed to machine-side and grid-side controllers are validated in the end.

#### *4.1. Replication of Static Features*

In order to obtain the static curves, the opening ratio of guide vane and the turbine speed are set manually. The opening ratio of guide vane is firstly set to a specific value. When the actuator reaches the steady state, the rotational speed is then adjusted. Figure 12 illustrates the emulated hydraulic features and the regression models obtained in Section 2: the flow rate in Figure 12a, the torque in Figure 12b, and the output hydraulic power in Figure 12c. Firstly, the results indicate that the hydraulic model emulated by the RTPE can correctly replicate the established regression models. Then, we can observe some features from Figure 12: the water flow rate *Q* will decrease when the rotation speed increases for a constant valve opening *γ*; the water flow rate and the output torque both will increase when we increase the opening ratio of guide vane *γ*; there always exists an optimal rotational speed that maximizes the gained hydraulic power for each given guide vane opening *γ*.

**Figure 12.** Static hydraulic features assessment, (**a**) water flow rate vs rotational speed, (**b**) hydraulic turbine torque vs rotational speed, and (**c**) output hydraulic power vs rotational speed.

#### *4.2. Dynamics Verification*

In the dSPACE environment, the execution time of each model part can be measured via atomic subsystems, which has been answered in [48]. In a real-time program, there is a variable named 'turnaroundTime' for each task, which can indicate how much time it takes to execute the complete task. The key dynamics parameters are provided in Table 3. Firstly, the execution time taken on the hydraulic model is fast enough compared to the actuator's dynamics and the sampling period. Then, the overall executional time of the complete simulation system (the computational time of hydraulic model, the executional time of employed control laws, the time on digital/analog and analog/digital conversion process, etc.) is shorter than its sampling period. The hydraulic dynamics depend on its capacity, adjusted range, and its device mechanism, which can range from several seconds to minutes [49,50]. The operation time of the guide vane actuator in this work is around 2∼4 s. The torque control of DCM with a settling time of around 4*ms* ensures that the emulated mechanical torque can efficiently replicate the hydraulic dynamics.


**Table 3.** Key parameters regarding dynamics.

The hydraulic dynamic features in the offline simulation in Matlab/Simulink and the PHIL real-time simulation are presented in Figures 13 and 14, respectively. A soft starting is employed in the beginning; the variable speed control starts at *t* = 40 s, with the guide vane opening of *γ* = [0.5, 0.6, 0.5]. The two figures show that the emulated hydraulic torque *Thdy*, the water flow rate *Q*, and the output hydraulic power *Phdy* in real time could achieve similar dynamic features as that of offline simulations. However, we notice some differences on the dynamic behaviour; the most relevant difference is dealing with the overshoot that affects the torque. Two possible reasons could have caused this difference. Firstly, our PHIL implementation is based on the existed flexible hybrid PHIL benchmark in the laboratory. The reason could be coming from the actuator which emulates the turbine, i.e., the DC motor of Figure 10, as its internal control is not accessible, we cannot change its internal controllers. Another reason could be the aging of the physical system. The offline simulation models (PMSG, converters, line filter, etc.) are defined by applying the parameters originally provided by the manufacturer. However, the hardware system under test has been used for many years. The parameters of the PMSG, the values of line inductance and resistance, etc. could have changed due to the aging factor. The parameters shifts could cause the differences on the control dynamics response. In fact, this is also one of the critical reasons why the control laws must be validated with the real-time PHIL benchmark. Moreover, compared to offline simulations, one of the simulated parts is replaced by the hardware devices in the PHIL system, and the real physical constraints are thus taken into account in the real-time simulation loop. This would help more precisely evaluate the whole electric performances.

**Figure 13.** Hydraulic dynamics process in the offline simulation.

**Figure 14.** Hydraulic dynamics process in the real-time simulation, *Phdy* the generated hydraulic power (100 w/div), *Thdy* the turbine torque (2 N·m/div), *<sup>Q</sup>* the water flow rate (0.05 m3/s/div).

#### *4.3. Control Laws Test*

Figure 15 presents the machine-side control performances under the built PHIL benchmark with/without introducing the dynamic model of guide vane opening. In the control design, a soft-starting process is considered in order to smoothly start the machine as follows:


Compared to Figure 15a, the actuator dynamics process could deteriorate both speed and current control performances as shown in Figure 15b. Thus, it makes sense to assess the effect of the induced dynamics on the electrical-side control performance.

The grid-side control results are presented in Figures 16 and 17. In the soft staring process, the DC-link voltage increases until reaching the steady state; no hydraulic power is sent to the grid until the guide vane actuator acts. In the Maximum Power Point Tracking (MPPT) process, the active power *Pg* and reactive power *Qg* (controlled to be zero) injected to the power grid change with the adjustment of guide vane opening; the DC-link voltage *Vdc* almost stays constant. Figure 17 shows the control performance of grid-injected currents. In the beginning, the injected sinusoidal current synchronizes with the grid voltage of phase A well. Then, the reactive current reference steps from 0 A to 2 A, and a phase shift occurs between the grid voltage and the injected current. The grid-injected currents are correctly controlled by using PR controllers.

(**a**) Dynamic guide actuator model non-included

(**b**) Dynamic guide actuator model included

**Figure 15.** Machine-side control performances with/without introducing guide vane opening dynamic model, the *q*-axis current *imq* (5 A/div), the *d*-axis current *imd* (5 A/div), the rotational speed *ω* (20 rad/s/div).

**Figure 16.** Grid-side control performance, the DC-link voltage *Vdc* (100 V/div), the active power *Pg* (200 W/div), the reactive power *Qg* (200 VarA/div).

**Figure 17.** Grid-injected current control, *iga* the grid current of phase A (5 A/div), *vga* the grid voltage of phase A (100 V/div), *igq* the reactive current (2 A/div), *igd* the active current (2 A/div).

#### **5. Conclusions**

In this paper, a reduced-scale hydraulic modelling approach is proposed for achieving efficient PHIL simulations. The static features are modelled based on the 'Hill Charts' measurements, which are analytically generated by using optimization routines. Furthermore, the nonlinear dynamics of the actuator model are introduced in order to evaluate the effects of the induced dynamics on the whole electric performance. Moreover, the reduced-scale models can be established for different laboratory conditions. The suggested modelling approach enabled positive features for its implementations on PHIL simulation systems: decent accuracy, light computational complexity, and high scalability.

A bulb turbine coupled to a PMSG with back-to-back PWM converters is chosen as an example for PHIL design and performance assessment. Some conclusions could be highlighted based on the experimental results: Firstly, the RTPE can correctly replicate the hydraulic static features. In addition, the computational time taken on the RTPE is rapid enough to ensure that the emulated turbine torque can efficiently track the hydraulic dynamics. The controls of variable speed operation and the grid-side integration are validated under the established PHIL testing benchmark in the end. The established benchmark can be used for advanced control techniques study of VS-HEPs. Furthermore, the flexible PHIL real-time test benchmark can be adapted to different hydraulic regimes such as pumped-storage hydropower plant for future research.

**Author Contributions:** B.G. wrote this article, designed the PHIL simulation system method, implemented the hardware platform, and conducted the experimental tests and validations. A.M. contributed to the mathematical model of hydraulic turbine. S.B. and M.A. supervised, reviewed, and checked the article manuscript. C.B. helped to confirm the descriptions of dSPACE and the physical system. J.P. helped to acquire the financial support for the project leading to this publication and revised the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is supported by PSPC Innov'hydro project, which brings together General Electric (GE) Renewable Energy, EDF (Électricité de France), Grenoble INP, and other key players in the hydroelectric sector, in France. The publication is financed by the School of Engineering, HES-SO, Valais, Switzerland.

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

#### **References**


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

© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

*Review*
