Zero-Dynamics Attack on Wind Turbines and Countermeasures Using Generalized Hold and Generalized Sampler

Abstract

Most wind turbines are monitored and controlled by supervisory control and data acquisition systems that involve remote communication through networks. Despite the flexibility and efficiency that network-based monitoring and control systems bring, these systems are often threatened by cyberattacks. Among the various kinds of cyberattacks, some exploit the system dynamics so that the attack cannot be detected by monitoring system output, the zero-dynamics attack is one of them. This paper confirms that the zero-dynamics attack is fatal to wind turbines and the attack can cause system breakdown. In order to protect the system, we present two defense strategies using a generalized hold and a generalized sampler. These methods have the advantage that the zeros can be placed so that the zero dynamics of the system become stable; as a consequence, the zero-dynamics attack is neutralized. The effects of the countermeasures are validated through numerical simulations and the comparative discussion between two methods is provided.

1. Introduction

Wind energy has been recognized as one of the major renewable energy sources for over two decades. As of 2019, the global wind power capacity reaches 60 GW, including onshore and offshore plants [1]. Nowadays, it is typical to install multiple turbines in the same area where wind energy is abundant, this is called a wind farm. To manage a wind farm consisting of several to hundreds of turbines efficiently, a network-based monitoring and control scheme is needed, such as supervisory control and data acquisition (SCADA) systems. By virtue of SCADA systems, it is possible to acquire the condition of the wind turbines installed over a wide area and manage them efficiently.

However, network-based control and monitoring systems, including SCADA systems, are frequently targeted by cyberattacks due to the presence of the network. Recently, several attempts of cyberattacks on network-controlled plants have been reported, including a German steel mill [2], Iranian and American nuclear facilities [3], Ukrainian power plant [4], among others. These cases imply that many systems using network-based control and monitoring systems, including the wind turbine systems, are not free from cyberattacks. In order to secure the networked control systems from possible cyberattacks, researches on cyberattacks and defense strategies have drawn attention, see, e.g., [5,6,7,8], surveys [9,10,11] and the references therein.

Cyberattacks on systems controlled through networks can be classified into two categories, model-free attacks and model-based attacks. The model-free attacks are basically accomplished in the form of overloading the system through superfluous requests or stealing information from network lines. A most well-known model-free attack is DoS (denial-of-service) [12]. An attacker running DoS may overload or possess the network and system resources by sending numerous requests to the target system, disturbing the normal operation of systems. DDoS (distributed denial-of-service) [13], which conducts DoS attacks in a distributed form, and PSDoS (power save denial-of-service) [14] methods that can conduct DoS attacks with relatively small attack resources, are variants of DoS. Eavesdropping [15] is another instance of model-free attacks, which is accomplished by occupying part of the network lines. The eavesdroppers collect information by intercepting data or signals that are sent over the network. A Replay attack [16] is performed by combining the above two attacks. An attacker running the replay attack steals data being transmitted and sends it to the system repeatedly, pretending that the requests are from the valid user.

The model-based attacks are another category of cyberattack that include covert false-data injection attacks [17], zero-dynamics attacks (ZDAs) [8], robust zero-dynamics attacks [18], and pole-dynamics attacks [19]. An attacker performing the covert false-data injection attack, targeting a system that has known dynamics, can make the remote monitoring system or controller recognize that the system is under normal operation, even though it is not. ZDA is another model-based attack that manipulates the input signal only with the knowledge of the system dynamics. The attacker of ZDA can make the internal states of the system become unbounded while making the impact of the attack hardly visible on the output. The pole-dynamics attack is an attack that a malicious signal is injected into the output of the system. One of the major features of the model-based attack is that it allows more sophisticated and clever attacks. The attack signals, generated by using knowledge of the target systems, lead the remote monitoring system to be mistaken as a normal operation, which implies that the attacks are stealthy.

In this paper, we are particularly interested in the zero-dynamics attack on wind turbine systems. ZDA is known as one of the most fatal cyberattacks and this is mainly because the attack exploits the system model and very hard to detect. This attack requires a high level of model knowledge and has the character of disruption resources in the attack space of [8]. Suppose that a dynamic system that is represented by a transfer function and is stabilized by a controller, e.g., Proportional–Integral–Differential (PID) type. If the zero dynamics (which corresponds to the zeros of the transfer function) is unstable, then there exists an input signal that drives some internal variable of the system to be unbounded while unnoticed by monitoring the system output. Exploiting this fact, one can construct an undetectable cyberattack by copying the zero dynamics of the system. If this attack is applied to the system, then, by stability, the internal variable will approach the attack signal while the change of the output can hardly be detected [8].

As one might have noticed, ZDA is ineffective if the system has stable zero dynamics, since the attack is generated from a copy of the zero dynamics and thus converges to zero, meaning that the attack signal diminishes. Unfortunately, this does not mean that this system is safe from ZDA. In fact, most of the modern control systems are operated by digital controllers and it is quite often the case that the sampled-data system of the original continuous-time system has unstable zero dynamics; if the system has a relative degree greater than two, then at least one zero is unstable [20]. This means that even if the continuous-time system has stable zero dynamics, it can happen under ZDA that the sampled output remains constant while its continuous counterpart or some internal states diverge. In the first part of this paper, we demonstrate that the sampled-data model of the wind turbine system has unstable zero dynamics and thus, is vulnerable to ZDA; a ZDA is constructed so that the generator angular velocity diverges but its sampled values remain almost constant. Thus, from the input/output data of the wind turbine collected by the SCADA system, the system appears to be operating normally until some variable of the system reaches its hardware limit.

Recognizing the lethality of ZDA, several strategies to protect the system have been developed [7,21,22,23,24]. Authors in [7] proposed defense strategies to defend against such an attack by modifying the system structures including the input gain matrix, output matrix, and the system matrix. For example, the input gain matrix can be modified by adding and removing actuators or by introducing a perturbation. In [21], a modulation-matrix-based detection method was proposed, which is another structure-modifying method. These methods focus on the fact that the modification of the system structure leads to a change in zero dynamics. If the attacker does not know about the change of the system and injects the attack signal generated based on the system before the modification, the effect of the attack can be detected. However, if the attacker knows the modified system’s information, a stealthy attack is still possible. Thus, the information on the modified system should be hidden.

Another defense strategy is to make all zeros of the system become stable. Since ZDA is effective to a system that has at least one unstable zero, the methods under this strategy focus on proactively blocking the threat of attacks. Provided that the zero dynamics is stable, the attack signal based on the dynamics converges to zero. In [25], authors proposed an attack neutralizing method by measuring the output several times during one sampling interval. Then, the lifted system with new measurements has no unstable zero.

Recently, a generalized hold (GH)-based zero-assignment method was introduced in [23], and a method with a generalized sampler (GS) was proposed in [24]. The GH is an interfacing device between discrete-time signals and continuous-time signals, generating a continuous-time signal based on the predetermined hold function (or weights). By using GH instead of zero-order hold (ZOH), the zeros of the plant can be located in the stable region [20]. However, the GH-based defense method may cause undesirable intersample behavior, since the signal generated by GH is typically uneven.

