Extension of the Langevin power curve analysis by separation per operational state

In the last few years, the dynamical characterization of the power output of a wind turbine by means of a Langevin equation has been well established. For this approach, temporally highly resolved measurements of wind speed and power output are used to obtain the drift and diffusion coefficients of the energy conversion process. These coefficients fully determine a Langevin stochastic differential equation with Gaussian white noise. We show that the dynamics of the power output of a wind turbine have a hidden dependency on the turbine's different operational states. Here, we use an approach based on clustering Pearson correlation matrices for different observables on a moving time window to identify different operational states. We have identified five operational states in total, for example, the state of rated power. Those different operational states distinguish non-stationary behavior in the mutual dependencies and represent different turbine control settings. As a next step, we condition our Langevin analysis on these different states to reveal distinctly different behaviors of the power conversion process for each operational state. Moreover, in our new representation, hysteresis effects which have typically appeared in the Langevin dynamics of wind turbines seem to be resolved. We assign these typically observed hysteresis effects clearly to the change of the wind energy system between our estimated different operational states.

1 Introduction

Wind turbines have become a significant source of renewable energy due to their environmentally friendly nature and potential to generate electricity (Ackermann and Söder, 2002; Hannan et al., 2023). Analyzing their energy conversion process (Yaramasu et al., 2015) can be challenging due to the intricate nature influenced by various factors, including wind speed, turbulence, and mechanical wear. Nonetheless, a precise understanding of the dynamics of wind turbines is critical to simulate and consequently optimize their energy output and detect malfunctions (Wächter et al., 2011).

Recently, the Langevin equation approach has been used to study the dynamics of wind turbines (Tabar, 2019; Mücke et al., 2015; Milan et al., 2010; Raischel et al., 2013; Lind et al., 2017). To capture the dynamics of the energy conversion process using this approach, highly resolved temporal measurements of wind speed and power output are employed to determine the drift and diffusion coefficients. These coefficients characterize both the deterministic and stochastic behaviors of the system. However, this approach assumes a quasi-stationary system and does not consider the potential impact of different operational states.

We employ k-means clustering, an established unsupervised machine method to a set of Pearson correlation matrices. This type of approach was first put forward in Münnix et al. (2012) for financial correlation matrices and is now also used in other fields, such as traffic science (Wang et al., 2021, 2022, 2023). The method was further shaped and extended in Heckens et al. (2020), Pharasi et al. (2020), Pharasi et al. (2021), and Heckens and Guhr (2022). Here, we employ recent progress on clustering structures for correlation matrices of wind turbines (Bette et al., 2023; Jungblut et al., 2025). By analyzing different observables over a moving time window, they identified various operational states, which distinguish non-stationary behavior in the mutual dependencies and represent different turbine control settings. The dynamics of these operational states have been studied in Bette et al. (2023), utilizing a Langevin ansatz to comprehend the deterministic behavior, enhancing the comprehension of the transition between these states.

This study employs the method of Bette et al. (2023) to identify and distinguish various operational states of the wind turbine, which are then used to condition the Langevin analysis. The analysis reveals unique behavior patterns in the power conversion process corresponding to each operational state and thus also successfully resolves hysteresis effects commonly observed in the power conversion process of wind turbines (Mücke et al., 2015; Lin et al., 2023).

2 Hysteresis effects in the power conversion process of wind turbines

Hysteresis effects are a significant phenomenon observed in the operational dynamics of wind turbines, particularly evident in high-frequency Supervisory Control and Data Acquisition (SCADA) data during the power conversion process. These effects arise primarily from two interrelated factors: the switching of operational states and the stochastic nature of wind.

In order to ensure the safe operation of wind turbines, control strategies are implemented which actively regulate the power output and trigger different operating states depending on the most recent history of the same wind speed. To avoid instability in the view of rapid fluctuations, robust control strategies induce delayed switching between a finite number of operating states.

