# Wind turbulence modeling for real-time simulation

Mohamed Hajjem, Stéphane Victor, Pierre Melchior, Patrick Lanusse, Lara Thomas

► To cite this version:

Mohamed Hajjem, Stéphane Victor, Pierre Melchior, Patrick Lanusse, Lara Thomas. Wind turbulence modeling for real-time simulation. Fractional Calculus and Applied Analysis, 2023, 26 (4), pp.1632-1662. 10.1007/s13540-023-00165-0. hal-04317726

HAL Id: hal-04317726

https://hal.science/hal-04317726v1

Submitted on 12 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Wind turbulence modeling for real-time simulation

Mohamed Hajjem¹,² · Stéphane Victor¹ ·

Pierre Melchior¹ · Patrick Lanusse¹ ·

Lara Thomas²

**Abstract** This paper proposes a method to design a real-time wind turbulence simulator. The objective of such a simulator is to make dynamical models more accurate in order to develop more robust controls, especially in the case of mechanical systems operating outdoors such as tracking antennas. The models used to generate a random wind speed are based on the wind spectral characteristics rather than time domain ones. Indeed, turbulence is a stochastic and non-stationary process corresponding to the short-term component of wind and therefore is difficult to model in time domain. Wind spectral characteristics are described by the power spectral density whose approximation is used in the real-time simulator to reproduce wind behavior. The von Kármán model is the most commonly used model to approximate the power spectral density of wind turbulence. However, this model was originally designed for aircraft, and therefore for high altitudes and moving systems. In our case, tracking antennas are used under different conditions: low altitude and slow moving systems. This makes the von Kármán model less precise in middle and high range frequency. Consequently, there is a need in more accurately modeling wind turbulence under these specific conditions. As von Karman's model expression uses fractional calculus, other models with the same range of parameter number are proposed to better approximate the wind power spectral density by using Cole-Cole and Davidson-Cole fractional models. The idea is to use shaping filters issued from these more precise fractional models to generate realistic random wind speed turbulence from a random white noise input. These filters are approximated by rational functions and their param-

<footnote>Mohamed Hajjem¹, Stéphane Victor¹,*, Pierre Melchior¹, Patrick Lanusse¹
Université de Bordeaux, IMS (UMR 5218 CNRS), 33405 Talence, France
E-mail: stephane.victor@ims-bordeaux.fr * corresponding author
E-mail: pierre.melchior@ims-bordeaux.fr, patrick.lanusse@ims-bordeaux.fr</footnote>

<footnote>Mohamed Hajjem², Lara Thomas²
Safran Data Systems, 33260 La Teste-de-Buch, France
E-mail: mohamed.hajjem@safrangroup.com, lara.thomas@safrangroup.com</footnote>
eters are tuned in the same way as those of the von Kármán model with the
same inputs (mean speed, standard deviation and length scale). If a rational
model is sought, a very high number of parameters (more than six times) is
needed to have the same precision as our proposed fractional Cole-Cole model
which has only four parameters. Finally, an accurate turbulence wind speed
generator is proposed.

**Keywords** Fractional models (primary) · Turbulence wind speed · Spectral analysis · Parameter estimation · Real-time simulation

# 1 Introduction

Many mechanical systems operating outdoors, such as parabolic antennas are continuously disturbed by wind turbulence. Therefore, it is essential to effectively reject wind effects in order to guarantee a high precision in path tracking. To achieve that, wind effects must be integrated into the system dynamical models with a realistic simulator capable of accurately reproducing the wind turbulence behavior, then, suitable control laws capable of rejecting them can be designed from these models. The study presented in this paper was conducted to improve the tracking performances of mobile antennas developed by Safran Data Systems [14]. Indeed, antennas are systems requiring a lot of precision, wind turbulence increases the errors and the signal loss risks, hence the need for a better compensation of wind effects on the system. The presented short horizon wind modeling can also be used in other fields. It is interesting for realistic dynamic behaviors simulations of other systems than antennas: wind turbines (blades, mechanical forces, ...) [12], wind response of solar panels, kites for traction or electric power generation [5] [4], sailing, etc.

The available wind studies are mostly focused on statistical and probabilistic aspects, especially for wind power analysis for wind energy conversion systems [17] [29]. For this type of modeling, the analyses are carried out over long periods of days and months and consequently only the long-term components are modeled. Therefore, these studies do not reproduce the wind impacts system dynamics, which mainly depends on turbulence and thus on the short-term component. An accurate model is then needed to well simulate realistic turbulence wind speed.

The short-term wind speed component, which corresponds to the turbulence, is a non-stationary stochastic process and thus is difficult to model in time domain. Therefore, the simulator design method presented in this paper is based on wind spectral characteristics. These characteristics will be described by models based on the well known von Kármán approximation [28] [19] and on the Davidson-Cole [7] and Cole-Cole [6] fractional functions. The spectral models are then used to design fractional order shaping filters capable of randomly generating wind speeds [9]. Davidson-Cole and Cole-Cole functions have been widely used to model other diffusive phenomenon such as heat transfers in lungs (see [10]). Havriliak-Negami functions can also be used to fit the power
spectral density (PSD) but they will be hard to simulate in time-domain (see for example [27]).

A first study was proposed in [13] where the approximation of PSD is pro-
posed by using simple fractional functions. Neither synthesis of shaping filters
nor any time domain simulations are provided. This article presents a complete
study with more complex and realistic approximations that takes into account
fractional transfer models which can be implemented in an antenna simulator
that includes wind turbulence disturbances. Moreover physical meaning is also
provided to the proposed model parameters in order to easily tune them for
realistic wind turbulence scenarios.

This paper is organized as follows. Section 2 will first present the dataset
that will be used to estimate the models. Section 3 is a state-of-the-art of
existing models to describe the wind, particularly the von Kármán model.
Section 4 presents the results of the von Kármán model identification and its
limitations. Then, Section 5 proposes improved fractional models based on
Cole-Cole and Davidson-Cole fractional models. Section 6 presents the design
of a shaping filter using the improved Cole-Cole model and the results of time
simulations. Finally, a conclusion is drawn in Section 7.

## 2 Experimental protocol

Wind spectral characteristics, on which the models will be based, will be used
as a reference for the identification. For this purpose, one-day wind speed
datasets recorded next to the city of Montpellier, France are used as a refer-
ence to validate the models. The sensor used to make the measurements is an
ultrasonic anemometer integrated into a Vaisala WXT520 weather transmitter
that measures wind speed and direction, temperature, air pressure, precipita-
tion, and relative humidity. The sensor has three identical transducers that
generate ultrasonic signals. The wind speed is calculated by measuring the
time it takes for the ultrasonic signal to travel from one transducer to the oth-
ers. The anemometer has a measurement range of 0 to 60 m/s and an accuracy
of ±3% (0 to 36m/s) and ±5% (36 to 60m/s).

In order to prevent having incoherent measurements, the anemometer is
configured to provide data every second with measurements that are averaged
over 20s. Such a configuration is imposed by Safran Data Systems in order to
detect long and strong gusts and to avoid any antenna damage. Indeed, the
Van der Hoven wind power spectrum [26] (see Figure 3) shows that above a
period of 20s, the turbulence power increases rapidly and therefore when this
power exceeds a certain limit, damage can occur and thus the antenna should
no longer be used and should remain in a safety position.

