Controller Development and Experimental Validation for a Vertical Axis Wind Turbine

: This paper reports experimental results, in a wind tunnel, of controllers for a vertical axis wind turbine (VAWT). Models and controllers from previous works were validated and improved through the acquisition of new data sets. The control methods used are PID, LQR and CRONE (a type of fractional order control, i.e., a controller using fractional derivatives). The simulation results led to the choice of the CRONE controller for implementation. Hardware and software were developed to perform tests in a wind tunnel, that proved the correct behavior of the controller whilst validation results exceeded those from simulation. These results are expected to contribute for the viability of VAWTs in urban environments since they can aid in the increasing demand of energy services in urban areas without increasing the carbon emission levels of the area.


Introduction
Due to concerns related to the environment, the interest in renewable energy sources has been steadily increasing in order to successfully transition the energy production to become more environmentally friendly, with much lower carbon emissions. Wind energy is one of these sources, being one of the highest average growth rate sectors, second only to solar power, among countries in the Organization for Economic Co-operation and Development (OECD) [1].
While the concept of a wind turbine has existed since the 19th century, its technology is under constant development as manufacturing and control methods are improved. At the same time, new problems in its implementation are being discovered, as unintended or unforeseen consequences are still being uncovered. A wind turbine's advantages include a renewable, sustainable and never ending resource, and the ability to provide power to remote locations where other sources of energy might be impractical or too expensive to reach or exploit. On the other hand, the wind itself fluctuates in intensity; thus, the power production is not constant, leading to the necessity of energy storage in case of energy shortages. Furthermore, due to its constant rotation, the blades might need constant maintenance and eventual replacement due to erosion. Another disadvantage of the commonly used wind turbines nowadays is their size, which in turn means its components are very large, requiring specialized transportation and complex logistics for installation and component replacement.
These wind turbines possess a horizontal rotational axis, thus they will be referred to as horizontal axis wind turbines (HAWT) from this point onward. Due to their size, wind farm installation near an urban area is impractical, thus a vertical axis wind turbine (VAWT) may be the solution for electricity production in such sites, closer to where energy demand is larger. These wind turbines, VAWTs, possess some advantages over the traditional HAWTs [2] such as: • Their natural independence of the wind flow direction translates into a lower number of components; • They possess a smaller size, thus requiring lower logistical effort in their installation and easing the required maintenance.
Yet VAWTs are not void of disadvantages such as producing less power when compared to the traditionally used HAWT, and they inherit some of HAWTs' disadvantages such as the blades' erosion problem. Nevertheless, interest in and development of VAWTrelated technology has been increasing due to the previously mentioned advantages over HAWTs.

State of the Art
A few studies were conducted to analyze the performance of VAWT installations in an urban context. In [3], a more traditional approach was performed, resorting to data acquisition through anemometers and a wind tunnel, while in [4] the data was acquired through a computational fluid dynamics (CFD) calculation of flow velocity and skew angle on the rooftop of several buildings. In [5], several configurations of a Darrieus-type VAWT are described, along with some of their characteristics, drawbacks and advantages.
VAWTs can be mainly categorized as two types: Darrieus or Savonius. While a Darrieus-type VAWT uses an angle of attack on its blade to spin, a Savonius type is simpler, as it spins using drag, also having a much lower self starting wind speed requirement. One advantage of a Savonius-type turbine over a Darrieus-type turbine is its low self-start wind speed requirement. Related to this drawback, a new methodology for development of new blade profiles for self-start Darrieus wind turbines is presented in [6], for symmetrical airfoils, and in [7] for asymmetrical airfoils. In [8], a hybrid Darrieus-Savonius turbine is proposed, although the expected performance was not achieved.
A new approach to a type of model called Double Multiple Streamtube (DMS) is presented in [9], aiming to present an easy tool for the study of complex-shaped VAWT blades, able to be integrated in computational modeling tools. Other methods for improving the characteristics of existing turbines without modifications on the turbine itself were studied. A study performed in [10] concluded that an increase in the rotational speed of the VAWT occurs with the installation of an omni-directional-guide-vane, in an urban context, improving the power output and diminishing the visual impact of wind turbines. An alternative to installing omni-directional-guide-vanes is a ducted turbine. In [11], a study was performed, predicting very high specific power output and having clear advantages, namely the visual impact and safety. Despite these advantages, it has a clear drawback, which is its directional sensitivity.
Regarding the implementation of control strategies, a control technique for a sensorless application to maximize the power output up to the rated wind speed of a VAWT using a permanent magnet synchronous generator (PMSG) is proposed in [12], whose results are presented on an equivalent circuit model in the Simulink environment. A model is proposed in [13], considering a variable speed wind turbine, two mass rotor, a PMSG and different power converters, providing a new control strategy based on fractional-order controllers. The model is considered more realistic and accurate, and simulation results show that the proposed control strategy improves disturbance attenuation and system robustness.
In [14], different power converters are studied, modeled and compared in simulation, considering more accurate dynamics regarding the turbine and its connection to the electric network. In [15], variable-speed wind turbines equipped with PMSG are modeled and studied, simulating a malfunction and the controller's response resorting to traditional and fractional PI, comparing their performance.
New blades for a VAWT are developed in [16] with the objective of improving the self-start characteristics without the aid of any additional component. These blades exhibit good performance for wind speeds below 25 m/s, without significant drawbacks such as instability and vibrations. From these blades, a VAWT prototype was developed, which was studied further in [17] with the identification of mathematical models and control development resorting to PID and LQR methods. Further controllers for the same prototype are given in [18], resorting to fractional control, more precisely using the CRONE control method. From the latter work came the need to validate the previous results through testing in a controlled environment such as a wind tunnel.
These works, whilst numerous, prove that technology regarding VAWT applications is still in development and far from being economically viable for their installation in urban areas, which are increasing the energy demands as time passes.