The strategy employing GS is a more recently introduced method that overcomes the shortcomings of intersample behavior of the GH-based defense method. The GS constructs a new output by taking a weighted average of multiple output measurements from one sampling interval [20,24], replacing the simple sampler. It is shown that the zeros can be assigned arbitrarily by using GS, hence, it can be used as a countermeasure against ZDA. Since the operation of GS does not affect the plant input, undesirable intersample behavior no longer occurs. The output of GS, however, may differ from that of a simple sampler and can be sensitive to the sensor noise.

Most recently, cyberattacks including ZDA and their countermeasures are developed for more general systems having multiple agents. In [26], the concept of ZDA is generalized to a cooperative attack on multiagent systems considering network switching and a defense strategy involving the design of a series of network switching and the Luenberger observer has been proposed. In addition, the detection of cyberattacks including the false-data injection attack and replay attack is studied and applied to DC microgrids in [27], where the Luenberger observer and the unknown input observer are used to estimate the states of an agent and its neighbors.

In this paper, we apply two security strategies based on GH and GS to wind turbine systems. Firstly, a GH is designed so that all the zeros of the sampled-data model of the wind turbine system are located inside the unit circle. It is shown through numerical simulations that the presence of the ZDA that has been designed using the unstable sampling zero is revealed shortly after the attack is injected. The effect of GH on the intersample behavior of the generator angular velocity is also discussed. Secondly, we shift the unstable sampling zero into the stable region by employing the GS, which reveals the presence of ZDA clearly. It is emphasized that the undesirable intersample behavior does not appear anymore. In addition to this advantage, the numerical simulations demonstrate that sensitivity to noise can be reduced substantially by properly choosing the zeros.

The rest of this paper is organized as follows. In Section 2, we briefly describe the wind turbine’s modeling. Section 3 addresses that the digitally controlled wind turbine system is vulnerable to ZDA. Section 4 briefly explains how to change the zeros of a system by using GH and GS. In Section 5, we apply the design of GH and GS against ZDA and demonstrate that ZDA can be neutralized via numerical simulations. Finally, Section 6 concludes the paper.

2. Dynamic Model of Wind Turbine

In this section, we briefly review the dynamics of a wind turbine. Key components of a wind turbine are the rotor blade, drivetrain, generator, and interface to the main grid. The rotor blade is modeled as a static function that produces the mechanical torque $T_a$ that is applied to the drivetrain. We model the drivetrain as a dynamic system with three state variables. The dynamics of the generator are not considered but it is assumed that it can apply generating torque $T_g$ to the drivetrain.

Based on the wind turbine model, we derive a linear model and its discrete-time approximation around an operating point. Through numerical simulations, it is seen that the closed-loop system converges to the steady state under constant wind speed.

2.1. Dynamic Model of Wind Turbine

The drivetrain can be modeled as a two-inertia system in which rotor blades and generator are combined in one axis through a gearbox. Since the moment of inertia of the rotor blade is very high, the influence of the moment of inertia on the gearbox can be neglected [28]. Therefore, the influence of the gearbox is reflected on the generator side, and the wind turbine is expressed as a two-inertia model. Assuming that the shaft is thin and long, a two-inertia system can be modeled as a system with torsion, connected by a spring and a damper. It is assumed that the rotor has inertia $J_r$ and the generator has inertia $J_g$. They are connected through a torsional spring with spring constant $K_{sh}$ and a torsional damper with damping constant $D_{sh}$, as illustrated in Figure 1. The aerodynamic torque applied to the drivetrain is denoted by $T_a$ and the generator torque is denoted by $T_g$. $T_{sh}$ stands for the torsional torque developed in the shaft. ${\omega }_r$ and ${\omega }_g$ are the angular velocities corresponding to rotor and generator, respectively. For other types of drivetrain models, see, e.g., [29,30,31,32,33].

The dynamics of the drivetrain are derived as [34]

$$\begin{array}{cc}\hfill T_a& =J_r{\dot{\omega }}_r+T_{sh}\hfill \\ \hfill T_{sh}& =J_g{\dot{\omega }}_g+T_g\hfill \\ \hfill T_{sh}& =Q_s+Q_d=K_{sh}({\theta }_r-{\theta }_g)+D_{sh}({\omega }_r-{\omega }_g).\hfill \end{array}$$

(1)

We define the state vector $x={[{\omega }_r,{\omega }_g,Q_s]}^{\top }$ and express the dynamics (1) in the state space as

$$\begin{array}{cc}\hfill \dot{x}& =A_{s}x+B_{s}u\hfill \\ \hfill y& =C_{s}x,\hfill \end{array}$$

(2)

where $u={[T_a,T_g]}^{\top }$ is the input; y is the system output; and $A_s$, $B_s$, $C_s$ are system matrices given by

$$\begin{array}{c}\hfill A_s=\left[\begin{array}{ccc}-D_{sh}/J_r& D_{sh}/J_r& -1/J_r\\ D_{sh}/J_g& -D_{sh}/J_g& 1/J_g\\ K_{sh}& -K_{sh}& 0\end{array}\right],B_s=\left[\begin{array}{cc}1/J_r& 0\\ 0& -1/J_g\\ 0& 0\end{array}\right],C_s=\left[\begin{array}{ccc}0& 1& 0\end{array}\right]\end{array}.$$

The turbine power $P_t$ of the wind turbine and aerodynamic torque $T_a$ applied to the wind turbine are given by

$$\begin{array}{cc}\hfill P_t& =\frac{1}{2}\rho \pi R^2{C}_p(\lambda ,\beta )V^3\hfill \\ \hfill T_a& =\frac{1}{2}\rho \pi R^3{C}_q(\lambda ,\beta )V^2,\hfill \end{array}$$

(3)

where $\rho$ is the air density, R is the radius of the blade, V is the wind speed, $C_p(\lambda ,\beta )$ is the power coefficient, and $C_q(\lambda ,\beta )$ is the torque coefficient. $C_p(\lambda ,\beta )$ and $C_q(\lambda ,\beta )$ depend on the tip-speed-ratio (TSR) $\lambda$ defined by $\lambda ={\omega }_{r}R/V$ and the pitch angle $\beta$ [31,32,35]; $C_p$ and $C_q$ are related as $C_q=C_p/\lambda$. We take a widely used model of $C_p$ given by

$$\begin{array}{c}\hfill C_p(\lambda ,\beta )=c_1\left(c_2\frac{1}{\mathrm{\Lambda }}-c_3\beta -c_4{\beta }^{c_5}-c_6\right)e^{-c_7\frac{1}{\mathrm{\Lambda }}},\end{array}$$

(4)

where the parameter $\frac{1}{\mathrm{\Lambda }}$ is defined as

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Lambda }}=\frac{1}{\lambda +0.08\beta }-\frac{0.035}{1+{\beta }^3},\end{array}$$