The moving average applied on the measurements increases the output
smoothness and attenuates rapid variations which modifies the high frequency
wind behavior. Therefore, in order to avoid the low pass filter effect of the
moving average and to be coherent with the original wind behavior, this study
will focus on spectral characteristics at frequencies below 0.1 Hz.
Wind turbulence modeling contains the short-term fluctuations of the wind speed, thus the signals will be studied over a period of 3 hours and not a whole day. Furthermore, the studied signals must have a quasi-constant average speed to avoid that long-term wind fluctuations interfere in the study of short-term ones. 28 signals were selected to be used as a reference for model identification and for comparison. Figure 1 presents one of the selected time signals of the measured wind speed.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Wind speed measurements during 3 hours (data used for identification)</figcaption></figure>

The wind simulator will be designed from models of wind turbulence in the frequency domain, therefore the reference data will be the PSD of the time-series data (Figure 1). Figure 2 represents the spectral distributions of the signal power, depending on the frequencies that compose them. The plain blue curve PSD, that corresponds to the time-series data of Figure 1, will be used for the parameter identification of the different approximation models.

# 3 State of the art

## 3.1 Wind speed components

The wind energy is distributed over a wide frequency band near the ground [3]. One of the first wind spectral analyses was made by Van Der Hoven in 1957 [26]. His results (see Figure 3) were obtained from measurements made at 100m height on the horizontal wind speed over a long period of time. They reveal several scales of fluctuations. The spectrum shows that the power spectral density can be divided into two components each with a maximum: one part corresponding to turbulence with short period variations, and the other part corresponding to long seasonal variations. The presence of three peaks can be
<figure><img src="image_3.png" /><figcaption>**Fig. 2** PSD of different wind speed measurements using Welch's method</figcaption></figure>

noticed. The turbulence peak has a period of about 76s and it is caused by gusts and surface effects. The diurnal peak has a period of about 12 hours and it corresponds to the variations in air mass circulation between day and night. Finally, the synoptic peak corresponds to a period of about 4 days and represents the variations due to large-scale atmospheric phenomena.

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Van der Hoven's wind spectrum [26]</figcaption></figure>

Wind speed can be considered as a stationary process in the energy gap ($15\text{min} < T < 3\text{h}$), contrary to turbulence ($T<10\text{min}$) where the speed follows a random and non-stationary process. Based on the Reynolds decomposition for fluids, the wind speed can be expressed as a sum of a steady component representing mean speed and a fluctuating component representing the turbu-
lence:

$$V(t) = \bar{V}(t) + V_{\text{turb}}(t) \qquad (3.1)$$

where $\bar{V}(t)$ is the long-term component and $V_{\text{turb}}(t)$ the turbulence component.

## 3.2 von Kármán model

The von Kármán model [28] [19] is a mathematical model of the power spectral densities of the linear and angular wind speed components. It approximates the PSD by a non-integer order function. It was developed by the United States Department of Defense based on on earlier work by von Kármán [8], for aircraft design and simulation applications. It is still widely used as in [18], [16] and [2]. The model can be used to design shaping filters which take white noise as input and return a stochastic signal with the spectral characteristics of the wind turbulence as output. The von Kármán model expresses the spectrum of the longitudinal, lateral and vertical components of the wind speed. The shape being similar for all three components, this study only focuses on von Karman's longitudinal model. The wind speed PSD according to the von Karman's longitudinal model is expressed as:

$$S(f) = \frac{4\sigma^2 \frac{L}{\bar{V}}}{\left(1 + 70.8\left(\frac{L}{\bar{V}}f\right)^2\right)^{5/6}}, \qquad (3.2)$$

where $f$ is the frequency, $\sigma$ is the turbulence standard deviation, $L$ is the length scale and $\bar{V}$ is the mean speed.

The von Kármán model is tuned by two parameters: the **length scale L** and the **turbulence intensity**.

**The length scale:** It represents the average size of a gust in appropriate directions and it is an important scaling factor in determining how rapidly gust properties vary in space. It depends on the surface roughness and the altitude. There are mathematical models to estimate the value of L. The most common is the one proposed by the Engineering Sciences Data Unit (ESDU) document [11]:

$$L_{ESDU} = 25z^{0.35}z_0^{-0.063} \qquad (3.3)$$

where $z$ is the altitude and $z_0$ is a surface roughness parameter. Typical values of $z_0$ are 0.01 for open sea and flat areas, 0.05 for open country, 0.3 for suburban and urban areas, and >1 for large city centers.

**The turbulence intensity:** It is expressed by:

$$I = \frac{\sigma}{\bar{V}}. \qquad (3.4)$$

It depends on the average wind speed but also on the altitude and the roughness of the surface as the length scale. Like the length scale, a mathematical
model exists to estimate wind turbulence intensity. The one proposed in the International Electrotechnical Commission (IEC) document [15] [22] is:

$$I_{IEC} = I_{ref}(0.75 + \frac{5.6}{\bar{V}}) \qquad (3.5)$$

where $I_{ref} = 0.16, 0.14$ and $0.12$ corresponding respectively to high, medium and low turbulence intensity.

The length scale and the turbulence intensity combination determines the shape of the PSD and thus the quality of the model fitting to the measured wind. The existing mathematical models to calculate these two parameters (ESDU, IEC standards, Danish Standard DS 472...) are simple approximations, so they are not always accurate. More complex models exist, however they often depend on other parameters that are difficult to measure. Therefore, accessing to the right combination for a good fit is often difficult.

### 3.3 Shaping filters

Approximating the PSD of wind turbulence is not sufficient for real-time simulation. For this purpose, filters capable of generating random wind speed must be designed. These filters are called shaping or modeling filters. They shape a random white noise at the input into a process that has the desired spectral characteristics. In our case, the spectral characteristics of wind turbulence are given by the approximated PSD. Thus, applying this PSD to a white noise thanks to shaping filters allows to have a random signal at the output with the same characteristics as a wind turbulence. The shaping filter with a frequency response $H(f)$ will transform a white noise $x[n]$ into a signal of interest $y[n]$, namely:

$$Y(f) = H(f)X(f) \qquad (3.6)$$

where $Y(f)$ and $X(f)$ respectively are the Fourier transforms (FT) of the output $y[n]$ and the input $x[n]$.

<figure><img src="image_5.png" /><figcaption>Fig. 4 Expression of the output PSD $S_{yy}(f)$ depending on the input PSD $S_{xx}(f)$ and the filter frequency response $H(f)$</figcaption></figure>

The white noise $x[n]$ at the input will have a constant PSD equal to its variance:

$$S_{xx}(f) = \sigma_x^2. \qquad (3.7)$$

The PSD of the output process $S_{yy}(f)$ is therefore expressed as follows:

$$S_{yy}(f) = S_{xx}(f)|H(f)|^2 = \sigma_x^2|H(f)|^2. \qquad (3.8)$$
Assuming that the input is a zero-mean and unitary variance white noise ($\sigma_x^2 = 1$), the output PSD can be expressed as:

$$S_{yy}(f) = |H(f)|^2 = H(jf)H^*(jf), \quad (3.9)$$

where * denotes the conjugate operator.

Therefore, the frequency response $H(f)$ of the shaping filter can be constructed from (3.9).