Objectives
Thus, to successfully achieve the objective of control validation, the following stages were accomplished in order:

•
Validate the models previously obtained in [17]  In summary, this paper follows the same line of thought as in [19] with a comparison between simulation and experimental results in order to extract conclusions regarding the system's performance.
This paper presents the Darrieus wind turbine addressed in Section 2. Section 3 presents and discusses the results of the above stages of the work. Section 4 provides possible actions to further improve the results obtained by this work.

Materials and Methods
The following overview of the generator sums up the extensive exposition in [17]. Most VAWTs consist of: with the latter being the main focus of this work. The control system is comprised of sensors, actuators, a central processing unit and controllers, and the flowchart, presented in Figure 1, represents a simplified version of the turbine's wind energy conversion system. While the prototype, together with its components, was developed in previous works as already mentioned, it nevertheless required some maintenance and some components needed replacement. The component acting as the central processing unit, an Arduino UNO, was replaced by an Arduino MEGA 2560 (Arduino, Monza, Italy), seen in Figure 2a, due to its higher processing capabilities compared to the previously used Arduino; the actuators were also replaced by better servos, with higher torque and precision, model Servo TowerPro MG995 (TowerPro, Singapore) 270 • , seen in Figure 2b at the base of the aluminum case. The sensors that will read the input data for the controller are an anemometer and an infrared sensor, responsible for transmitting the wind velocity and turbine's angular velocity, respectively, as seen in Figure 3. These sensors will send their respective digital signals to the Arduino, which in turn will transform the digital signal into the respective wind speed and turbine angular speed values. The anemometer's digital signal conversion occurs with every switch closure, with the rate of 1 switch closure/second = 2.4 km/h [20], while the conversion rate of the infrared sensor's digital signal takes the form of the following equation: since the turbine possesses 5 blades, thus the same number of arms, and ∆t represents the time between consecutive interruptions. The last components that remain to be mentioned are the rheostats. These are simple rotational rheostats, with a rotation angle of [0, 270] • and resistance values of [0, 200] Ω. With these components, the acquisition of new data sets was made possible, with the aim of using the data to validate or create new models resorting to system identification methods. In order to save the data from the sensors, as are read by the Arduino, a laptop was connected to the Arduino through an USB port. Running in the laptop, through Simulink, was a file, as seen in Figure 4, which took the data sets and recorded it in Matlab R2020a's workplace. The two larger subsystems are responsible for the real time conversion of the digital signal read by the Arduino to the values of wind speed intensity and turbine angular speed, corresponding to the top and bottom subsystems, respectively.
In order to minimize possible errors, a procedure was adopted for the acquisition of data sets. This procedure can be seen in the Appendix A. Tables 1 and 2 present lists of the values of resistance and wind speed imposed when collecting the experimental data, together with the corresponding measured values. Errors are small, with small average value (−0.26 Ω and 1.5 × 10 −3 m/s) and small standard deviation (0.92 Ω and 0.013 m/s); nominal values can consequently be used.
An example of a single data set, corresponding to a resistance of 200 Ω, is presented in Figure 5. It can be seen that the sensor's readings are significantly affected by noise; thus, a filter must be applied before the data can be used. The choice fell on a Butterworth filter with cut-off frequency and order determined from the performance index known as root mean square deviation (RMSD): where y t is the unfiltered data, z t is the Butterworth filtered data, with a certain cut-off frequency and order, and T is the number of times the value is observed. As the index is to be applied for different wind velocity values, a normalized root mean square deviation (NRMSD) is calculated in order to be able to compare the data uniformly. The NRMSD is given by: where y t is the average value of the unfiltered data.     Table 2.
By calculating the value of NRMSD as a function of the cut-off frequency, from 0.20 rad/s to 1.00 rad/s with increments of 0.05 rad/s, for several different filter's orders, the following results are obtained, presented by Figure 6, for a wind velocity value of 9 m/s. It is worth mentioning that these results are from the data set presented in Figure 5 and that the results for other data sets are very similar to the one presented in Figure 6, although the results get worse when applied to the extreme values of the wind velocity presented in Table 2, i.e., 6.6 m/s and 12 m/s. Through these results, the desired specifications for the Butterworth filter are: • Cut-off frequency = 0.40 rad/s; • Order = 2 and 3, since both are used for model identification purposes.

