function windData = WindAnalysisFcn(filename)
%% Wind Turbine Data Analysis
% This demo analyzes wind data measured on a meteorological observation
% tower to see if the location is a good prospect for a wind turbine. Data
% is from three different wind sensors at 80m. Temperature is also recorded
% at 3m height. Data is logged every hour.
% Copyright 2013 The MathWorks, Inc.
%% Handle No-Input Case
if nargin == 0
% If user doesn't supply a file name, request one
[fname, pname] = uigetfile('*.txt');
filename = fullfile(pname, fname);
else
% Check to make sure input is a string
validateattributes(filename, {'char'}, {'row'})
% Check to make sure file exists
if ~exist(filename, 'file')
error('File "%s" cannot be found', filename);
end
end
%% Read in Turbine Data from Text File
% Function autogenerated from import tool
[~, fname, ~] = fileparts(filename);
[time,velS1,velS2,velS3,tempS] = importWindDataFcn(filename);
%% Average Wind Speed for Different Sensors
velAvg = mean([velS1, velS2, velS3], 2);
%% Determine Icing Conditions
% Remove any readings effected by icing which results in extremely
% inaccurate values biased towards zero. Icing conditions are when
% tavg < tIce and vavg < vIce.
tIce = 2;
vIce = 1;
% Comparing to the critical values
idxTIce = tempS < tIce;
idxVIce = velAvg < vIce;
idxIce = idxTIce & idxVIce ;
% Remove values related to icing
time(idxIce) = [];
tempS(idxIce) = [];
velAvg(idxIce) = [];
%% Distribution of Wind Speeds at Hub Height
% Use a weibull distribution to fit the distribution of wind speeds which
% is known to often give a good fit to wind speed distributions.
% Plotting a histogram of wind speeds
dv = 0.5;
vbins = 0:dv:ceil(max(velAvg));
% Plotting a histogram of wind speeds
figure; hist(velAvg, vbins);
xlabel('wind velocity (m/s)'); ylabel('count');
title(['Wind Speed Distribution from ', fname]);
% Calculate probability of wind speed range
nelements = hist(velAvg, vbins);
probvbins = nelements/sum(nelements); % Probability of a given velocity range
%% Defining the Turbine Power Curve
% To calculate average turbine power and capacity factor, we need to make
% some assumptions regarding the wind turbine model and its power curve.
% We will assume a 1MW wind turbine and the following turbine power
% curve.
% Turbine power curve coefficents
prated = 1e6; % wind turbine rated power (W)
vin = 2; % cut-in speed (m/s)
vr = 14; % rated output speed (m/s)
vout = 25; % cut-out speed (m/s)
% Calculating power curve
powervbins = prated*(vbins.^2 - vin^2)/(vr^2 - vin^2);
powervbins(vbins <= vin) = 0;
powervbins(vbins > vout) = 0;
powervbins(vbins >= vr & vbins <= vout) = prated;
%% Calculating Average Turbine Power and Capacity Factor
% Capacity factor is ratio of the actual output of a turbine over a period
% of time and its potential output if it had operated at full capacity the
% entire time. Typical capacity factor range from 20-50% depending on
% location and wind turbine.
% Integrate power at given velocity * velocity probability distribution
% function over range of possible velocities
% Calculate average power by summing the products of vel probability and power
avgPower = sum(probvbins .* powervbins); % (W)
% Calculating capacity factor (average power / rated power)
cf = avgPower / prated;
disp(['Assumed wind turbine rated power (MW): ', num2str(prated/1e6)]);
disp(['Averaged turbine power (kW): ', num2str(avgPower/1e3)]);
disp(['Capacity factor (%): ', num2str(cf*100)]);
%% Combine into a Structure
% We will combine the information into a structure for easy data
% management.
windData.date = time;
windData.temp = tempS;
windData.vavg = velAvg;
windData.avgPower = avgPower;
windData.cf = cf;
windData.filename = fname;