The switching of operational states, which occurs in response to varying wind conditions, leads to a nonlinear response in turbine performance. For instance, when wind speed fluctuates, turbines may transition between modes such as cut-in, rated, and cut-out states. This switching can result in distinct performance paths during wind acceleration and deceleration. Consequently, the system exhibits hysteresis, where the output power does not return along the same trajectory when wind conditions reverse, leading to discrepancies in expected versus actual power generation.

In our analysis, we specifically condition the dynamics on the operational state rather than examining the transitions between them. This approach is crucial, as it reveals that the dynamics of the power conversion process are intrinsically linked to the operational state in which the turbine is functioning. Without conditioning on operational states, the analysis of the dynamics of the power conversion process yields a mixed average of dynamics across different states, effectively diluting the unique characteristics associated with each operational condition. By separating the data according to operational states, we can derive the distinct dynamics pertinent to each state. This granularity enables us to identify multiple stable fixed points at various wind speeds, corresponding to different operational states. Importantly, this methodology clarifies the hysteresis effects induced by the switching of operational states. By isolating each state, we can observe how these stable fixed points shift in response to changes in wind conditions, thereby resolving the complexities of hysteresis that may otherwise obscure our understanding.

The random fluctuations of the wind can trigger noise-induced transitions, causing the system to shift from one stable state to another, even when the drift term indicates stability in the original state.

3 Dataset

The data utilized in this study are sourced from the SCADA system of a Vestas V90 turbine located in the Thanet offshore wind farm. These measurements were recorded at approximately 5 s intervals throughout the year 2017. To ensure consistent time stamps and a stable frequency, the data were aggregated by averaging over 10 s intervals. It is important to note that if no measurements were obtained within the original 5 s interval, the aggregated dataset may contain missing data during the corresponding 10 s interval. In this study we rescaled the ActivePower and the WindSpeed values.

The dataset under analysis comprises six variables, namely the following:

• ActivePower – generated active power;

• CurrentL1 – generated current (chosen from one of the three phases due to no deviations in the data);

• RotorRPM – rotations per minute of the rotor;

• GeneratorRPM – rotations per minute of the high-speed shaft at the generator;

• BladePitchAngle – blade pitch angle of the blades;

• WindSpeed – wind speed.

Our expectation for the V90 turbine is a shift in control strategy as the wind speed changes. This shift includes transitioning from a low wind speed regime with variable rotation speed and increasing power to an intermediate regime with constant rotation and increasing power and finally to a rated region with constant rotation and constant power. The selection of the above six variables allows us to effectively track these operational state changes.

4 Correlation matrix states

In order to track the operational state of a turbine, we employ a method presented in Bette et al. (2023). Pearson correlation matrices are calculated for moving windows with non-overlapping time intervals called epochs to obtain a time series of correlation matrices. These are clustered to find structurally different operational states and thereby a time series S(t), which offers us the current operational state.

For this calculation we use all variables presented in Sect. 3. Each of these variables is represented by a time series X_l(t), where $l=\mathrm{1},\mathrm{\dots },\mathrm{6}$ represents the different variables, and $t=\mathrm{1},\mathrm{\dots },T_{\text{end}}$ is the time variable. To capture non-stationarity we separate the whole time series into disjoint intervals λ of length T=30 min.

The 30 min time span represents a compromise. The choice of T depends on the specific system being analyzed, as it must balance two competing factors: accuracy and statistical uncertainty. A larger T improves statistical reliability by averaging over more data points, reducing noise. However, longer time windows can obscure short-term dynamics and fail to capture rapid changes in the system. Conversely, smaller T values provide higher temporal resolution but may introduce greater statistical uncertainty due to fewer data points. In our approach, we select T based on a trade-off that best captures the system's dynamics while maintaining a reasonable level of statistical accuracy. Given that external factors, such as wind, can change on timescales ranging from minutes to hours, shorter epochs are required to capture the non-stationarity. Such trade-offs are common when working with time series of correlation matrices. ε represents the starting time of an interval λ.