and $c_1,\cdots ,c_7$ are constants that depend on system parameters [34,36]. Since $C_p(\lambda ,\beta )$ is related to the wind-rotor aerodynamic characteristics [37], the numerical values of $c_1,\cdots ,c_7$ depend on the wind turbine under consideration. In this paper, the parameters of the T-100 of Argolabe. S.L. Engineering Company are used, where $C_{p,max}$ is 0.4728 and $\lambda$ is 6 [38]. The $C_p(\lambda ,\beta )$ curve is extracted from Matlab/Simulink simulations with zero pitch angle ($\beta =0^{\circ }$). Through simulations, the parameters of $c_1,\cdots ,c_7$ are chosen as $c_1=0.29$, $c_2=115$, $c_3=0.5$, $c_4=0$, $c_5=0$, $c_6=6$, and $c_7=13.1$.

2.2. Discrete-Time Linear Model

Most modern control systems are controlled by digital devices. One of the well-established procedures to develop controllers is based on discretization of the system model. In this subsection, we first linearize the wind turbine system (2) with $T_a$ given by (3) around a point of the rated power operation and then discretize under the assumption that ZOH at the actuator side and simple sampler (SS) at the sensor side are used. It is emphasized that this discrete linear model containing unstable zeros can be used to generate an undetectable cyberattack, which will be discussed in the later part of this paper.

It is observed that the aerodynamic torque $T_a$ given in (3) is a nonlinear function of ${\omega }_r$ and V, and it can be linearized around an operating point $({\overline{\omega }}_r,\overline{V})$ (the rated power point) [30,34] as

$$\begin{array}{cc}\hfill {\widehat{T}}_a& =k_{{\omega }_r}({\overline{\omega }}_r,\overline{V}){\widehat{\omega }}_r+k_V({\overline{\omega }}_r,\overline{V})\widehat{V},\hfill \end{array}$$

(5)

where ${\overline{\omega }}_r$ and $\overline{V}$ are the values of ${\omega }_r$ and V at the operating point, ${\widehat{\omega }}_r={\omega }_r-{\overline{\omega }}_r$, $\widehat{V}=V-\overline{V}$, and ${\widehat{T}}_a=T_a-{\overline{T}}_a$ (${\overline{T}}_a$ is the aerodynamic torque at the operating point). The gains $k_{{\omega }_r}$ and $k_V$ are given by

$$\begin{array}{cc}\hfill k_{{\omega }_r}({\overline{\omega }}_r,\overline{V})& ={\left.\frac{\partial T_a}{\partial {\omega }_r}\right|}_{({\overline{\omega }}_r,\overline{V})}={\left.\frac{1}{2}\rho \pi R^4\overline{V}\frac{\partial C_q}{\partial \lambda }\right|}_{({\overline{\omega }}_r,\overline{V})}\hfill \\ \hfill k_V({\overline{\omega }}_r,\overline{V})& ={\left.\frac{\partial T_a}{\partial V}\right|}_{({\overline{\omega }}_r,\overline{V})}={\left.\frac{1}{2}\rho \pi R^3\overline{V}\left(2C_q-\lambda \frac{\partial C_q}{\partial \lambda }\right)\right|}_{({\overline{\omega }}_r,\overline{V})}.\hfill \end{array}$$

Substituting (5) to the dynamics (2), a linearized state-space model is obtained as

$$\begin{array}{cc}\hfill \dot{\widehat{x}}& ={\widehat{A}}_s\widehat{x}+{\widehat{B}}_s\widehat{u}+{\widehat{B}}_V\widehat{V}\hfill \\ \hfill \widehat{y}& =C_s\widehat{x},\hfill \end{array}$$

(6)

where $\widehat{x}=[{\widehat{\omega }}_r,{\widehat{\omega }}_g,{\widehat{Q}}_s]$; $\widehat{u}={\widehat{T}}_g$; and the matrices ${\widehat{A}}_s$, ${\widehat{B}}_s$, ${\widehat{B}}_V$, and $C_s$ are given by

$$\begin{array}{cc}\hfill {\widehat{A}}_s& =\left[\begin{array}{ccc}(k_{{\omega }_r}-D_{sh})/J_r& D_{sh}/J_r& -1/J_r\\ D_{sh}/J_g& -D_{sh}/J_g& 1/J_g\\ K_{sh}& -K_{sh}& 0\end{array}\right],{\widehat{B}}_s=\left[\begin{array}{c}0\\ -1/J_g\\ 0\end{array}\right]\hfill \\ \hfill {\widehat{B}}_V& =\left[\begin{array}{c}{k}_V/J_r\\ 0\\ 0\end{array}\right],C_s=\left[\begin{array}{ccc}0& 1& 0\end{array}\right].\hfill \end{array}$$

From this linear model, the transfer function from ${\widehat{T}}_g$ to ${\widehat{\omega }}_g$ is computed as

$$\begin{array}{c}\hfill G\left(s\right)=\frac{{\widehat{\omega }}_g\left(s\right)}{{\widehat{T}}_g\left(s\right)}=-\frac{b_2{s}^2+b_{1}s+b_0}{a_3{s}^3+a_2{s}^2+a_{1}s+a_0},\end{array}$$

(7)

where

$$\begin{array}{cc}\hfill a_3& =J_r{J}_g,a_2=D_{sh}(J_r+J_g)-k_{{\omega }_r}({\overline{\omega }}_r,\overline{V})J_g\hfill \\ \hfill a_1& =K_{sh}(J_r+J_g)-k_{{\omega }_r}({\overline{\omega }}_r,\overline{V})D_{sh},a_0=-k_{{\omega }_r}({\overline{\omega }}_r,\overline{V})K_{sh}\hfill \\ \hfill b_2& =J_r,b_1=D_{sh}-k_{{\omega }_r}({\overline{\omega }}_r,\overline{V}),b_0=K_{sh}.\hfill \end{array}$$

From the linearized model (6), we derive a discrete-time model for the purpose of controller design. Let $T_s$ be the sampling time and suppose that the control input $u_k$ is determined by a digital controller and the actual control input $\widehat{u}\left(t\right)$ is generated via ZOH so that $\widehat{u}\left(t\right)={\widehat{u}}_k$ for $k{T}_s\le t<(k+1)T_s$, where k is a non-negative integer. The measured output that is transmitted to the controller is sampled at each sampling time $k{T}_s$, we call this sampling device a simple sampler, and this sampled output is denoted by ${\widehat{y}}_k$. Assuming the wind turbine is at steady state under constant wind speed, we have the following discrete-time linear model

$$\begin{array}{cc}\hfill {\widehat{x}}_{k+1}& ={\widehat{A}}_{\mathrm{d}}{\widehat{x}}_k+{\widehat{B}}_{\mathrm{d}}{\widehat{u}}_k\hfill \\ \hfill {\widehat{y}}_k& =C_{\mathrm{d}}{\widehat{x}}_k,\hfill \end{array}$$

(8)