In the case of wind turbulence modeling, the desired PSD is given by the expression of the von Kármán model, and therefore can be used to obtain the transfer function of a shaping filter, namely:

$$S_{VK}(f) = \frac{4\sigma^2 \frac{L}{V}}{\left(1 + 70.8\left(\frac{L}{V}f\right)^2\right)^{5/6}} = \frac{K}{\left(1 + (\tau f)^2\right)^{5/6}}, \quad (3.10)$$

where $K = 4\sigma^2 \frac{L}{V}$ and $\tau = \sqrt{70.8} \frac{L}{V}$.

The shaping filter issued from the von Kármán PSD approximation can be expressed as:

$$H_{VK}(jf) = \frac{\sqrt{K}}{(1 + j\tau f)^{5/6}}. \quad (3.11)$$

This frequency response (3.11) verifies the equation (3.9):

$$\begin{aligned} H_{VK}(jf)H_{VK}^*(jf) &= \frac{\sqrt{K}}{(1 + j\tau f)^{5/6}} \frac{\sqrt{K}}{(1 - j\tau f)^{5/6}} \\ &= \frac{K}{(1 + (\tau f)^2)^{5/6}} \\ &= S_{VK}(f). \end{aligned} \quad (3.12)$$

The transfer function stemming from the von Kármán PSD model is consequently:

$$H_{VK}(s) = \frac{\sqrt{K}}{\left(1 + \frac{\tau}{2\pi} s\right)^{5/6}}, \quad (3.13)$$

where $s$ symbolizes the Laplace variable.

This transfer function can be used later in dynamical models of mechanical systems such as outdoor parabolic antennas to simulate wind-induced disturbances.

# 4 von Kármán model study

## 4.1 Model identification

The nonlinear simplex [25] is used to estimate the parameters of the von Kármán model allowing a good fit to the PSD of the available data set. The
nonlinear simplex method proposed in [25] and [23] is an extension of the classic simplex algorithm to handle nonlinear optimization and inequality constraints.

The expression of the von Kármán model (3.2) that will be estimated is:

$$ \hat{S}_{VK}(f) = \frac{K}{(1 + (\tau f)^2)^{5/6}}, \qquad (4.1) $$

where $K = 4\sigma^2 \frac{L}{V}$ and $\tau = \sqrt{70.8}\frac{L}{V}$.

The nonlinear simplex algorithm will estimate the parameters $K$ and $\tau$ by minimizing the cost function corresponding to the normalized sum of the square errors between the measured PSD, $S$, and the estimated model, $\hat{S}$, as:

$$ J = \frac{1}{N} \sum_{i=1}^{N} \left| S_{dB}(f_i) - \hat{S}_{dB}(f_i, K, \tau) \right|^2 \qquad (4.2) $$

where $N$ is the number of samples and $f_i$ are the frequencies defined in the chosen frequency range.

The optimization is carried out in dB in order to maintain a better accuracy in the considered frequency range. Indeed, the PSD modulus is tending towards zero, thus leading to suboptimal optimization.

The estimated von Kármán model obtained on Matlab® is presented in Figure 5. Figure 5 shows that the von Kármán model fits the measured wind

<figure><img src="image_6.png" /><figcaption>Fig. 5 Comparison between the estimated von Kármán PSD (blue) and the measured wind PSD (green)</figcaption></figure>

overall, especially in low frequencies where a good estimation of the gain is observed. However, for $f > 0.01Hz$, the slope of the model seems to be different from the one of the measurements. Therefore, the von Kármán model does not perform very well in middle and high frequency. However, according to the spectrum of Van der Hoven (Figure 3), the frequency range of wind turbulence
is between 0.01Hz and 0.1Hz which corresponds to the range where the model is less accurate. This can be explained by the fact that the von Kármán model was originally designed for aircraft, and consequently for high altitudes and moving systems. In our case, the data set is used under different conditions: low altitude and slow moving systems. Therefore, there is a need for more precisely wind turbulence models under these specific conditions.

## 4.2 Model parameters

The von Kármán model allows to well estimate the gain in low frequency, however the slope in high frequency is not sufficiently accurate. In this section, a particular attention is paid to the estimated parameters K and τ for different datasets. Figure 6 shows the estimated parameters K and τ.

<figure><img src="image_7.png" /><figcaption>Fig. 6 The estimated parameters K and τ of the von Kármán model for different datasets</figcaption></figure>

The values of the parameters depend on the mean speed and on the reference standard deviation σ which are assumed known and provided. Then, there is the length scale which is a more complicated input. It is supposed to be a constant since the measurements are made at the same place, but it is almost impossible to know its exact value because it depends on several factors, and at best it can be approximated by mathematical models (3.3). According to the von Kármán model, both model parameters K and τ depend on this length scale L and can be expressed as:

$$K = 4\sigma^2 \frac{L}{\bar{V}} \qquad (4.3)$$

and

$$\tau = \sqrt{70.8} \frac{L}{\bar{V}} = 8.4 \frac{L}{\bar{V}}. \qquad (4.4)$$
The approximated value of the length scale from the ESDU model is:

$$L_{ESDU} = 25z^{0.35}z_{0}^{-0.063} = 109.8m, \qquad (4.5)$$

where $z = 40m$ and $z_0 = 0.05$.

Extracting the length scale from (4.3) and (4.4) respectively from the estimated parameters K and $\tau$ gives:

$$L_K = \frac{K\bar{V}}{4\sigma^2} \qquad (4.6)$$

and

$$L_{\tau} = \frac{\tau \bar{V}}{19.5}. \qquad (4.7)$$

Figure 7 shows $L_K$ and $L_\tau$ the length scale value calculated respectively from the estimated parameters K and $\tau$. $L_K$ directly stems from the von Kármán model gain and therefore, the same coefficient 4 appears in the denominator. On the other hand, $L_\tau$, needs to be adjusted with a coefficient of 19.5 in order to be close to $L_K$.

<figure><img src="image_8.png" /><figcaption>Fig. 7 Comparison between $L_K$ the length scale value from K and $L_\tau$ obtained from $\tau$</figcaption></figure>

The obtained values of $L_K$ and $L_\tau$ vary between 80m and 130m (see Figure 7 for different datasets) with an average value of 106.4m for $L_K$ and 109.2m for $L_\tau$, therefore, the two parameters K and $\tau$ do not provide aberrant values but rather values quite close to the approximated $L_{ESDU}$. $L_K$ and $L_\tau$ have also the same evolution, which confirms that K and $\tau$ are related as in the true expression of the von Kármán model.

The estimated gain provides directly accurate values of the length scale, modifying the gain expression of the von Kármán model was not needed. This confirms that the von Kármán model well estimates the static gain as it was observed in Figure 5. Nevertheless, the coefficient $\sqrt{70.8}$ was too low
to provide coherent values of the length scale. This means that the estimated time constant $\tau$ was too high. This can be explained by the fact that the slope of the von Kármán model fixed by the non-integer order 5/6 does not fit that of the measured PSD and it seems lower. Therefore, the cutoff frequency should be lower to have a better fit.

## 5 Fractional model approximations