Next, we normalize each time series in every interval to zero mean and a standard deviation of 1 by

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{\displaystyle }{\mathit{\mu }}_l\left(t\right)=〈X_l\left(t^{\prime }\right)|\mathit{\epsilon }\le t^{\prime }<\mathit{\epsilon }+T〉,l=\mathrm{1},\mathrm{\dots },L,t\in \mathit{\lambda }\text{(2)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{\mathit{\sigma }}_l\left(t\right)=& {\displaystyle }\sqrt{〈\left(X_l\right(t^{\prime })-{\mathit{\mu }}_l(t^{\prime }\left){)}^{\mathrm{2}}\right|\mathit{\epsilon }\le t^{\prime }<\mathit{\epsilon }+T〉}\\ {\displaystyle }& {\displaystyle }l=\mathrm{1},\mathrm{\dots },L,t\in \mathit{\lambda }\end{array}\text{(3)}& {\displaystyle }& {\displaystyle }{\displaystyle }{G}_l\left(t\right)={\displaystyle \frac{X_l\left(t\right)-{\mathit{\mu }}_l\left(t\right)}{{\mathit{\sigma }}_l\left(t\right)}},l=\mathrm{1},\mathrm{\dots },K,\end{array}$$

with μ_l(t) and σ_l(t) being the mean value and the standard deviation of variable l in the interval λ with a starting time ε. By arranging the variable time series in each epoch in a L×T data matrix,

$$\begin{array}{}\text{(4)}& \mathbf{G}\left(\mathit{\lambda }\right)=\left[\begin{array}{ccc}{G}_{\mathrm{1}}\left(\mathit{\epsilon }\right)& \mathrm{\dots }& G_{\mathrm{1}}(\mathit{\epsilon }+T-\mathrm{1})\\ \mathrm{⋮}& & \mathrm{⋮}\\ G_k\left(\mathit{\epsilon }\right)& \mathrm{\ddots }& G_k(\mathit{\epsilon }+T-\mathrm{1})\\ \mathrm{⋮}& & \mathrm{⋮}\\ G_K\left(\mathit{\epsilon }\right)& \mathrm{\dots }& G_K(\mathit{\epsilon }+T-\mathrm{1})\end{array}\right].\end{array}$$

We calculate the correlation matrix in the interval λ:

$$\begin{array}{}\text{(5)}& \mathbf{C}\left(\mathit{\lambda }\right)={\displaystyle \frac{\mathrm{1}}{T}}\mathbf{G}\left(\mathit{\lambda }\right){\mathbf{G}}^{†}\left(\mathit{\lambda }\right).\end{array}$$

Here, G^†(λ) denotes the transpose of G(λ). Each matrix element C_ij(λ) is the Pearson correlation coefficient between the variables i and j in the interval λ.

We apply hierarchical k-means clustering to find recurring states in our system. The algorithm is a divisive clustering that splits by applying standard k-means with k=2. In each step, the cluster with the largest internal distance to its own center is split. Hence, we must define a distance $d(\mathit{\lambda },{\mathit{\lambda }}^{\prime })$ between the correlation matrices for intervals λ and λ^′:

$$\begin{array}{}\text{(6)}& d(\mathit{\lambda },{\mathit{\lambda }}^{\prime })=\sqrt{\sum _{i,j}\left(C_{ij}\right(\mathit{\lambda })-C_{ij}({\mathit{\lambda }}^{\prime }){)}^{\mathrm{2}}}=\left|\right|\mathbf{C}\left(\mathit{\lambda }\right)-\mathbf{C}\left({\mathit{\lambda }}^{\prime }\right)\left|\right|.\end{array}$$

The center of cluster s is calculated as the element-wise mean ${\stackrel{\mathrm{‾}}{C}}_{ij}\left(s\right)=\langle C_{ij}\left(\mathit{\lambda }\right)\left|S\right(\mathit{\lambda })=s)\rangle$. S(t) is the function which results in the cluster s for any time stamp t assigned by the algorithm to the interval λ that contains t. A more detailed description of the clustering procedure is found in Bette et al. (2023).