Results
This section presents the results from the comparison of previously obtained models with recently identified models resorting to the newly acquired data sets. The comparison of the controllers, the implementation and testing results are also presented.

New Models
Following the same structure as the previously presented models, the new models will be a sum of two terms, one with respect to the wind velocity variation ∆v(s), which is the disturbance variable, and the other to the imposed resistance variation, ∆R c (s), which will be the control variable, whilst ∆ω(s) is the output variable. Matlab's system identification toolbox and command systemIdentification were used for this purpose.

3.2.1.
∆ω(s) ∆v(s) Identified Model Starting with the first term, Table 3 presents the best results from the toolbox. The comparison value will be the fit range of each model, which determines the agreement between the model's response and the measured output, with 0 indicating a poor fit where the model output is the same as the mean value of the measured output [21]. From Table 3, and comparing the Bode diagram and step response of each model, it is concluded that models with best results were the following: • For a 2nd order model: resistance 50 Ω, 2 poles, 1 zero, and filter order 2; r50_2_1_2(s) = 0.682s + 0.01027 s 2 + 0.5484s + 0.007979 , • For a 3rd order model: resistance 40 Ω, 3 poles, 1 zero, and filter order 2; r40_3_1_2(s) = 5.707s + 0.003067 s 3 + 15.87s 2 + 3.675s + 0.002303 . 3.2.2.

∆ω(s)
∆R c (s) Identified Model Whilst the identification of the previous term was performed resorting to the acquired data sets directly, the data sets could not be used directly for this term as a relationship between the resistance's variation and turbine's angular speed had to be created. Data manipulation was required before the use of Matlab R2020a's system identification toolbox could be used.
For the purpose of creating new data sets resorting to the available ones, the following steps were taken: 1.
From the available data, segments where the wind velocity is stable were identified; 2.
With the segments identified, their respective turbine angular speed is removed; 3.
Segments that possess the same wind speed value are joined, creating a continuous line with respect to the desired resistance variation; 4.
The turbine's angular speed is normalized in order to remove the wind speed contribution; 5.
In order to remove the initial signal rise, as it is an undesired dynamic, the minimum resistance value was set as 20 Ω.
The experimental procedure did not include changing the resistance while keeping the wind speed constant because the servos were not available; in fact, identification results were necessary to choose part of the instrumentation of the turbine.
With these steps, new data sets are created where the angular speed varies in function of the resistance's variation. An example of such data sets can be seen in Figure 7, corresponding to the constant wind velocity of 9 m/s. One last step was required for the purpose of simulating a real variation speed for the rheostat. In order to remove the instant changes of values, a slope was created between each resistance value change. The slope's value is given by: where 200 270 refers to the physical limitations of the rheostat itself and 1 60 is the assumed servo speed. This rotational speed for the servo was taken from the technical data sheet of the component [22], according to the value of operating speed with 4.8 V and assuming a load imposed by the rheostat. The final data sets are presented by Figure 8.  From Table 4, and through the comparison of the models' Bode diagram and step response, the models with the best results were: • For the 2nd order model: wind velocity 9 m/s, filter order 2, 2 poles, 1 zero; v9_n2_2_1(s) = 2.399 × 10 −5 s + 1.059 × 10 −6 s 2 + 0.0216s + 0.002461 , • For the 3rd order model: wind velocity 10 m/s, filter order 3, 3 poles, and 2 zeros; v10_n3_3_2(s) = 4.832 × 10 −5 s 2 + 5.513 × 10 −6 s + 2 × 10 −7 s 3 + 0.105s 2 + 0.006388s + 0.0002622 .