where ${\widehat{x}}_k=\widehat{x}\left(k{T}_s\right)$ is the state vector, ${\widehat{A}}_{\mathrm{d}}=e^{{\widehat{A}}_s{T}_s}$, ${\widehat{B}}_{\mathrm{d}}={\int }_0^{T_s}{e}^{{\widehat{A}}_s(T_s-\tau )}{\widehat{B}}_{s}d\tau$, and $C_d=C_s$. It is noted that this model is an exact discretization of the continuous-time model (6).

We use a PI-type digital controller with sampling time $T_s$ whose discrete-time transfer is given by

$$\begin{array}{c}\hfill C_{\mathrm{PI}}\left(z\right)=K_P+K_I\frac{T_s}{z-1},\end{array}$$

where $K_P$ and $K_I$ are the proportional gain and integral gain, respectively. The discrete control input signal ${\widehat{u}}_k$ is computed from the relation ${\widehat{u}}_k={\mathcal{Z}}^{-1}(C_{\mathrm{PI}}\left(z\right)(0-\widehat{Y}(z))$ where $\widehat{Y}\left(z\right)$ is the z-transform of ${\widehat{y}}_k$.

2.3. Wind Turbine Simulation Model and Its Behavior under Normal Condition

We take a small size wind turbine whose system parameters, summarized in

Table 1. Model parameters of a small-sized wind turbine.

${\mathit{J}}_{\mathit{r}}$	${\mathit{J}}_{\mathit{g}}$	${\mathit{K}}_{\mathit{s}\mathit{h}}$	R	${\mathit{\lambda }}_{\mathit{o}\mathit{p}\mathit{t}}$	${\mathit{C}}_{\mathit{p},\mathit{m}\mathit{a}\mathit{x}}$
30,375 kg $·{\mathrm{m}}^2$	151 kg $·{\mathrm{m}}^2$	0.31$\times {10}^6$ N·m/rad	11.25 m	6	47.24%

, are taken from [33] and the data-sheet of the T-100 wind turbine of Argolabe S.L. Engineering Company [38]. It is a horizontal axis wind turbine and is designed for distributed generation and (or) electric self-consumption applications, connected to a power grid. According to the data-sheet, the shaft of the rotor and the generator are connected by a gearbox, but by including the gearbox in the inertia of the generator, it can be modeled as a two-inertia model connected by one shaft with torsion, as shown in Figure 1.

Figure 2 is the result of the turbine behavior simulation obtained using the simulation tool Matlab/Simulink. As shown in Figure 2a, wind speed is simulated with the numerical data described under the scenario that the wind speed changes from 0 to 10.5 m/s and reaches 10.5 m/s at 100 s. Figure 2b shows the generated power calculated using (3), it is observed that the generated power in the turbine changes according to the wind speed (V) and power coefficient ($C_p$) and that the generated power reaches its maximum after about 200 s. The maximum power is about 133 kW. Figure 2c shows the response of the power coefficient ($C_p$), which reaches 0.4724 after about 200 s. The generator angular velocity (${\omega }_g$) and rotor angular velocity (${\omega }_r$) are plotted in Figure 2d. According to the specifications [38], the rated speed of the rotor is 5.6 rad/s and the gear ratio is 22.2. Using this information, the generator speed ${\omega }_g$ is plotted considering the gear ratio, of which the rated value is 124.32 rad/s. A PI controller is used to control the angular velocity of the generator and the gains are chosen as $K_p$ = 400 and $K_I$ = 600. In Figure 2d, as the wind speed changes—as shown in Figure 2a—it is observed that the generator angular velocity reaches 124.32 rad/s after 200 s so that the rotor angular velocity reaches its rated value.

It is noted that the steady state behavior of this wind turbine is obtained under the condition that no cyberattack has been injected. In what follows, we consider the case where the system is under a ZDA that is obtained from this steady state condition (i.e., the information on the operating point) and show that the presence of the attack cannot be detected for a considerably long time interval.

3. Zero-Dynamics Attack on Wind Turbine

Thanks to the development of communication and computing technology, modern wind turbine systems are generally operated by remote control and monitoring systems such as SCADA systems [39,40,41]. Operating wind power generators using the SCADA system enables efficient monitoring and control over wide areas such as wind farms [42,43]. An operator of the SCADA system can comprehensively manage the system based on the signals delivered through the network such as rotor and generator speed, wind speed, and generator torque; see, e.g., [39,40] for details. For the reasons discussed above, the network-based monitoring (or controlling) systems bring an advantage in terms of efficient management of the system. Through a communication network, multiple wind turbines installed in a wide area can be managed in one place, which can result in reduced manpower and immediate state monitoring [44,45].

However, there are also problems with network usage, such as vulnerabilities to cyberattack. If a large area is connected through communication lines, the number of paths through which an attacker with a malicious purpose can inject signals into the system may increase [46,47,48]. In fact, many studies on cyberattacks using these vulnerabilities have been conducted recently [8,49], such as data integrity attack [6], false-data injection attack [50], ZDA [7], etc. Among these cyberattacks, ZDA is one of the sophisticated cyberattacks based on system dynamics [23,24].

Consider a wind turbine system that operates as shown in Figure 3. The turbine receives control input from a remote controller and feeds the system state back to the controller through a network. ZDA becomes possible by occupying the network line between the controller and target system.

Assume that the attacker who has taken control of the network line knows the model information of the system and can add the attack signal $a_k$ to the controller output, as shown in Figure 4. Then, using the system model, the attacker can generate a sophisticated attack signal that enables the discrete-time output signal to pretend that the system is operating normally but make the internal state actually unbounded [51,52]. The mathematical explanation for the zero-dynamics attack is as follows:

From (8), the dynamics of the wind turbine under ZDA becomes

$$\begin{array}{cc}\hfill {\widehat{x}}_{k+1}& ={\widehat{A}}_{\mathrm{d}}{\widehat{x}}_k+{\widehat{B}}_{\mathrm{d}}({\widehat{u}}_k+a_k)\hfill \\ \hfill {\widehat{y}}_k& =C_{\mathrm{d}}{\widehat{x}}_k,\hfill \end{array}$$

(9)

where $a_k$ is the attack signal of ZDA. The attack signal $a_k$ is generated from a dynamic system which is identical to the zero dynamics of the system. We recall that the zero dynamics of the system (9) can be identified by rewriting it in the normal form [53] given by

$$\begin{array}{cc}\hfill {\eta }_{k+1}& =S_{\mathrm{d}}{\eta }_k+P_{\mathrm{d}}{\xi }_k,{\eta }_k\in {\mathbb{R}}^2\hfill \\ \hfill {\xi }_{k+1}& ={\psi }_{\mathrm{d}}^{\top }{\eta }_k+{\varphi }_{\mathrm{d}}{\xi }_k+g_{\mathrm{d}}({\widehat{u}}_k+a_k),{\xi }_k\in \mathbb{R}\hfill \\ \hfill {\widehat{y}}_k& ={\xi }_k,\hfill \end{array}$$

(10)

where the dynamics of ${\xi }_k$ explains how the input $u_k$ directly affects the system output (${\widehat{y}}_{k+1}$ explicitly depends on $u_k$) and that of ${\eta }_k$ describes the internal behavior of the system. The dynamics ${\eta }_{k+1}=S_{\mathrm{d}}{\eta }_k$ is called the zero dynamics and the eigenvalues of $S_{\mathrm{d}}$ correspond to the zeros of the system (9).

The attack signal $a_k$ of ZDA is generated from a dynamic system given by

$$a_k=-\frac{1}{g_{\mathrm{d}}}{\psi }_{\mathrm{d}}^{\top }{z}_k,z_{k+1}=S_{\mathrm{d}}{z}_k,$$

(11)

where $z_k\in {\mathbb{R}}^2$ is the state of attack generator. It is noted that the attack is constructed using the system parameters such as $S_{\mathrm{d}}$, ${\psi }_{\mathrm{d}}$, and $g_{\mathrm{d}}$.

Now, we investigate the behavior of the closed-loop system under ZDA. Firstly, the controller $C_{\mathrm{PI}}\left(z\right)$ represented in the state space

$$\begin{array}{cc}\hfill c_{k+1}& =c_k+e_k,e_k=-{\widehat{y}}_k\hfill \\ \hfill {\widehat{u}}_k& =k_I{c}_k+k_P{e}_k,\hfill \end{array}$$

has been designed so that the closed-loop system under no ZDA is stable, as demonstrated in Section 2.3, namely, the matrix ${\widehat{A}}_{\mathrm{CL}}$ shown below is Schur.

$${\widehat{A}}_{\mathrm{CL}}=\left[\begin{array}{ccc}{S}_{\mathrm{d}}& P_{\mathrm{d}}& 0\\ {\psi }_{\mathrm{d}}^{\top }& {\varphi }_{\mathrm{d}}-g_{\mathrm{d}}{k}_P& g_{\mathrm{d}}{k}_I\\ 0& -1& 1\end{array}\right].$$

(12)

When the attack $a_k$ generated by (11) is injected into system (10) (it is equivalent to system (9)), we can derive

$$\left[\begin{array}{c}{\eta }_{k+1}\\ {\xi }_{k+1}\\ c_{k+1}\\ z_{k+1}\end{array}\right]=\left[\begin{array}{cccc}{S}_{\mathrm{d}}& P_{\mathrm{d}}& 0& 0\\ {\psi }_{\mathrm{d}}^{\top }& {\varphi }_{\mathrm{d}}-g_{\mathrm{d}}{k}_P& g_{\mathrm{d}}{k}_I& -{\psi }_{\mathrm{d}}^{\top }\\ 0& -1& 1& 0\\ 0& 0& 0& S_{\mathrm{d}}\end{array}\right]\left[\begin{array}{c}{\eta }_k\\ {\xi }_k\\ c_k\\ z_k\end{array}\right],$$

from which one has

$$\left[\begin{array}{c}{\eta }_{k+1}-z_{k+1}\\ {\xi }_{k+1}\\ c_{k+1}\end{array}\right]={\widehat{A}}_{\mathrm{CL}}\left[\begin{array}{c}{\eta }_k-z_k\\ {\xi }_k\\ c_k\end{array}\right].$$

Since ${\widehat{A}}_{\mathrm{CL}}$ given in (12) is Schur, there exist $\kappa >0$ and $\lambda$ with $\left|\lambda \right|<1$ such that

$$∥\left[\begin{array}{c}{\eta }_k-z_k\\ {\xi }_k\\ c_k\end{array}\right]∥\le \kappa {\lambda }^k∥\left[\begin{array}{c}{\eta }_0-z_0\\ {\xi }_0\\ c_0\end{array}\right]∥.$$

(13)

This relation implies that under a ZDA, the internal state ${\eta }_k$ converges to the state of ZDA, while other states ${\xi }_k$ and $c_k$ converge to zero.

It is remarkable that the relation (13) holds regardless of the stability of the zero dynamics. Hence, if $S_{\mathrm{d}}$ is unstable, then the internal state ${\eta }_k$ diverges whenever $z_0$ depends on an unstable eigenvector of $S_{\mathrm{d}}$, while this behavior cannot be observed by monitoring the signal $({\xi }_k,c_k)$.

Unfortunately, the wind turbine system (8) has an unstable zero when the sample time $T_s$ belongs to some region. Figure 5a shows the locus of the zeros with respect to the sampling time from $T_s=0.001$ s to $T_s=0.1$ s. When $T_s=0.001$ s, two zeros are located near 1, as depicted by the crosses. When $T_s$ increases and then becomes $0.08$ s, $z_2$ is located outside the unit circle. The blue and red circles on the real axis in Figure 5a indicate the location of two zeros when $T_s=0.1$ s. As shown in Figure 5b, if $T_s\in [0.08,0.13]$ s, at least one zero is located outside the unit circle.

Suppose that $T_s$ belongs to the region where the wind turbine has an unstable zero. When the attack (11) is injected, the rotor angular velocity ${\omega }_r$ and the (spring) torsional torque $Q_s$ will diverge while the generator angular velocity ${\omega }_g$ and the controller state $c_k$ converge to zero. This can be interpreted so that the attack intentionally moves two components (${\omega }_r$ and $Q_s$) of the operating point but leaves ${\omega }_g$ at the normal operating point and deceives the controller as if all components remain unchanged.

Under the situation described in Figure 4, the effect of ZDA on the wind turbine system is presented. Since the system is of nonminimum phase when $T_s=0.1$ s, the attack signal diverges for an appropriately chosen initial condition as shown in Figure 6, and it is expected that the internal states also diverge. Suppose that the hacker injected $a_k$ through the communication network at $t=0$. Let ${\widehat{\omega }}_{g,\mathrm{Th}}$ be a threshold that the steady state value of ${\widehat{\omega }}_g$ should not exceed (the dotted line in Figure 7a and ${\widehat{\omega }}_{g,\mathrm{Th}}=0.18$ rad/s), in other words, if $|{\widehat{\omega }}_g|>{\widehat{\omega }}_{g,\mathrm{Th}}$, the monitoring system determines that a fault has occurred or an attack has been injected. As can be seen in Figure 7a, the continuous-time output $y\left(t\right)={\omega }_g\left(t\right)$—the generator speed—becomes unbounded, while the discrete-time signal $y_k={\omega }_{g,\mathrm{ZOH}}$ that is transmitted to the controller and the SCADA system remain almost unchanged, indicating that the wind turbine still operates normally.

Meanwhile, the attack also affects the generated power and internal state. Figure 7b,c show the responses of the internal states ${\omega }_r$ and $Q_s$, respectively, and it is observed that ${\omega }_r$ and $Q_s$ become unbounded. In addition, the generated power under ZDA decreases, as shown in Figure 7d.

4. Two Countermeasures against Zero-Dynamics Attack

ZDA becomes effective when the discrete-time system has an unstable zero and, as discussed in Section 3, it can happen even if the continuous-time system has stable zero dynamics. Among several countermeasures to ZDA, we introduce two strategies that share the same idea, shifting zeros. These approaches are based on the fact that the zeros of the discrete-time system can be arbitrarily assigned if the ZOH is replaced by a generalized hold (GH) or if a generalized sampler (GS) is used instead of SS [20]. Applications of these ideas to security problems are reported in [23,24], and in this section, we apply these approaches to wind turbine systems.

4.1. Generalized-Hold-Based Strategy

GH has been introduced in [20] and involves a function $h_g\left(t\right)$, the so-called hold function that is defined as a piecewise continuous function $h_g$ so that the actual input applied to the system is given by

$$\begin{array}{c}\hfill \widehat{u}\left(t\right)=\sum _{k=-\infty }^{\infty }{h}_g(t-k{T}_s){\widehat{u}}_k.\end{array}$$

If a GH having a hold function $h_g\left(t\right)$ is used instead of ZOH, the sampled-data model of wind turbine system (6) under constant wind speed $\overline{V}$ (so that $\widehat{V}=0$) and ZDA becomes

$$\begin{array}{cc}\hfill {\widehat{x}}_{k+1}& ={\widehat{A}}_{\mathrm{d}}{\widehat{x}}_k+{\widehat{B}}_g({\widehat{u}}_k+a_k)\hfill \\ \hfill {\widehat{y}}_k& =C_{\mathrm{d}}{\widehat{x}}_k,\hfill \end{array}$$

(14)

where ${\widehat{A}}_{\mathrm{d}}=e^{{\widehat{A}}_s{T}_s}$, ${\widehat{B}}_g={\int }_0^{T_s}{e}^{{\widehat{A}}_s(T_s-\tau )}{\widehat{B}}_s{h}_g\left(\tau \right)d\tau$, and $C_{\mathrm{d}}=C_s$. The discrete-time transfer function from the generator torque ${\widehat{u}}_k(={\widehat{T}}_g)$ (equivalently from the attack $a_k$) to generator angular velocity ${\widehat{y}}_k(={\widehat{\omega }}_g)$, denoted by $G_{\mathrm{d}}\left(z\right)$, is then given by

$$\begin{array}{c}\hfill G_{\mathrm{d}}\left(z\right)=C_{\mathrm{d}}{(zI-{\widehat{A}}_{\mathrm{d}})}^{-1}{\widehat{B}}_g.\end{array}$$

(15)

It is emphasized that since $({\widehat{A}}_s,{\widehat{B}}_s)$ is controllable, one can always find a hold function $h_g$ so that the zeros of $G_{\mathrm{d}}\left(z\right)$ can be placed anywhere in the complex plane [20,23].

Let $z_{d,1}$ and $z_{d,2}$ be the desired zeros located inside the unit circle and $k_d$ be a gain. The problem is to find ${\widehat{B}}_g$ such that the transfer function $G_{\mathrm{d}}\left(z\right)$ has desired zeros and gain, i.e., the following identity holds

$$\begin{array}{c}\hfill G_{\mathrm{d}}\left(z\right)=C_{\mathrm{d}}{(zI-{\widehat{A}}_{\mathrm{d}})}^{-1}{\widehat{B}}_g=k_d\frac{(z-z_{d,1})(z-z_{d,2})}{\det (zI-{\widehat{A}}_{\mathrm{d}})}=:G_{\mathrm{d}}^{*}\left(z\right).\end{array}$$

Let $\det (zI-{\widehat{A}}_{\mathrm{d}})=z^3+d_2{z}^2+d_{1}z+d_0$. Then, $G_{\mathrm{d}}^{*}\left(z\right)$ can be realized in the control canonical form [54] given by

$$\begin{array}{cc}\hfill {\overline{x}}_{k+1}& =A_{\mathrm{con}}{\overline{x}}_k+B_{\mathrm{con}}{\overline{u}}_k\hfill \\ \hfill {\overline{y}}_k& =C_{\mathrm{con}}{\overline{x}}_k,\hfill \end{array}$$

(16)

where

$$\begin{array}{cc}\hfill A_{\mathrm{con}}& =\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -d_0& -d_1& -d_2\end{array}\right],B_{\mathrm{con}}=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right],C_{\mathrm{con}}=\left[\begin{array}{ccc}{k}_d{z}_{d,1}{z}_{d,2}& -k_d(z_{d,1}+z_{d,2})& k_d\end{array}\right].\hfill \end{array}$$