Fractional calculus has extended the possibilities of time and frequency identification by releasing the limit of the integer order. Two elementary fractional functions, the Cole-Cole [6] [24] [23] and Davidson-Cole [7] [24] [23] functions, are generally used to model complex systems with a reduced number of parameters. In fact, the von Kármán model which uses a non-integer order, produces a filter whose transfer function (3.13) is in the form of a Davidson-Cole function, with a differentiation order equal to 5/6. Therefore, these two types of fractional functions, Cole-Cole and Davidson-Cole, will be relied on to improve the accuracy of the von Kármán model over the frequency band of interest range of wind turbulence.

### 5.1 Davidson-Cole based model

The von Kármán model does not well fit the measured data in high frequency. Indeed, the slope obtained with von Kármán is fixed by the non integer order 5/6. However, according to the real data, this constraint can be released. Consequently, it is necessary to estimate the correct order of the model which may not be of 5/6. In order to improve the accuracy of the model, the fixed order will be replaced by an unknown parameter $\nu$ which will be estimated by the nonlinear simplex method, along with $K$ and $\tau$, thus there is an additional parameter to estimate compared to the von Kármán model.

The new model will be called Davidson-Cole (DC) model and its transfer function has the following expression:

$$H_{DC}(s) = \frac{\sqrt{K}}{\left(1 + \frac{\tau}{2\pi}s\right)^{\nu}} \qquad (5.1)$$

Using the equation (3.9), the PSD approximated by the Davidson-Cole model can be expressed as follows:

$$\hat{S}_{DC}(f) = \frac{K}{(1 + (\tau f)^2)^{\nu}} \qquad (5.2)$$

The Davidson-Cole model keeps the same general structure as the von Kármán model.

The estimated PSD from the identified Davidson-Cole (DC) model is presented in Figure 8.
<figure><img src="image_9.png" /><figcaption>**Fig. 8** Comparison between the estimated von Kármán model (blue), DC model (red) and the measured wind PSD</figcaption></figure>

Remark 1: According to Van der Hoven wind spectrum (see Figure 3), the wind turbulence component is characterized in the frequency range $[\frac{1}{10\text{min}}, \frac{1}{5\text{s}}] = [5.55 \times 10^{-3}\text{Hz}, 0.2\text{Hz}]$, meaning that the wind turbulence has more power inside this specific frequency range. In other words, the estimation at low frequencies ($f < 5.55 \times 10^{-3}\text{Hz}$) cannot be optimal and consequently, is not necessary.

Therefore, the DC model has difficulties to estimate the static gain in low frequency.

It can be seen that the high frequency slope of the DC model well corresponds to the measurements and therefore meets the expectations. However, the gain in low frequency seems to be a little bit low in comparison with the measured data PSD and the von Kármán model.

The main difference between the DC model and the von Kármán model is the differentiation order. Figure 9 shows the estimated values of $\nu$ for different datasets.

<figure><img src="image_10.png" /><figcaption>**Fig. 9** The estimated parameter $\nu$ of the Davidson-Cole model</figcaption></figure>

The estimated values of $\nu$ varies between 1.30 and 1.48, so it is very different from 5/6 $\approx 0.83$ which is the fixed order of the von Kármán model.
This confirms that the von Kármán model is not adapted to model the wind turbulence in our conditions and especially in high frequency.

In order to compare the quality of the estimated models, the Akaike information criterion (AIC) [1] will be used. The AIC estimates the error between the model and the data but it includes a term penalizing the number of parameters. It deals with the trade-off between the goodness of fit of the model and the simplicity of the model. The normalized AIC (nAIC) which will be used has the following expression:

$$nAIC = N \log(\det(\frac{1}{N}E^tE)) + \frac{2n_p}{N}, \qquad (5.3)$$

where N is the number of samples, E is the vector of errors expressed as:

$$E_i = S_{dB}(f_i) - \hat{S}_{dB}(f_i, K, \tau), \quad \forall i = 1, \dots, N, \qquad (5.4)$$

and $n_p$ is the number of parameters.

<figure><img src="image_11.png" /><figcaption>Fig. 10 Comparison between the nAICs of the identified models: von Kármán (blue) and DC (red) for 28 different data sets</figcaption></figure>

Figure 10 shows that the Davidson-Cole model has a lower AIC than the von Kármán model for all of the datasets. This result confirms that a DC model is more suited and that the differentiation order is other than 5/6 as provided by the von Kármán model.

## 5.2 Cole-Cole based model

In order to generate random wind speed for time simulation, filters with fractional transfer functions need to be approximated by rational functions. Cole-Cole functions are simpler to synthesis than Davidson-Cole functions. Therefore, a wind turbulence model can be constructed from Cole-Cole functions
and compared to the von Kármán and the DC models. The Cole-Cole model (CC) has the following transfer function:

$$H_{CC}(s) = \frac{\sqrt{K}}{1 + \left(\frac{\tau}{2\pi} s\right)^{\nu}}. \qquad (5.5)$$

Its resulting PSD will be used for parameter identification to fit the spectral characteristics of the wind turbulence. The PSD expression is deduced from $H_{CC}(s)$ using equation (3.9):

$$ \begin{aligned} \hat{S}_{CC}(f) &= H_{CC}(jf)H_{CC}(jf)^* \\ &= \frac{\sqrt{K}}{1 + (j\tau f)^\nu} \frac{\sqrt{K}}{1 + (-j\tau f)^\nu} \qquad (5.6) \\ &= \frac{K}{(1 + (j^\nu + (-j)^\nu)(\tau f)^\nu + (\tau f)^{2\nu})} \end{aligned} $$

The expression of the Cole-Cole model PSD, $\hat{S}_{CC}(f)$, is then:

$$ \hat{S}_{CC}(f) = \frac{K}{1 + 2 \cos\left(\nu \frac{\pi}{2}\right) (\tau f)^\nu + (\tau f)^{2\nu}}. \qquad (5.7) $$

There are 3 parameters to estimate $K$, $\tau$ and $\nu$ like the DC model, but the structure is different. Figure 11 presents the PSD obtained from the CC model compared to the previous models.

<figure><img src="image_12.png" /><figcaption>Fig. 11 Comparison between the estimated von Kármán model (blue), DC model (red), CC model (pink) and the measured wind PSD</figcaption></figure>

Like the DC model, the CC model provides a good fit in high frequency (see Remark 1 p.13). The slope of the model follows well the PSD measured. However, the gain in low frequency seems to be underestimated. It is even
lower with the DC model. Therefore, like the DC model, the CC model has good performance in high frequency but a poor approximation in middle and low frequencies.

<figure><img src="image_13.png" /><figcaption>Fig. 12 Comparison between the nAICs of the identified models: von Kármán (blue), DC (red) and CC (pink) for 28 different data sets</figcaption></figure>

Moreover, the CC model has a higher AIC than the DC model (see Figure 12). This difference in performance is reflected in the gain at low frequencies, since at high frequencies their accuracy seems similar. Consequently, with the CC model, filters can be synthesized more directly but there is a loss in quality for the same number of parameters.

## 5.3 Improved Cole-Cole based model

The Davidson-Cole and Cole-Cole models well estimate the high frequency slope, however they do not well fit in low and middle frequencies. They underestimate the static gain and thus the amplitude of the turbulence around the mean wind speed and as a result, they reduce their impact on the system dynamics. Therefore, there is a need to improve the low and middle frequency estimation, and thus the proposed fractional models have to be modified. The improvement of the CC model is preferred, to facilitate the approximation to rational functions later.