We visualize the different states of the power output for our dataset in Fig. 1.

Figure 1 Visualization of the data and the operational states: scatter plot with color-coded markers based on the operational state and contour lines based on the density (see Appendix A) to show the borders of each operational state.

An open question is the possible existence of operational states of complete wind farms, which could be investigated following a similar methodology. The results could potentially support the operation of wind farms, namely wind farm control.

While we focused on linear correlations in this study, which offer a simple and effective first-order approach, we acknowledge that the system's nonlinearity may suggest the potential for more advanced methods. Techniques such as mutual information or Granger causality could provide deeper insights into the relationships between variables.

5 Estimation of the Kramers–Moyal coefficients

In the context of stochastic differential equations, the ensemble dynamics of a stochastic process can be described using the Kramers–Moyal coefficients. The first coefficient represents the drift, which indicates the average tendency of the system to move in a certain direction. This drift reflects trends that would be considered the backbone in a deterministic model approach. The second coefficient corresponds to the diffusion, capturing the variance of the random fluctuations or noise in the system. This diffusion coefficient quantifies how much the possible states of the process can spread out over time due to these random influences.

We start with the traditional approach to model the power conversion process,

$$\begin{array}{}\text{(7)}& \dot{P}\left(t\right){|}_{u\left(t\right)=u}=D_{\mathrm{P}}^{\left(\mathrm{1}\right)}\left(P\right(t),u)+\sqrt{D_{\mathrm{P}}^{\left(\mathrm{2}\right)}\left(P\right(t),u)}\cdot \mathrm{\Gamma }\left(t\right),\end{array}$$

of a wind turbine in terms of a stationary Langevin equation (Risken, 1996; Wächter et al., 2011; Raischel et al., 2013; Tabar, 2019). Here, the power output P(t) is modeled as a 1-dimensional stationary stochastic process for a fixed wind speed u. We assume a Gaussian-distributed, delta-correlated noise Γ(t) with zero mean and a variance of 2.

Analytically, the nth order conditional moments of the power output,

$$\begin{array}{}\text{(8)}& M_{\mathrm{P}}^{\left(n\right)}(P,u,\mathit{\tau })=\langle \left({\mathrm{\Delta }}_{\mathit{\tau }}P\right(t){)}^n\rangle {|}_{P\left(t\right)=P,u\left(t\right)=u},\end{array}$$

can be derived with expectation value of the increments to the power of n ${\mathrm{\Delta }}_{\mathit{\tau }}P(t{)}^n=(P(t+\mathit{\tau })-P\left(t\right){)}^n$ over the time step τ at the specific state (P,u) (Risken, 1996; Tabar, 2019).

With the nth conditional moments, the nth Kramers–Moyal coefficient,

$$\begin{array}{}\text{(9)}& D_{\mathrm{P}}^{\left(n\right)}(P,u)=\underset{\mathit{\tau }\to \mathrm{0}}{lim}{\displaystyle \frac{M_{\mathrm{P}}^{\left(n\right)}(P,u,\mathit{\tau })}{n\mathrm{!}\cdot \mathit{\tau }}},\end{array}$$

can be calculated (Risken, 1996; Tabar, 2019).

We consider a 2-dimensional dataset (P, u) with N data points with a uniform sampling interval ${\mathit{\tau }}_{\mathrm{s}}=\mathrm{1}/f_{\mathrm{s}}$. Furthermore, we define ${\mathit{\tau }}_m=m\cdot {\mathit{\tau }}_{\mathrm{s}}$, where m∈ℕ. With this definition, we can calculate the increments of the power output over a time lag τ_m with

$$\begin{array}{}\text{(10)}& {\mathrm{\Delta }}_{{\mathit{\tau }}_m}{P}_i=P_{i+m}-P_i.\end{array}$$