Equating the Markov parameters of the two transfer functions, one has

$$\begin{array}{cc}\hfill C_{\mathrm{d}}{\widehat{A}}_{\mathrm{d}}^k{\widehat{B}}_g& =C_{\mathrm{con}}{A}_{\mathrm{con}}^k{B}_{\mathrm{con}},k=0,1,\cdots ,\hfill \end{array}$$

from which ${\widehat{B}}_g$ is determined by

$${\widehat{B}}_g={\left[\begin{array}{c}{C}_{\mathrm{d}}\\ C_{\mathrm{d}}{\widehat{A}}_{\mathrm{d}}\\ C_{\mathrm{d}}{\widehat{A}}_{\mathrm{d}}^2\end{array}\right]}^{-1}\left[\begin{array}{c}{C}_{\mathrm{con}}\\ C_{\mathrm{con}}{A}_{\mathrm{con}}\\ C_{\mathrm{con}}{A}_{\mathrm{con}}^2\end{array}\right]B_{\mathrm{con}},$$

(17)

where the invertibility is assured by the observability of $({\widehat{A}}_{\mathrm{d}},C_{\mathrm{d}})$.

As discussed in [20,23], one candidate of GH is to use a piecewise constant function given by

$$\begin{array}{c}\hfill h_g\left(t\right)=h_i,\frac{(i-1)T_s}{N}\le t<\frac{i{T}_s}{N},i=1,\cdots ,N,\end{array}$$

