Code covered by the BSD License  

Highlights from
Example files for "Programming with MATLAB" Webinar

image thumbnail

Example files for "Programming with MATLAB" Webinar

by

 

Example files for "Programming with MATLAB" Webinar, first delivered October 3, 2013.

WindAnalysisFcn.m
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;

Contact us