To estimate the nth conditional moment $M_{\mathrm{P}}^{\left(n\right)}(P,u,{\mathit{\tau }}_m)$, we employ the Nadaraya–Watson estimator with a 2-dimensional kernel $K_{a,b}(x,y)=k_a\left(x\right)\cdot k_b\left(y\right)$ (Nadaraya, 1964; Epanechnikov, 1969; Silverman, 1986). This kernel can be represented as the product of two 1-dimensional kernels. For the first conditional moment, we can calculate the weighted average of the power output increments using

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\displaystyle }{\widehat{M}}_{\mathrm{P}}^{\left(n\right)}(P,u,{\mathit{\tau }}_m)=& {\displaystyle }{\sum }_{i=\mathrm{1}}^{N-m}({\mathrm{\Delta }}_{{\mathit{\tau }}_m}{P}_i{)}^n\\ {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{K_{a,b}\left(\frac{P_i-P}{h_{\mathrm{P}}},\frac{u_i-u}{h_{\mathrm{u}}}\right)}{{\sum }_{j=\mathrm{1}}^{N-m}{K}_{a,b}\left(\frac{P_j-P}{h_{\mathrm{P}}},\frac{u_j-u}{h_{\mathrm{u}}}\right)}}.\end{array}\end{array}$$

These weights in our specific case are determined by the states P and u, as well as the kernel functions k_a(x) and k_b(x), along with the bandwidths for power output (h_P) and wind speed (h_u).

There are plenty of different kernel functions which are useful for different scenarios.

The Epanechnikov, Gaussian, and rectangular kernels are three commonly used kernel functions in non-parametric estimation and smoothing techniques (Epanechnikov, 1969). Each of these kernels has distinct properties that impact their use and the resulting estimation or smoothing outcomes. We summarize the three most commonly used kernels in Table 1.

Table 1 Description of different kernel functions.

The choice between these kernels depends on the specific characteristics of the data and the desired properties of the estimation or smoothing procedure (Wied and Weißbach, 2012). The Epanechnikov kernel,

$$\begin{array}{}\text{(12)}& k_{\text{E}}\left(x\right)=\left\{\begin{array}{ll}\mathrm{1}-x^{\mathrm{2}},& \left|x\right|\le \mathrm{1}\\ \mathrm{0},& \left|x\right|>\mathrm{1},\end{array}\right.\end{array}$$

is often favored when robustness and efficiency are important and when a localized smoothing effect is desired. The Gaussian kernel is popular for its smoothness and computational efficiency. The rectangular kernel is suitable when simplicity and computational efficiency are prioritized, but it may not handle outliers or extreme values as the other kernels. In this study, we use a Epanechnikov kernel function (Eq. 12).

At least as important as the kernel function is the related bandwidth (Nadaraya, 1965; Jones et al., 1996; Scott, 2010; Silverman, 1986). For the analysis of large structures (macro-scale structures), large bandwidths should be used. However, with larger bandwidths, the small structures (micro-scale structures) are no longer visible. For a more comprehensive understanding of bandwidth selection for estimating Kramers–Moyal coefficients, we recommend consulting the study (Wiedemann et al., 2024). As for this study, 2- and 3-dimensional estimation is necessary; a rigorously optimized bandwidth such as in Wiedemann et al. (2024) is unfortunately not applicable. We therefore used the bandwidths according to the IEC 61400-12-1 (IEC, 2005). The bandwidth h_u for the wind speed is 1 m s⁻¹, and the bandwidth for the power h_P is 100 kW. These bandwidths should be adjusted based on the given dataset (larger bandwidths for a smaller dataset, smaller bandwidths for a larger dataset). For the dataset we used, we found that these bandwidths, in conjunction with the Epanechnikov kernel, yield fully reliable results.

We assume that the nth conditional moments $M^{\left(n\right)}(P,u,\mathit{\tau })$ are linear for small time steps τ. We can estimate the Kramers–Moyal coefficients,

$$\begin{array}{}\text{(13)}& {\widehat{D}}_{\mathrm{P}}^{\left(n\right)}(P,u)={\displaystyle \frac{\mathrm{1}}{M}}\sum _{m=\mathrm{1}}^M{\displaystyle \frac{{\widehat{M}}^{\left(n\right)}(P,u,{\mathit{\tau }}_m)}{n\mathrm{!}\cdot {\mathit{\tau }}_m}},\end{array}$$

by averaging the nth conditional moments divided by the used (small) time step τ_m times n factorial (Tabar, 2019). It can also make sense to give smaller m values higher weight. We employed M=3 for this particular dataset, yielding meaningful outcomes.

We can determine the fixed points P₀(u) of the system (Wächter et al., 2011). These fixed points correspond to values of P at which the drift term becomes zero:

$$\begin{array}{}\text{(14)}& {\widehat{D}}_{\mathrm{P}}^{\left(\mathrm{1}\right)}(P_{\mathrm{0}},u)=\mathrm{0},\end{array}$$

indicating an equilibrium state.

Furthermore, in order to assess the stability of these fixed points, we examine the derivative of the drift at the fixed point. If the derivative is negative, it signifies that the fixed point is stable:

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}P}}{\widehat{D}}_{\mathrm{P}}^{\left(\mathrm{1}\right)}(P_{\mathrm{0}},u)<\mathrm{0}.\end{array}$$