where $h_i$ are constant gains and N is the number of subintervals. It can be shown that the gains $h_i$ and the vector ${\widehat{B}}_g$ are related as

$${\widehat{B}}_g=\sum _{l=1}^N{h}_l{\int }_{\frac{(l-1)T_s}{N}}^{\frac{l{T}_s}{N}}{e}^{{\widehat{A}}_s(T_s-\tau )}{\widehat{B}}_{s}d\tau$$

(18)

and this can be rewritten as

$${\widehat{B}}_g=\left[\begin{array}{cccc}{A}_{\mathrm{d},N}^{N-1}{B}_{\mathrm{d},N}& \cdots & A_{\mathrm{d},N}{B}_{\mathrm{d},N}& B_{\mathrm{d},N}\end{array}\right]h=:{\mathcal{C}}_{\mathrm{d},N}h,$$

where $h={[h_1,\cdots ,h_N]}^{\top }$ and

$$\begin{array}{c}\hfill A_{\mathrm{d},N}=e^{{\widehat{A}}_s\frac{T_s}{N}},B_{\mathrm{d},N}={\int }_0^{\frac{T_s}{N}}{e}^{{\widehat{A}}_s(\frac{T_s}{N}-\tau )}{\widehat{B}}_{s}d\tau .\end{array}$$

The hold gains are then computed as