To improve the model, a second cell of a fractional order $2\nu$ is added to the CC fractional function of the PSD, which will be denoted as Cole-Cole×2 (CC×2):

$$H_{CC \times 2}(s) = \frac{\sqrt{K}}{\left(1 + \left(\frac{\tau_1}{2\pi} s\right)^\nu\right)\left(1 + \left(\frac{\tau_2}{2\pi} s\right)^{2\nu}\right)}. \quad (5.8)$$
Based on the PSD calculation (5.6)-(5.7) for the CC model, the PSD obtained from the CC×2 model is expressed as follows:

$$ \hat{S}_{CC \times 2}(f) = \frac{K}{D_1(f)D_2(f)} \qquad (5.9) $$

where:

$$ D_1(f) = 1 + 2 \cos\left(\nu\frac{\pi}{2}\right)(\tau_1 f)^{\nu} + (\tau_1 f)^{2\nu} \qquad (5.10) $$

and

$$ D_2(f) = 1 + 2 \cos(\nu\pi)(\tau_2 f)^{2\nu} + (\tau_2 f)^{4\nu}. \qquad (5.11) $$

Remark 2: The fractional order on $D_2$ could have been independent of the one of $D_1$, $\nu_1$. In this case, an additional parameter $\nu_2$ would be added to the optimization parameter vector. After several trials, it has been noticed that this fractional order is very close to $2\nu_1$. Therefore, it has been chosen to propose only one model such a given in relations (5.9)-(5.11).

The interest of adding a second cell of fractional order is to bring a second slope which will allow the model to have a better estimation in the low and middle frequency ranges. Nevertheless, the model has one more parameter to estimate and its PSD structure is more complex. Figure 13 presents the PSD obtained from the CC×2 model compared to the previous proposed fractional models.

<figure><img src="image_14.png" /><figcaption>Fig. 13 Comparison between the estimated von Kármán model (blue plain round), DC model (red dash-dotted diamond), CC model (pink dashed round), CC×2 model (purple square) and the measured wind PSD (green plain)</figcaption></figure>

The second slope brought by the second cell is clearly visible for $f < 0.01Hz$. Indeed in this frequency range, the PSD of the model CC×2 does not have a constant gain as with the DC and CC models, but a second slope which provides a better fit as desired.
To confirm the quality of the CC×2 model observed in Figure 13, the nAIC of the different models presented in this study will be compared (see Figure 14).

<figure><img src="image_15.png" /><figcaption>**Fig. 14** Comparison between the nAICs of the identified models: von Kármán (blue plain round), DC model (red dash-dotted diamond), CC model (pink dashed round) and CCx2 model (purple square) for 15 different data sets</figcaption></figure>

Figure 14 confirms that the CC×2 model fits the best the wind data, since its AIC is the lowest. Table 1 summarizes the results observed in figures 10 and 14.

<table><caption>**Table 1** Comparison between the mean nAIC of the different models</caption><thead><tr><td>Model</td><td>Number of parameters</td><td>Mean nAIC</td></tr></thead><tbody><tr><td>von Kármán</td><td>2</td><td>2.027</td></tr><tr><td>Davidson-Cole</td><td>3</td><td>1.200</td></tr><tr><td>Cole-Cole</td><td>3</td><td>1.283</td></tr><tr><td>Cole-Cole×2</td><td>4</td><td>1.162</td></tr></tbody></table>

The CC×2 model fits the wind measurements well over the entire turbulence frequency range. Even if it is a little more complex than the other models, it is the one chosen for the time simulation and to be integrated into the dynamic model of our system.
# 6 Wind model design and time simulation

## 6.1 Model parameter tuning

In this section, a closer look is brought for the estimated parameters K, v, τ₁ and τ₂ of the CC×2 model, and then a method to tune them is proposed if there is no access to reference data. The estimated values of the CC×2 are presented in Figure 15.

<figure><img src="image_16.png" /><figcaption>Fig. 15 The estimated parameters K, v, τ₁ and τ₂ of the Cole-Cole×2 model</figcaption></figure>

As in the study of the von Kármán model in Section 4.2, the length scale will be the link between the static gain K and the time constants τ₁ and τ₂. As both models, CC×2 and von Kármán ones, provide a good fit at low frequencies with a well approximated static gain, the CC×2 model gain K is assumed to have the same expression as the von Kármán gain:

$$K = \frac{4\sigma^2L}{\bar{V}}. \qquad (6.1)$$

Then, Lₖ the value of the length scale is deduced from the estimated parameter K:

$$L_k = \frac{K\bar{V}}{4\sigma^2}. \qquad (6.2)$$
The von Kármán model has a single time constant $\tau$ from which $L_{\tau}$ is deduced, while the CC×2 model has two time constants $\tau_1$ and $\tau_2$.

$$L_{\tau} = \frac{\overline{V\tau}}{19.5}. \qquad (6.3)$$

Therefore, the von Kármán model provides the inspiration to find a relationship between $\tau_1$, $\tau_2$ and the length scale that will be called $L_{12}$:

$$L_{12} = \frac{\overline{V\sqrt{\tau_1\tau_2}}}{4.7}. \qquad (6.4)$$

To verify these two relationships, the values of $L_K$ and $L_{12}$ are compared for different datasets to the approximated value of the length scale $L_{ESDU} = 109.8m$.

<figure><img src="image_17.png" /><figcaption>**Fig. 16** Comparison between $L_K$ the length scale value from (6.2) and $L_{12}$ obtained from (6.4)</figcaption></figure>

The obtained values of $L_K$ and $L_{12}$ vary between 90m and 140m, therefore, the estimated model parameters do not provide aberrant values of length scale but rather values quite close to the approximated $L_{ESDU}$. It can be noticed that $L_K$ and $L_{12}$ have the same evolution, which confirms that K and time constants are related as in the von Kármán model.

It is also necessary to find a simple relationship between the time constants. In this way, $\tau_1$ can be calculated using the mean speed and the length scale and then $\tau_2$ can be deduced from $\tau_1$. The ratio between the two quantities is between 2.6 and 4.6 with an average of 3.6 (Figure 17). Therefore, $\tau_2$ can be approximated by:

$$\tau_2 = \frac{\tau_1}{3.6}. \qquad (6.5)$$
This is an optimistic approximation but it allows to have consistent values for the model in the absence of a reference dataset. The equation to calculate $\tau_1$ from the inputs becomes:

$$ \tau_1 = 8.9 \frac{L}{V}. \qquad (6.6) $$

<figure><img src="image_18.png" /><figcaption>Fig. 17 Ratio between $\tau_1$ and $\tau_2$</figcaption></figure>

Therefore, the three equations (6.1), (6.5) and (6.6) allow to calculate the parameters of the CC×2 model from the mean speed, the standard deviation and the length scale. For the order $\nu$, the mean value: $\nu = 0.516$ is chosen.

## 6.2 Rational filter design