The derivative of the drift at the fixed point plays a crucial role in understanding the stability of the fixed point as well as providing valuable insights into the mean reversal time.

When studying the stability of a fixed point, we are interested in how the system responds to small perturbations of its equilibrium state. The derivative of the drift provides information about the local behavior of the system near the fixed point.

As said before, if the derivative of the drift evaluated at the fixed point is negative, it indicates that the fixed point is stable. In this case, any small disturbances from the equilibrium will eventually dampen out, and the system will return to its steady state. On the other hand, if the derivative is positive, the fixed point is unstable, and even the slightest perturbations will cause the system to diverge from the equilibrium. Additionally, the derivative of the drift defines the mean reversal time of a system.

We make the assumption that the operational states S(t) of the wind turbine can only take discrete values, specifically $S\left(t\right)\in [\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5}]$ (which is related to the five identified clusters). Furthermore, we consider that both the drift and diffusion coefficients depend on the turbine's operational state. By incorporating this additional condition, we can reformulate the Langevin equation for the power conversion process using

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\displaystyle }\dot{P}\left(t\right){|}_{u\left(t\right)=u,S\left(t\right)=S}=& {\displaystyle }{D}_{\mathrm{P}}^{\left(\mathrm{1}\right)}\left(P\right(t),u,S)\\ {\displaystyle }& {\displaystyle }+\sqrt{D_{\mathrm{P}}^{\left(\mathrm{2}\right)}\left(P\right(t),u,S)}\cdot \mathrm{\Gamma }\left(t\right).\end{array}\end{array}$$

The numerical approach can be derived in a similar manner as before. The only distinction is that we employ a 3-dimensional kernel $K_{a,b,c}(x,y,z)=k_a\left(x\right)\cdot k_b\left(y\right)\cdot k_c\left(z\right)$. Due to the discrete values of the operational state, we can utilize a dedicated Boolean kernel function.

$$\begin{array}{}\text{(17)}& k_{\text{Bool}}\left(x\right)=\left\{\begin{array}{ll}\mathrm{1}& x=\mathrm{0}\\ \mathrm{0}& x\ne \mathrm{0}\end{array}\right.\end{array}$$

We use

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }{\widehat{M}}_{\mathrm{P}}^{\left(n\right)}(P,u,S,{\mathit{\tau }}_m)={\sum }_{i=\mathrm{1}}^{N-m}({\mathrm{\Delta }}_{{\mathit{\tau }}_m}{P}_i{)}^n\\ {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{K_{a,b,\text{bool}}\left(\frac{P_i-P}{h_{\mathrm{P}}},\frac{u_i-u}{h_{\mathrm{u}}},S_i-S\right)}{{\sum }_{j=\mathrm{1}}^{N-m}{K}_{a,b,\text{bool}}\left(\frac{P_j-P}{h_{\mathrm{P}}},\frac{u_j-u}{h_{\mathrm{u}}},S_j-S\right)}}\end{array}\end{array}$$