Models Comparison
A comparison between the previously obtained models, taken from [17], given by (4) and (5), and the recent models from the previous sections given by: ∆ω(s) = 0.682s + 0.01027 s 2 + 0.5484s + 0.007979 ∆v(s) + 2.399 × 10 −5 s + 1.059 × 10 −6 s 2 + 0.0216s + 0.002461 ∆R c (s), (11) ∆ω(s) = 5.707s + 0.003067 s 3 + 15.87s 2 + 3.675s + 0.002303 ∆v(s) + 4.832 × 10 −5 s 2 + 5.513 × 10 −6 s + 2 × 10 −7 s 3 + 0.105s 2 + 0.006388s + 0.0002622 ∆R c (s), (12) was then performed, imposing the same input signal, and comparing the outputs with values experimentally known. Starting with 2nd order models, the response of each model's terms can be seen in Figure 9. As can be seen, the new model possesses a response closer to the experimental output, represented by the black color, in both plots. This is specially visible in the model's second term, ∆ω ∆R c (s), where the previous model has a response very different from the expected. Regarding the 3rd order models, the results are similar, with the new model's output closer to the expected values, as seen in Figure 10. These results are confirmed in Table 5 where the performance index RMSD, given by (2), is calculated through a simulation resorting to Simulink. It must be mentioned that the simulation results could not perfectly replicate the signals given by the system identification toolbox; thus, the fit range is also mentioned, as it is a performance index calculated by the toolbox, and these results also confirm the results given by RMSD.
Because the new models have a far better performance, new controllers were developed resorting to PID, LQR and CRONE methods, to replace those developed for the models in [17], that had been tested in simulation. As the performance of the 3rd order model exceeds that of the 2nd order model, without significant additional complexity in calculations, only the former model was used for control development.

Control Development
The proposed control objectives are the same as in [17], namely: Maximum settling time of 10 s (t s [2%] ≤ 10 s); • Constant value for the tip speed ratio (TSR), λ = 0.2353, which is the value imposed as optimal in [18].

PID Controller
The transfer function (TF) of the controller is given by: where parameters were provided by Matlab R2020a's PID Tuner application [23], according to the desired response's specifications imposed by the user, set to a response time of 1.97 s and a damping coefficient value of 0.9.
Introducing these values into the PID block on the Simulink file and imposing a step signal as an input, the obtained results are seen in Figure 11

LQR Controller
The LQR controller minimizes a performance index J given by where x is the state vector of the state-space representation of the plant, u is the input vector, and Q is a matrix determining the relative weight of x and u in the minimization. Several values for Q were used, but the results turned out to be systematically worse than those of PID's, with controllers not providing acceptable values of settling time and overshoot, and struggling instead to achieve the desired TSR value. Table 6 presents several examples and the system's characteristics when applied in simulation.