$$\begin{array}{c}\hfill h=C_{d,N}^{†}{\widehat{B}}_g.\end{array}$$

(19)

For more details on the derivation, see [23].

4.2. Generalized-Sampler-Based Approach

By GS, we mean a device that generates a discrete-time signal ${\breve{y}}_k$ from a continuous-time signal $y\left(t\right)$ in a way that N measurements (i.e., generator angular velocity) $y(\frac{1}{N}{T}_s+(k-1)T_s),y(\frac{2}{N}{T}_s+(k-1)T_s),\cdots ,y\left(k{T}_s\right)$ are taken from the sampling interval $((k-1)T_s,k{T}_s]$ and a weighted average of them is computed as

$$\begin{array}{c}\hfill {\breve{y}}_k=\sum _{i=1}^N{w}_{i}y\left(\frac{i}{N}{T}_s+(k-1)T_s\right),\end{array}$$

(20)

where $w_1,\cdots ,w_N$ are weights for GS.

Similar to the case of GH, we can rewrite the system (6) under constant wind speed and ZDA as

$$\begin{array}{cc}\hfill {\widehat{x}}_k& ={\widehat{A}}_{\mathrm{d}}{\widehat{x}}_{k-1}+{\widehat{B}}_{\mathrm{d}}({\widehat{u}}_{k-1}+a_{k-1})\hfill \\ \hfill {\breve{y}}_k& ={\breve{C}}_{\mathrm{d}}{\widehat{x}}_{k-1}+{\breve{D}}_{\mathrm{d}}({\widehat{u}}_{k-1}+a_{k-1}),\hfill \end{array}$$

(21)

where

$$\begin{array}{cccc}\hfill {\widehat{A}}_{\mathrm{d}}& =e^{{\widehat{A}}_s{T}_s},\hfill & \hfill {\widehat{B}}_{\mathrm{d}}& ={\int }_0^{T_s}{e}^{{\widehat{A}}_s(T_s-\tau )}{\widehat{B}}_{s}d\tau \hfill \\ \hfill {\breve{C}}_{\mathrm{d}}& =\sum _{i=1}^N{w}_i{C}_{\mathrm{d}}{e}^{{\widehat{A}}_s\frac{i}{N}{T}_s},\hfill & \hfill {\breve{D}}_{\mathrm{d}}& =\sum _{i=1}^N{w}_i{C}_{\mathrm{d}}{\int }_0^{\frac{i}{N}{T}_s}{e}^{{\widehat{A}}_s(\frac{i}{N}{T}_s-\tau )}{\widehat{B}}_{s}d\tau .\hfill \end{array}$$

From (21), we can compute the transfer function from ${\widehat{u}}_k$ to ${\widehat{y}}_k$ as

$$\begin{array}{c}\hfill G_{\mathrm{d}}(z)=z^{-1}({\breve{C}}_{\mathrm{d}}{(zI-{\widehat{A}}_{\mathrm{d}})}^{-1}{\widehat{B}}_{\mathrm{d}}+{\breve{D}}_{\mathrm{d}}).\end{array}$$

(22)

Note that ${\breve{C}}_{\mathrm{d}}$ and ${\breve{D}}_{\mathrm{d}}$ contain the sampler weights $w_1$, ⋯, $w_N$ of GS, which are design parameters. If the weights are chosen appropriately, it is expected that the numerator of the transfer function (22) can be chosen as desired. In fact, this is true under mild assumptions [24].

Let $z_{d,1},z_{d,2}$, and $z_{d,3}$ be the desired zeros whose magnitudes are less than 1. We want to find the weights of GS such that the transfer function $G_{\mathrm{d}}\left(z\right)$ becomes identical to

$$\begin{array}{c}\hfill G_{\mathrm{d}}^{*}\left(z\right)=k_d{z}^{-1}\frac{(z-z_{d,1})(z-z_{d,2})(z-z_{d,3})}{\det (zI-{\widehat{A}}_{\mathrm{d}})},\end{array}$$

(23)

where $k_d$ is a high-frequency gain. To proceed, let $c_0$, $c_1$, $c_2$ be such that $(z-z_{d,1})(z-z_{d,2})(z-z_{d,3})=z^3+c_2{z}^2+c_{1}z+c_0$. We first find ${\breve{C}}_{\mathrm{d}}$ and ${\breve{D}}_{\mathrm{d}}$, then determine the weights $w_i$. As in the case of GH, we realize (23) in the control canonical form given by

$$\begin{array}{cc}\hfill {\overline{x}}_k& =A_{\mathrm{con}}{\overline{x}}_{k-1}+B_{\mathrm{con}}{\overline{u}}_{k-1}\hfill \\ \hfill {\overline{y}}_k& =C_{\mathrm{con}}{\overline{x}}_{k-1}+D_{\mathrm{con}}{\overline{u}}_{k-1},\hfill \end{array}$$

where $A_{\mathrm{con}}$ and $B_{\mathrm{con}}$ are identical to those of (16), and

$$\begin{array}{c}\hfill C_{\mathrm{con}}=\left[\begin{array}{ccc}{k}_d(c_0-d_0)& k_d(c_1-d_1)& k_d(c_2-d_2)\end{array}\right],D_{\mathrm{con}}=k_d.\end{array}$$

From the fact that two transfer functions $G_{\mathrm{d}}\left(z\right)$ and $G_{\mathrm{d}}^{*}\left(z\right)$ are identical if and only if

$$\begin{array}{cc}\hfill {\breve{D}}_{\mathrm{d}}& =D_{\mathrm{con}}\hfill \\ \hfill {\breve{C}}_{\mathrm{d}}{\widehat{A}}_{\mathrm{d}}^k{\widehat{B}}_{\mathrm{d}}& =C_{\mathrm{con}}{A}_{\mathrm{con}}^k{B}_{\mathrm{con}},k=0,1,\cdots ,\hfill \end{array}$$

we have, from the controllability of $({\widehat{A}}_{\mathrm{d}},{\widehat{B}}_{\mathrm{d}})$,

$$\begin{array}{cc}\hfill {\breve{C}}_{\mathrm{d}}& =C_{\mathrm{con}}\left[\begin{array}{ccc}{B}_{\mathrm{con}}& A_{\mathrm{con}}{B}_{\mathrm{con}}& A_{\mathrm{con}}^2{B}_{\mathrm{con}}\end{array}\right]{\left[\begin{array}{ccc}{\widehat{B}}_{\mathrm{d}}& {\widehat{A}}_{\mathrm{d}}{\widehat{B}}_{\mathrm{d}}& {\widehat{A}}_{\mathrm{d}}^2{\widehat{B}}_{\mathrm{d}}\end{array}\right]}^{-1}\hfill \\ \hfill {\breve{D}}_{\mathrm{d}}& =k_d.\hfill \end{array}$$