to estimate the nth conditional moment at a specific state (P, u, S). With these conditional moments, we are able to obtain the Kramers–Moyal coefficients in a similar manner to that shown above.

6 Stochastic analysis of the power conversion process

In this section, we elucidate the outcomes of our investigation into the wind turbine power conversion process using the Kramers–Moyal coefficients, considering scenarios both with and without separation per operational state. Our primary focus centers on analyzing the drift and diffusion values governing the power output of a wind turbine. To deepen our understanding, we extend the analysis to include the computation of fixed points, their associated stability, and the diffusion values at these fixed points, revealing the nuanced dynamics intrinsic to diverse operational states.

The calculated drift values of the power output, as depicted in Fig. 2, reveal a familiar pattern observed in prior studies without operational state separation (Wächter et al., 2011; Mücke et al., 2015; Milan et al., 2010). The top-left plot illustrates typical behavior in the power conversion process. Significant differences emerge when comparing drift maps for distinct operational states. Notably, a clear contrast is evident when analyzing state 2 against state 4, particularly at rescaled wind speeds (u) of approximately 0.4−0.5. Similar variations are observed when comparing state 4 and state 5.

Figure 2 This figure illustrates the drift maps of the power conversion process, categorized by turbine states. The drift values are depicted using a color-coded scheme, with a corresponding color bar provided on the right-hand side for reference.

To deepen the analysis, we calculate stable fixed points and their derivatives. Figure 3 depicts stable fixed points per wind speed, highlighting disparities, particularly in regions characterized by rescaled u values around 0.4−0.6 across different operational states. Multiple stable fixed points are identified for a given wind speed, with states 1 and 2 displaying relative similarity. In contrast, significant differences are observed in other states, confirming the presence of hysteresis effects within the system dynamics (Mücke et al., 2015; Lin et al., 2023). The absence of multiple fixed points per wind speed without operational state separation is attributed to the choice of a relatively high bandwidth during the estimation process for Kramers–Moyal coefficients. This, coupled with the use of a kernel function and the distribution of operational states, may have led to a more aggregated representation of the system dynamics.

Figure 3 Rescaled stable fixed points P₀ of the power output in relation to the rescaled wind speed. The fixed points are categorized based on the operational states, with each state distinguished by a distinct color.

We further explore the stability of these fixed points through the derivatives if the drift at the fixed points. Figure 4 illustrates the derivatives of the drift at stable fixed points per wind speed. Negative values signify stable fixed points, with larger absolute values indicating a shorter mean reversal time towards the fixed point. Comparing derivatives for different operational states reveals similarities for states 1, 2, 3, and 4 across rescaled u values of approximately 0.0−0.6. However, a significant change occurs for state 4 around u≈0.6, aligning it with state 5. Notably, at $u\approx \mathrm{0.5}-\mathrm{0.6}$, clear differences emerge, with values for state 5 consistently lower than those for other states.

Figure 4 Derivative of the drift at the stable fixed points of the power output in relation to the rescaled wind speed. The derivatives are categorized based on the operational states, with each state distinguished by a distinct color.

We extend our analysis from the deterministic parts of the behavior the calculation of the diffusion coefficients. The results are presented in Fig. 5. Without operational state separation, diffusion values are generally smaller near fixed points than towards the edges. Small diffusion values are observed at rated wind speeds (u>0.5) and lower power values (P<0.4). Differences between the diffusion values for different states are identified across various wind speeds, with states 1 and 2 exhibiting similarities, particularly around u between 0.4 and 0.6.

Figure 5 Diffusion maps of the power conversion process, categorized by turbine states. The diffusion values are depicted using a color-coded scheme, with a corresponding color bar provided on the right-hand side for reference.