Fractional operators are usually approximated by high-order rational functions in a limited frequency band. The fractional model and its rational approximation must have the same dynamical behavior in this limited frequency band. The most commonly used approximation of the fractional integro-differentiator $s^{\nu}$ (with $\nu$ a non-integer order) in the limited frequency band $[\omega_A, \omega_B]$ is the recursive distribution of zeros and poles or the Oustaloup's approximation [20, 21]. However, the rational approximation of the fractional operators in $[\omega_A, \omega_B]$ induces deterioration around the limit frequencies $\omega_A$ and $\omega_B$. This deterioration is usually reduced by extending the frequency band over which the approximation is performed from $[\omega_A, \omega_B]$ to $[\omega_B, \omega_h]$. In the case of wind turbulence, the frequency band $[\omega_A, \omega_B]$ will be $2\pi[10^{-4}, 10^{-1}]$. Therefore, the rational approximation will be carried out on $2\pi[10^{-6}, 10^1]$. $s^{\nu}$ is approximated by:

$$ s^{\nu} \rightarrow s^{\nu}_{[\omega_A, \omega_B]} = C_0 \left( \frac{1 + \frac{s}{\omega_b}}{1 + \frac{s}{\omega_h}} \right)^{\nu} \approx C_0 \prod_{k=1}^{N} \frac{1 + \frac{s}{\omega'_k}}{1 + \frac{s}{\omega_k}} \qquad (6.7) $$

where:

$$ \omega'_1 = \eta^{1/2}\omega_b, \quad \omega_1 = \alpha\omega'_1, \qquad (6.8) $$

$$ \omega'_{k+1} = \alpha\eta\omega'_k, \quad \omega_{k+1} = \alpha\eta\omega_k, \qquad (6.9) $$
$$ \nu = \frac{\log \alpha}{\log(\alpha\eta)}, \quad \alpha = \left(\frac{\omega_h}{\omega_b}\right)^{\frac{\nu}{N}}, \quad \eta = \left(\frac{\omega_h}{\omega_b}\right)^{\frac{1-\nu}{N}} \qquad (6.10) $$

and

$$ C_0 = \left| \frac{1 + j \frac{1}{\omega_h}}{1 + j \frac{1}{\omega_b}} \right|^\nu = \left( \frac{\omega_h}{\omega_b} \right)^\nu \left( \frac{1 + \omega_b^2}{1 + \omega_h^2} \right)^{\nu/2} . \qquad (6.11) $$

<figure><img src="image_19.png" /><figcaption>Fig. 18 Fractional integrator $s^{-0.5}$ and its approximation $s_{[0.1,10]}^{-0.5}$ based on Oustaloup's recursive poles and zeros approximation</figcaption></figure>

Remark 3: Time-domain simulation of fractional systems is an extensively studied topic in the literature. Note that any time-domain simulation algorithm of fractional systems can be used, at the discretion of the user. However, fractional systems should be correctly simulated with negligible simulation errors.

## 6.3 Time simulation

In this section, a time simulation is carried out from the CC×2 model. The reference wind data are those shown in Figure 1. The CC×2 model identified directly from the reference is a black box model, contrary to from a grey box CC×2 model tuned with the relations (6.1), (6.5) and (6.6) from Section 6.1. Their simulation results will be compared in order to check the validity of the grey box model. The reference signal has a mean speed $\bar{V} = 6.6$ m/s and a standard deviation $\sigma = 1.92$. The estimation of the CC×2 parameters using the nonlinear simplex method gives the following values:

$$ K = 301.09, \quad \tau_1 = 179.17s, \quad \tau_2 = 50.13s, \quad \nu = 0.518. \qquad (6.12) $$
If the parameters of the model are calculated from the known inputs $\bar{V} = 6.6\,\text{m/s}$, $\sigma = 1.92$ and with $L$ and $\nu$ fixed at the mean value of the estimations $L = 120\,\text{m}$ and $\nu = 0.516$, the model will have the following parameter values:

$$K = 269.94, \quad \tau_1 = 161.74\,\text{m}, \quad \tau_2 = 44.93\,\text{s}. \qquad (6.13)$$

Figure 19 shows the PSD obtained from the black box identified and the tuned grey box CC×2 models.

<figure><img src="image_20.png" /><figcaption>Fig. 19 Comparison between the measured wind PSD (green), the black box CC×2 (purple square) and the grey box CC×2 (red)</figcaption></figure>

It can be seen that the grey box model is very close to the one identified and more importantly, it still well fits the measured wind PSD. This shows that the approximations chosen to tune the model provide fairly consistent values.

The fractional transfer functions of the shaping filters can be calculated from the model parameters. The black box model has the following transfer function:

$$H_{BlackBox}(s) = \frac{17.35}{48.72s^{1.554} + 8.59s^{1.036} + 5.67s^{0.518} + 1}. \qquad (6.14)$$

For the grey box model, the transfer function is:

$$H_{GreyBox}(s) = \frac{16.43}{40.7s^{1.548} + 7.62s^{1.032} + 5.34s^{0.516} + 1}. \qquad (6.15)$$

Rational transfer functions are then obtained from these two fractional functions using the Oustaloup's approximation. Finally, these rational functions allow to make real-time simulations by putting as an input a random white noise with a zero mean, a standard deviation $\sigma_{noise} = 1$, a sampling period of 1s and a duration of 3h. This white noise is then used as input for the filters.
The filters then should return signals with the same spectral characteristics as the wind turbulence. These signals have a zero mean, so the desired wind speed need to be added at the end. Figure 20 shows the time domain outputs of both models, identified and tuned, and the error between them. The mean value $\bar{V}$ = 6.6m/s, were added to the signals. Then, to check the accuracy of the signals and whether they model wind turbulence, their spectral characteristics and therefore their PSDs are compared with that of the reference signal (see Figure 21).

<figure><img src="image_21.png" /><figcaption>**Fig. 20** Comparison between the output signal from the black box CC×2 (purple square) and the grey box one (red) and the error between them</figcaption></figure>

The PSDs of the signals generated by the wind models from a random white noise have the same shape as the PSD of the wind measurements. This confirms that the signals of the Figure 20 well reproduce the evolution of wind speed, and consequently they can be used in a dynamical model to simulate the disturbances caused by the wind. It can be noticed also that the grey box model output signal is very similar to the one generated by the directly identified model. The error between them has a maximum value of 0.4m/s. This confirms the validity of the grey box model.
<figure><img src="image_22.png" /><figcaption>Fig. 21 Comparison between the measured wind PSD (green), the black box CC×2 generated PSD (blue) and the grey box CC×2 PSD (red)</figcaption></figure>

## 6.4 Comparison with rational models

In this section, it is proposed to compare our improved Cole-Cole CC×2 model with rational models. Let us consider the normalized Akaike information criterion (nAIC), such as described in Section 5.1, and the fitness criterion (Fit) defined in percentage as:

$$ \text{Fit} = \left( 1 - \frac{\|S_{dB} - \hat{S}_{dB}\|^2}{\|S_{dB}(f_i)\|^2} \right) \times 100, \qquad (6.16) $$

where $S_{db}$ denotes the vector of the measured PSD and $\hat{S}_{dB}$ the vector of the PSD stemming from the approximated model. Table 2 summarizes the comparative results of the proposed rational approximations, that can be observed in Figure 22.

It can clearly be seen that the CC×2 model which has only 4 parameters, provides the best results as the Fit criterion is the closest to 100%, and the nAIC criterion is the lowest.

If two approximations are carried out by using the Oustaloup approximation, such as described in Section 6.2, the resulting approximations have the following form :

$$ H_{N_{\text{cells}}}^{\text{Oustaloup approx}} = H_0 \prod_{k=1}^{N} \frac{1 + \frac{s}{\omega_k'}}{1 + \frac{s}{\omega_k}}, \qquad (6.17) $$

where $N_{cells}$ denotes the number of cells used per decade. Indeed, the higher $N_{cells}$, the better the approximation. Note that the frequency band is set to $[10^{-6} \text{Hz}, 1 \text{Hz}]$, therefore $N = 7$. For $N_{cells} = 2$, the corresponding approximation model, $H_2^{\text{Oustaloup approx}}$, provides a Fit criterion of 91.51% and a nAIC which is almost twice as the one obtained with the CC×2. To achieve the
same performances as the $CC \times 2$ fractional model, the corresponding approximation model, $H_4^{Oust approx}$ should have $N_{cells} = 4$ thus leading to a rational transfer function that has 26 parameters.

If the user would like to use a rational model with only 4 parameters, such as the $CC \times 2$ model, the optimization with the same nonlinear simplex method provides:

$$H_4^{Rat approx} = H_0 \frac{1 + \frac{s}{\omega_2}}{\left(1 + \frac{s}{\omega_1}\right) \left(1 + \frac{s}{\omega_3}\right)}, \quad (6.18)$$

and the resulting performances of such a model are really poor as its Fit criterion equals 73.17% and its nAIC reaches almost thrice as the one obtained with the $CC \times 2$.

<table><caption>**Table 2** Comparison with classic rational models</caption><thead><tr><th>Model</th><th>Number of parameters</th><th>Fit (%)</th><th>nAIC</th></tr></thead><tbody><tr><td>Cole-Cole×2</td><td>4</td><td>96.473</td><td>1.193</td></tr><tr><td>Oust. Approx $H_2^{Oust approx}$</td><td>14</td><td>91.514</td><td>2.083</td></tr><tr><td>Oust Approx $H_4^{Oust approx}$</td><td>26</td><td>96.502</td><td>1.214</td></tr><tr><td>Rat. Approx $H_4^{Rat approx}$</td><td>4</td><td>73.166</td><td>3.222</td></tr></tbody></table>

<figure><img src="image_23.png" /><figcaption>**Fig. 22** Comparison between the measured wind PSD (green), the black box $CC \times 2$ generated PSD (blue), and rational models</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Fig. 23 Scheme of the wind speed generator</figcaption></figure>

## 6.5 Wind speed generator

The real-time wind speed simulator is built from the shaping filter, stemmed from the wind model, and from the different equations to tune its parameters. The scheme in Figure 23 illustrates the functioning of the wind speed generator. The model parameters $K$, $\tau_1$ and $\tau_2$ are calculated from the mean speed $\bar{V}$, the length scale $L$ and the turbulence intensity $I_{ref}$ which are constants and chosen by the user. The order $\nu$ is also fixed. The generator input is a white noise with a standard deviation $\sigma_{noise} = 1$. The filter output is the short-term component of the wind speed $V_{turb}(t)$ representing the turbulence speed. The constant mean speed $\bar{V}$ is then added to obtain the wind speed. This signal can be injected later into a system dynamical model to simulate disturbance.

The Matlab® script showing how to compute the CC×2 model transfer function for the wind speed generator is provided in 8. All simulations have been carried out with the CRONE toolbox developed in Matlab®, which is dedicated to fractional calculus, fractional system simulation and system identification with fractional models. The new version of the CRONE toolbox is an object-oriented version with several classes defined for fractional models (LTI, explicit form, implicit form, ZPK, state-space representation, etc.). This CRONE toolbox is freely available at http://archive.ims-bordeaux.fr/CRONE/toolbox/.

As an example, by providing the following parameter values to the wind simulator: $\bar{V}$ = 10m/s, $I_{ref}$ = 0.12, $L$ = 110m and $\nu$ = 0.51, Figure 24 shows the generated time signal of the wind turbulence.

Figure 25 illustrates the generated PSD and the PSD corresponding to the tuned CC×2 model. This example shows that our generator is able to accurately reproduce wind turbulence under low altitude conditions.

## 7 Conclusion

Wind speed can be decomposed into a long-term component representing the mean speed and a short-term component corresponding to turbulence. The mean speed is considered constant and therefore it is the turbulence speed that
<figure><img src="image_25.png" /><figcaption>**Fig. 24** An example of a time domain signal generated by the wind speed generator</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 25** The PSD of the Figure 24 wind speed signal</figcaption></figure>

should be modeled. Wind spectral characteristics are relied on to model wind fast fluctuations because turbulence is a stochastic process and it is difficult to model it in time domain. The von Kármán model is the reference in terms of spectral modeling of wind turbulence. It is widely used and it remains very effective, especially for aircraft. It is accurate in low frequency but its fixed differentiation order prevents optimal modeling in high and middle frequencies. Therefore, improved models based on the von Kármán model structure have been developed and tested. As the von Kármán model uses fractional calculus, the new models have the Davidson-Cole and Cole-Cole forms of fractional functions. They allow a better fit in high frequency, because they have no constraints on the differentiation order. However, the accuracy of the fit is reduced at low and middle frequencies, and the static gain is underestimated.
Consequently, another fractional model has been introduced. It is based on a Cole-Cole two-cell fractional function, enabling the PSD to have a second slope. This model is the one that fits the best but has one more parameter to estimate. A method for tuning its parameters has been proposed for the case where there is no access to reference data. The CC×2 model was then used to synthesize a fractional transfer function of a shaping filter capable of generating a random turbulence speed. The time simulation results showed that the signals generated by the model have the same spectral characteristics as the measured wind. A wind speed generator, based on a grey box CC×2 model, was proposed. The model parameters are tuned from known constants: the mean speed, the length scale and the turbulence intensity. The generator takes as an input a random white noise and returns a signal with the same spectral characteristics as wind turbulence. The real-time wind simulator will be integrated into dynamical models of tracking parabolic antennas, in order to develop control laws capable of countering the wind effects effectively and thus making the antennas more precise. In this study, the wind mean speed was considered constant because only short-term fluctuations were studied, however the mean speed varies slowly depending on the atmospheric conditions and the season. It could be interesting to integrate a long period wind mean speed generator into the real-time wind generator.

## Conflict of interest

The authors have no relevant financial or non-financial interests to disclose. The authors declare that they have no conflict of interest.

# 8 Appendix: Matlab script for the wind turbulence generator

Figure 26 provides the Matlab script in order to get the Cole-Cole×2 transfer function approximation. Such a rational filter can then be used for the generation of the wind turbulence such as illustrated in Figure 23.

## References

1. Akaike, H.: A new look at the statistical model identification. *IEEE Transactions on Automatic Control* **19**(6), 716-723 (1974). DOI:10.1109/TAC.1974.1100705

2. Banerjee, A., Chakraborty, T., Matsagar, V., and Achmus, M.: Dynamic analysis of an offshore wind turbine under random wind and wave excitation with soil-structure interaction and blade tower coupling. *Soil Dynamics and Earthquake Engineering* **125**, 105699 (2019). DOI:10.1016/j.soildyn.2019.05.038

3. Burton, T., Jenkins, N., Sharpe, D., and Bossanyi, E.: Wind Energy Handbook. John Wiley & Sons (2011)

4. Cadalen, B., Fabien, G., Lanusse, P., Sabatier, J., and Parlier, Y.: Modelling and control of a tethered kite in dynamic flight. *Journal of Sailing Technology* (2018). DOI:10.5957.jst.2018.03
1.% Model inputs
2. L = 110; %Length scale (m)
3. Vm = 80/3.6; %Mean speed (m/s)
4. I_ref = 0.12; %Turbulence intensity reference (low: 0.12, mean:0.14, % high: 0.16)
5. Ts = 1;
7. %% Tuning model parameters
8. I = I_ref.* (0.75+5.6/Vm);
9. sigma = I*Vm;
10. K = 4.* (sigma)^2.*L./Vm;
11. Tau1 = 8.*L./Vm;
12. Tau2 = Tau1./3;
13. nu = 0.51;
14.
15. %% Fractional transfer function
16. num=sqrt(K);
17. den1=frac.poly_exp([ (Tau1/2/pi)^nu 1],[nu 0]);
18. den2=frac.poly_exp([ (Tau2/2/pi)^ (2*nu) 1],[2*nu 0]);
19.
20. sys = frac.tf(num,den1)*frac.tf(1,den2);
21
22. %% Approximated rational transfer function
23. set(sys,'N',10); %Cell number
24. set(sys,'band', [1e-4,10]); %Frequency range + 2 decades
25. set(sys,'sim','Oustaloup');
26
27. H_CCx2 = frac2int(sys);

<caption>Fig. 26 Matlab script showing how to compute the Cole-Cole×2 transfer function</caption>

5. Cadalen, B., Lanusse, P., Sabatier, J., Fabien, G., and Parlier, Y.: Robust control of a tethered kite for ship propulsion. In European Control Conference 2018, Limassol, Cyprus, 3065–3070 (2018). DOI:10.23919/ECC.2018.8550274

6. Cole, K. S., Cole, R. H.: Dispersion and absorption in dielectrics I. Alternating current characteristics. Journal of Chemical Physics 9(4), 341–351 (1941). DOI:10.1063/1.1750906

7. Davidson, D. W. and Cole, R. H.: Dielectric relaxation in glycerol, propylene glycol, and N-propanol. Journal of Chemical Physics, 19(12), 1484–1490 (1951). DOI:10.1063/1.1748105

8. de Kármán, T. and Howarth, L.: On the statistical theory of isotropic turbulence. Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences 164(917), 192–215 (1938). DOI:10.1098/rspa.1938.0013

9. Diop, A. D., Ceanga, E., Rétiveau, J.-L., Méthot, J.-F., and Ilinca, A.: Real-time three-dimensional wind simulation for windmill rig tests. Renewable Energy 32(13), 2268–2290 (2007). DOI:10.1016/j.renene.2006.04.011

10. Duhé, J.-F., Victor, S., Melchior, P., Abdelmounen, Y., and Roubertie, F.: Modeling thermal systems with fractional models: human bronchus application. Nonlinear Dynamics 108, 579–595 (2022). DOI:10.1007/s11071-022-07239-3

11. ESDU: 85020-Characteristics of atmospheric turbulence near the ground. Part II: single point data for strong winds (neutral atmosphere). Technical report, ESDU (1985)

12. Feytout, B.: Commande CRONE appliquée à l'optimisation de la production d'une éolienne. PhD Thesis, Université Bordeaux I (2013)

13. Hajjem, M., Victor, S., Melchior, P., Lanusse, P., and Thomas, L.: Wind fractional modeling by spectral analysis. In 2021 9th International Conference on Systems and Control (ICSC), 436–441 (2021). DOI:10.1109/ICSC50472.2021.9666651

14. Hajjem, M., and Lanusse, P., and Victor, S., and Melchior, P., and Thomas, L.: Robust control of 2-axis tracking antennas using the CRONE approach. In 1st IFAC Workshop on Control of Complex Systems COSY 2022, Vol. 55(40), 169–174, IFAC PapersOnLine (2022). DOI:10.1016/j.ifacol.2023.01.067
15. IEC: Wind turbines – Part 1: Design requirements. 61400-1. Standard, International Electrotechnical Commission (2005)

16. Lei, J., Huang, P., Zhang, L., Yuan, Y., Deng, W., Mao, S., and Zhang, J.: Experimental study on flow characteristics of a large-scale open jet wind tunnel for outdoor pool fire research. Journal of Wind Engineering and Industrial Aerodynamics **210**, 104522 (2021). DOI:10.1016/j.jweia.2021.104522

17. Leithead, W. E., De la Salle, S., and Reardon, D.: Role and objectives of control for wind turbines. In IEE Proceedings C-Generation, Transmission and Distribution, volume 138, pages 135–148. IET (1991). DOI:10.1049/ip-c.1991.0017

18. Mahgoub, Y. and El-Badawy, A.: Sliding mode control of directly excited structural dynamic model of twin-tailed fighter aircraft. Journal of the Franklin Institute **358**(18), 9721–9740 (2021). DOI:10.1016/j.jfranklin.2021.10.017

19. Moorhouse, D. J. and Woodcock, R. J.: Background information and user guide for MIL-F-8785c, military specification-flying qualities of piloted airplanes. Technical report, Air Force Wright Aeronautical Labs Wright-Patterson AFB OH (1982)

20. Oustaloup, A. : La dérivation non entière. Hermès (1995)

21. Oustaloup, A., Lebron, F., Mathieu, B., and Nanot, F. M.: Frequency-band complex noninteger differentiator: characterization and synthesis. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications **47**(1), 25–39 (2000). DOI:10.1109/81.817385

22. Ren, G., Liu, J., Wan, J., Li, F., Guo, Y., and Yu, D.: The analysis of turbulence intensity based on wind speed data in onshore wind farms. Renewable Energy **123**, 756–766 (2018). DOI:10.1016/j.renene.2018.02.080

23. Sommacal, L.: Synthèse de la fonction d'Havriliak-Negami pour l'identification temporelle par modèle non entier et modélisation du système musculaire. PhD Thesis, Université Bordeaux 1 (2007)

24. Sommacal, L., Melchior, P., and Oustaloup, A.: Havriliak-Negami function for thermal system identification. In Proceedings of the American Control Conference, pages 1316–1321 (2008). DOI:10.1109/ACC.2008.4586675

25. Subrahmanyam, M. B.: An extension of the simplex method to constrained nonlinear optimization. Journal of Optimization Theory and Applications, **62**(2), 311–319 (1989). DOI:10.1007/BF00941060

26. van der Hoven, I. : Power spectrum of horizontal wind speed in the frequency range from 0.0007 to 900 cycles per hour. Journal of Atmospheric Sciences **14**(2), 160–164 (1957). DOI:10.1175/1520-0469(1957)014<0160:PSOHWS>2.0.CO;2

27. Victor, S., Melchior, P., Pellet, M., and Oustaloup, A.: Lung thermal transfer system identification with fractional models. IEEE Transactions on Control Systems Technology **28**(1), 172–182 (2020). DOI:10.1109/TCST.2018.2877606

28. von Kármán, T.: Progress in the Statistical Theory of Turbulence. Proceedings of the National Academy of Sciences **34**(11), 530-9 (1948). DOI:10.1073/pnas.34.11.530

29. Welfonder, E., Neifer, R., and Spanner, M.: Development and experimental identification of dynamic models for wind turbines. Control Engineering Practice **5**(1), 63–73 (1997). DOI:10.1016/S0967-0661(96)00208-0