With ${\breve{C}}_{\mathrm{d}}$ and ${\breve{D}}_{\mathrm{d}}$ obtained above, it follows from the relation between the weights and (${\breve{C}}_{\mathrm{d}}$, ${\breve{D}}_{\mathrm{d}}$) that

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{\breve{C}}_{\mathrm{d}}& {\breve{D}}_{\mathrm{d}}\end{array}\right]=w\left[\begin{array}{cc}{C}_{\mathrm{d}}{e}^{{\widehat{A}}_s\frac{1}{N}{T}_s}& C_{\mathrm{d}}{\int }_0^{\frac{1}{N}{T}_s}{e}^{{\widehat{A}}_s(\frac{1}{N}{T}_s-\tau )}{\widehat{B}}_{s}d\tau \\ C_{\mathrm{d}}{e}^{{\widehat{A}}_s\frac{2}{N}{T}_s}& C_{\mathrm{d}}{\int }_0^{\frac{2}{N}{T}_s}{e}^{{\widehat{A}}_s(\frac{2}{N}{T}_s-\tau )}{\widehat{B}}_{s}d\tau \\ ⋮& ⋮\\ C_{\mathrm{d}}{e}^{{\widehat{A}}_s{T}_s}& C_{\mathrm{d}}{\int }_0^{T_s}{e}^{{\widehat{A}}_s(T_s-\tau )}{\widehat{B}}_{s}d\tau \end{array}\right]=:wM\end{array},$$

and the weights are computed as

$$\begin{array}{c}\hfill w=\left[\begin{array}{cc}{\breve{C}}_{\mathrm{d}}& {\breve{D}}_{\mathrm{d}}\end{array}\right]M^{†},\end{array}$$

(24)

where $M^{†}$ is the pseudo-inverse of M. For more details, see [24].

5. Evaluation of Countermeasures against ZDA

In this section, we apply the theory given in Section 4 and demonstrate that the two countermeasures that can shift the zeros into the stable region effectively reveal the presence of ZDA. Although both approaches work well in ideal situations, they also face challenges arising from practical issues such as nonlinearities and measurement noise. These issues are also discussed through intensive numerical simulations.

Following the procedure described in Section 4, a GH is designed so that the discrete-time system has zeros at $z_{d,1}=0.1$ and $z_{d,1}=0.9$. We use a piecewise constant hold function with three subintervals ($N=3$) and the hold gain is $h={[2.155,-0.577,1.422]}^{\top }$. The sampling time is given by $T_s=0.1$ s. Figure 8 shows the behavior of the system under the ZDA (11) that has been designed using the system parameters $S_{\mathrm{d}}$, $g_{\mathrm{d}}$, and ${\psi }_{\mathrm{d}}$ assuming that ZOH and SS are used to interface analog and digital signals. The attack is injected when the system is at a steady state (rated power point). It is seen that the generator angular velocity ${\omega }_g\left(t\right)$ starts oscillating with increasing magnitude and this is captured by the sampled output ${\omega }_g\left(k{T}_s\right)$ (denoted as ${\omega }_{g,\mathrm{GH}}$ in the figure) when GH is used, while the sampled output under ZOH remains almost unchanged. It is noted that the signal ${\omega }_{g,\mathrm{ZOH}}$ also diverges as time goes to infinity but very slowly compared to ${\omega }_{g,\mathrm{GH}}$, which means that it is practically meaningless to use ${\omega }_{g,\mathrm{ZOH}}$ as a monitoring signal for the purpose of attack detection. This comes from the nonlinearity of the wind turbine; the diverging attack signal makes the state variables escape the region where the linear approximation is valid.

It is well known that even though a GH can shift the zeros to desired locations, it may induce a violent transient between sampling instants [20]. Typically, this can happen when the pattern associated to the GH has a large transition. For example, consider the GH designed above and suppose that the control input generated by the GH using $u_k={\overline{T}}_g$, shown in Figure 9a, is applied to the wind turbine. Then, as can be seen in Figure 9b, the sampled output ${\omega }_{g,\mathrm{GH}}$ seems to converge to a constant, but its continuous-time counterpart oscillates severely, and this can happen even if no attack is injected.

The possibly undesirable intersample behavior can be avoided by using GS instead of SS at the output side and using ZOH at the input side. To demonstrate this, we follow the design described in Section 4 to obtain a GS with $N=4$ and $w=[-28.252,58.5731,-70.906,41.5846]$ so that the zeros are placed at $z_{d,1}=0.1$, $z_{d,2}=0.9$, and $z_{d,3}=0$. Figure 10 shows that the presence of ZDA can be detected by monitoring the signal ${\omega }_{g,\mathrm{GS}}$, which is the output of the GS. In the simulation, the same ZDA injected in the case of GH is used. It is emphasized that the behavior of the internal states are the same as the case with ZOH and SS shown in Figure 7, and it is free of violent intersample behavior possibly induced by a GH.

In practice, the noise of measurements is always present and if the weights of GS are very large, then the noise will be amplified, leading to a false alarm. To demonstrate this, suppose that the measurement ${\omega }_g\left(t\right)$ is contaminated by noise and consider two designs of GS: GS1 with $w=[-28.252,58.5731,-70.906,41.5846]$ and $z_{d,1}=0.1$, $z_{d,2}=0.9$, $z_{d,3}=0$; GS2 with $w=[-8.924,18.983,-22.903,13.844]$ and $z_{d,1}=0.1$, $z_{d,2}=0.7$, $z_{d,3}=0$. Figure 11 shows the effect of the measurement noise under the same setting of Figure 10. The sampled output of GS ${\breve{y}}_k$ is denoted by ${\omega }_{g,\mathrm{GS}}$ in the figure. In Figure 11a, some of the sampled outputs of GS ${\omega }_{g,\mathrm{GS}1}$ exceed the threshold although no attack signal is injected, leading to a false alarm. On the contrary, Figure 11b shows that GS2 with relatively smaller weights is less affected by measurement noise.

Since the sampler weights depend on the location of desired zeros, we numerically investigate how they are related. With N and $z_{d,3}$ fixed as $N=4$ and $z_{d,3}=0$, the desired zeros $z_{d,1},z_{d,2}$ are selected from $(-1,1)\times (-1,1)$ and the corresponding sampler weights are determined. The result is shown in Figure 12a. In addition, by doing the simulation with the designed GS, we count the number of false alarms, see Figure 12b. It is observed that the number of false alarms is roughly proportional to the norm of sampler weights, and this explains why GS2 is less sensitive to measurement noise.

6. Conclusions

In this paper, we study the security problem on wind turbines that are controlled and monitored through the communication network. It is shown that at the rated power point, the linearized discrete-time model of the wind turbine has an unstable zero for a range of sampling periods, which means that wind turbines that are digitally controlled are vulnerable to zero-dynamics attacks. In order to increase security against ZDA, two countermeasures based on generalized hold and generalized sampler have been proposed with detailed design procedures, and through numerical simulations, it is shown that these approaches make ZDA ineffective. Practical issues such as nonlinearities and measurement noise are discussed in detail.

We are currently working on a robust and optimal design of the proposed strategies, which are challenging research topics. Validation of the proposed approaches using more realistic models or software and simultaneous design of two components are also interesting future research topics.