In dynamical systems, the presence of a diffusion term can significantly impact the location of stable fixed points identified through the drift term. While the drift term typically defines the deterministic dynamics of the system, the diffusion term accounts for random fluctuations and noise, which can induce transitions between these fixed points.

When noise is introduced, it may alter the effective landscape of the system, potentially shifting the positions of stable fixed points. This phenomenon can lead to noise-induced transitions, where the system, under the influence of stochastic perturbations, may escape from one stable state and transition to another, even if the drift term suggests stability in the original state.

Further analysis involves the calculation of diffusion values at the stable fixed points conditioned on wind speeds, visualized in Fig. 6. Notably, here, we compare the diffusion values for the same wind speeds but different power values. States 1 and 2 exhibit similarity, while significant differences are observed across all other states, especially for u between 0.4 and 0.6. The diffusion values for state 5 are smaller than those for other states at u smaller than 0.6.

Figure 6 Diffusion values at the fixed points of the power output in relation to the Rescaled wind speed. The diffusion values are categorized based on the operational states, with each state distinguished by a distinct color.

The diffusion analysis further underscores distinctions in the behavior of state 5 compared to other states, particularly at wind speeds below 0.6. The smaller diffusion values for state 5, coupled with the reduced derivatives at the fixed points, contribute to diminished fluctuations around these stable fixed points of the power time series of state 5 in comparison to other states at wind speeds below 0.6. In contrast, state 2 exhibits higher diffusion values at the fixed points for wind speeds between 0.4 and 0.55 than the other states. Having high diffusion values at the fixed points and similarities in derivative values with states 3 and 5 results in higher fluctuations for these wind speeds when compared to all other states.

7 Conclusions

In this study, we successfully applied a method to estimate the dynamics of the power conversion process while taking into account different operational states, identified using a correlation matrix algorithm (Bette et al., 2023) to take non-stationarity into account. Our analysis revealed distinct dynamics associated with each operational state in the power conversion process, emphasizing the significant influence of these states on the overall behavior of the system.

We successfully resolved hysteresis effects within the power conversion process. When separating per operational state, distinct fixed points per wind speed are visible. Without accounting for states, these are averaged out into one fixed point per wind speed. The presented analysis also allows us to identify differences in the dynamic behavior of states. State 5, representing rated power production, displayed a much more stable behavior, with fewer fluctuations than other states. This even remained true for wind speed values where state 5 overlaps with other states.

The results clearly show that it is possible to enhance existing methods by considering the operational states described. The analysis concept does not need to change much but rather only takes the automatically detected operational state as a distinction parameter for multiple subanalyses with the original method.

Appendix A: Operational state contour lines

To represent the distribution of operational states visually, we utilize kernel density estimation given by

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{\displaystyle }\widehat{f}(u_{\mathrm{0}},P_{\mathrm{0}},S_{\mathrm{0}})=& {\displaystyle }{\sum }_{i=\mathrm{1}}^N{K}_{\text{gauss},\text{gauss},\text{bool}}\\ {\displaystyle }& {\displaystyle }\cdot \left({\displaystyle \frac{P_i-P_{\mathrm{0}}}{h_{\mathrm{P}}}},{\displaystyle \frac{u_i-u_{\mathrm{0}}}{u_h}},S_i-S_{\mathrm{0}}\right).\end{array}\end{array}$$

Contour lines are then generated using the formula:

$$\begin{array}{}\text{(A2)}& \stackrel{\mathrm{̃}}{f}(u_{\mathrm{0}},P_{\mathrm{0}},S_{\mathrm{0}})=\left\{\begin{array}{ll}\mathrm{1},& \mathit{\rho }(u_{\mathrm{0}},P_{\mathrm{0}},S_{\mathrm{0}})\ge {\mathit{\rho }}_{\mathrm{0}}.\\ \mathrm{0},& \text{else}\end{array}\right.\end{array}$$

Here, ρ₀ is a predefined threshold (set at 20).