CRONE Controller
There are three generations of CRONE controllers (from the French acronym Commande Robuste d'Ordre Non-Entier), having in common that the justification of their design is based upon fractional derivatives [24,25]. The resulting controller, however, has a TF that only resorts to integer order derivatives, i.e., a TF, which is a rational function in s. Third generation CRONE controllers are used when the model of a plant is known with uncertainty, as is the case of the wind turbine [26]. The purpose of the controller is to ensure a desired behavior for the nominal plant, and minimize performance changes if model parameters vary within the expected ranges of uncertainty. For the development of this controller, a Matlab toolbox created by the CRONE Research Group of Université de Bordeaux was used.
To develop this type of controller, extreme plants must be defined, corresponding to the maximum effects of uncertainty in the model. For that purpose, two types of uncertainties were considered, parametric and pole/zero location. Through the calculation of NRMSD, given by (3), it was found that the parametric uncertainty possessed better results, as the corresponding dispersion was more localized, and the resulting controller more robust. The extreme plants are given by: l (s) = 4.107 × 10 −5 s 2 + 4.686 × 10 −6 s + 1.7 × 10 −7 s 3 + 0.08925s 2 + 0.00543s + 0.0002229 , G 3 u (s) = 5.557 × 10 −5 s 2 + 6.34 × 10 −6 s + 2.3 × 10 −7 s 3 + 0.1207s 2 + 0.007346s + 0.0003015 , (17) corresponding to the lower and upper bound extreme plants, respectively.
With the nominal and extreme plants determined, the toolbox can be used to develop CRONE controllers, with the choice of desired specifications introduced by the user [27]. The developed controllers are presented in Table 7, where the value of "Desired t s " is chosen by the user and defines the settling time t s [5%]. The results show a significant discrepancy between the desired and the actual settling time, which may be due to the incapacity of the plant to respond as fast as desired, when the specification t s is low. Taking into account the results obtained, the chosen CRONE controller is the one with t s set for 700 s, given by: The corresponding response can be seen in Figure 12.

Controller's Performance Comparison
From the responses for an input signal in the form of a step, with the characteristics given by Table 8, it is clear that the LQR controller possesses the worst performance; as such it will be omitted from other input signals [28]. Changing the input signal from a step signal to a ramp with slope 1/120, in order to simulate a unitary increase in wind velocity in 2 minutes, the results are as shown in Figure 13. In order to simulate a real input signal, the last input signal is corrupted with additive white noise, and the results become those in Figure 14. As can be observed, the PID controller has a poorer response for ramp signals, whether corrupted with noise or not, whilst the CRONE controller maintains the rheostat's behavior. Thus, this was the controller chosen for implementation.

System Discretization
The chosen method of discretization was Tustin's, in order to preserve stability in conversion (which is not ensured by simpler methods, like rectangular approximations), and decrease the resulting bias due to the discretization. For the purpose of choosing a sample time value, its effect on the controller's behavior was studied, with the results presented in Figure 15. A sampling time of t s = 10 s is seen to include the bandwidth of the controller. Time response simulations show that the overshoot and settling time of the rheostat would be more significant with lower sampling times, so there is no advantage in choosing a shorter sample time. As such, that value of the sampling time was chosen, resulting in discrete TFs given by: which were then applied by the Arduino.

Controller Testing
In order to test the controller, a control problem must be defined first. For that purpose, the system's overall TSR is studied, which is given by the following equation: where r = 0.173 m, defining the radius of the turbine, v is the wind velocity intensity, and ω defines the turbine's angular velocity from each term, v or resistance R c . Through simulation results, the system's minimum and maximum TSR values were found, as functions of the wind velocity, as presented by Figure 16, which also shows the expected actuation range of the wind tunnel. As the maximum TSR distribution is significantly different from the minimum distribution, it was decided to set the reference value as the system's TSR mean value in order to impose a reaction from the controller with either an increase or decrease in wind velocity intensity.
The results of the tests in the wind tunnel are presented in Figure 17, and compared with simulation results, for the same wind velocity input signal, in Figure 18.
From Figure 17, there are two instances of interest: the first occurring around t ≈ 900 s and the second at about t ≈ 1750 s. These instances correspond to a decrease and an increase in the perceived wind velocity intensity by the anemometer, respectively, which prove that the controller works as intended.
For instance, for a decrease in the perceived wind velocity intensity, the mean TSR is increased, but due to the delay caused by the discrete TF's sampling time, the reference TSR is still the prior to the change of wind velocity, which is now lower than the desired TSR. Thus, the controller lowers the rheostat's value, corresponding to the sudden fall visible at t ≈ 900 s in Figure 17b. The same occurs for an increase of wind velocity intensity, with a sudden increase in the rheostat's value seen at t ≈ 1750 s. Comparing the results in the simulation from the results obtained in the wind tunnel, it can be observed, from Figure 18, that the results of simulation are worse than the results from the actual controller's implementation, with higher overshoot and settling time. This may be due to a faster implementation of the algorithm in the Arduino, when compared to the numerical approximation used by Simulink.

Conclusions
With this work, the technology readiness level is increased to 4. Although the controller's implementation is considered a success, even surpassing the expected performance shown by simulation results, it might be further improved. For instance, with the presence of a reliable controller and some programming, it is now possible to obtain dedicated data sets for the second term of the prototype, i.e., the term of ∆ω ∆R c (s), instead of using the manipulated data.
Other minor improvements are also possible to occur with simple maintenance of the turbine, such as the repair of superficial flaws of the blades, correct alignment of the blades and the correct alignment of the generator's rotational axis.
Author Contributions: G.R.: methodology, software, investigation, data curation, writing-original draft preparation, visualization; D.V.: conceptualization, methodology, validation, investigation, data curation, writing-review and editing, supervision; R.M.: conceptualization, methodology, validation, investigation, data curation, writing-review and editing, supervision. All authors have read and agreed to the published version of the manuscript.

Acknowledgments:
The authors express their gratitude to Agostinho Fonseca for the availability and aid provided for this work. This work was supported by FCT, through IDMEC, under LAETA project UIDB/50022/2020; and by FCT, through ICT, project UIDB/04683/2020.

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