www.pudn.com > timeseries.rar > timeseries.m
%File name : nprogram.m
%Description : This file reads the data from its source into their respective matrices prior to
% performing wavelet decomposition.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Clear command screen and variables
clc;
clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% user desired resolution level (Tested: resolution = 2 is best)
level = menu('Enter desired resolution level: ', '1',...
'2 (Select this for testing)', '3', '4');
switch level
case 1, resolution = 1;
case 2, resolution = 2;
case 3, resolution = 3;
case 4, resolution = 4;
end
msg = ['Resolution level to be used is ', num2str(resolution)];
disp(msg);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% user desired amount of data to use
data = menu('Choose amount of data to use: ', '1 day', '2 days', '3 days', '4 days',...
'5 days', '6 days', '1 week (Select this for testing)');
switch data
case 1, dataPoints = 48; %1 day = 48 points
case 2, dataPoints = 96; %2 days = 96 points
case 3, dataPoints = 144; %3 days = 144 points
case 4, dataPoints = 192; %4 days = 192 points
case 5, dataPoints = 240; %5 days = 240 points
case 6, dataPoints = 288; %6 days = 288 points
case 7, dataPoints = 336; %1 weeks = 336 points
end
msg = ['No. of data points to be used is ', num2str(dataPoints)];
disp(msg);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Menu for data set selection
select = menu('Use QLD data of: ', 'Jan02',...
'Feb02', 'Mar02 (Select this for testing)', 'Apr02', 'May02');
switch select
case 1, demandFile = 'DATA200601_QLD1';
case 2, demandFile = 'DATA200602_QLD1';
case 3, demandFile = 'DATA200603_QLD1';
case 4, demandFile = 'DATA200604_QLD1';
case 5, demandFile = 'DATA200605_QLD1';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Reading the historical DEMAND data into tDemandArray
selectedDemandFile=[demandFile,'.csv'];
[regionArray, sDateArray, tDemandArray, rrpArray, pTypeArray] ...
= textread(selectedDemandFile, '%s %q %f %f %s', 'headerlines', 1, 'delimiter', ',');
%Display no. of points in the selected time series demand data
[demandDataPoints, y] = size(tDemandArray);
msg = ['The no. of points in the selected Demand data is ', num2str(demandDataPoints)];
disp(msg);
%Decompose historical demand data signal
[dD, l] = swtmat(tDemandArray, resolution, 'db2');
approx = dD(resolution, :);
%Plot the original demand data signal
figure (1);
subplot(resolution + 2, 1, 1); plot(tDemandArray(1: dataPoints))
legend('Demand Original');
title('QLD Demand Data Signal');
%Plot the approximation demand data signal
for i = 1 : resolution
subplot(resolution + 2, 1, i + 1); plot(approx(1: dataPoints))
legend('Demand Approximation');
end
%After displaying approximation signal, display detail x
for i = 1: resolution
if( i > 1 )
detail(i, :) = dD(i-1, :)- dD(i, :);
else
detail(i, :) = tDemandArray' - dD(1, :);
end
if i == 1
subplot(resolution + 2, 1, resolution - i + 3); plot(detail(i, 1: dataPoints))
legendName = ['Demand Detail ', num2str(i)];
legend(legendName);
else
subplot(resolution + 2, 1, resolution - i + 3); plot(detail(i, 1: dataPoints))
legendName = ['Demand Detail ', num2str(i)];
legend(legendName);
end
i = i + 1;
end
%Normalising approximation demand data
maxDemand = max(approx'); %Find largest component
normDemand = approx ./ maxDemand; %Right divison
maxDemandDetail = [ ];
normDemandDetail = [, ];
detail = detail + 4000;
for i = 1: resolution
maxDemandDetail(i) = max(detail(i, :));
normDemandDetail(i, :) = detail(i, :) ./maxDemandDetail(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Reading the historical historical PRICE data into rrpArray
selectedPriceFile = [demandFile, '.csv'];
[regionArray, sDateArray, tDemandArray, rrpArray, pTypeArray] ...
= textread(selectedDemandFile, '%s %q %f %f %s', 'headerlines', 1, 'delimiter', ',');
%Display no. of points in Price data
[noOfDataPoints, y] = size(rrpArray);
msg = ['The no. of points in Price data data is ', num2str(noOfDataPoints)];
disp(msg);
%Decompose historical Price data signal
[dP, l] = swtmat(rrpArray, resolution, 'db2');
approxP = dP(resolution, :);
%Plot the original Price data signal
figure (2);
subplot(resolution + 3, 1, 1); plot(rrpArray(2: dataPoints))
legend('Price Original');
title('Price Data Signal');
%Plot the approximation Price data signal
for i = 1 : resolution
subplot(resolution + 3, 1, i + 1); plot(approxP(2: dataPoints))
legend('Price Approximation');
end
%After displaying approximation signal, display detail x
for i = 1: resolution
if( i > 1 )
detailP(i, :) = dP(i-1, :)- dP(i, :);
else
detailP(i, :) = rrpArray' - dP(1, :);
end
if i == 1
[B,A]=butter(1,0.65,'low');
result =filter(B,A, detailP(i, 1: dataPoints));
subplot(resolution + 3, 1, resolution - i + 4);plot(result(i, 2: dataPoints))
legendName = ['low pass filter', num2str(i)];
legend(legendName);
subplot(resolution + 3, 1, resolution - i + 3); plot(detailP(i, 2: dataPoints))
legendName = ['Price Detail ', num2str(i)];
legend(legendName);
else
subplot(resolution + 3, 1, resolution - i + 3); plot(detailP(i, 2: dataPoints))
legendName = ['Price Detail ', num2str(i)];
legend(legendName);
end
i = i + 1;
end
%Normalising approximation Price data
maxPrice = max(approxP'); %Find largest component
normPrice = approxP ./ maxPrice; %Right divison
maxPriceDetail = [ ];
normPriceDetail = [, ];
detailP = detailP + 40;
for i = 1: resolution
maxPriceDetail(i) = max(detailP(i, :));
normPriceDetail(i, :) = detailP(i, :) ./maxPriceDetail(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Function here allows repetitive options to,
% 1) Create new NNs, 2) Retrain the existing NNs,
% 3) Perform load demand forecasting and 4) Quit session
while (1)
choice = menu('Please select one of the following options: ',...
'CREATE new Neural Networks',...
'Perform FORECASTING of load demand', 'QUIT session...');
switch choice
case 1, scheme = '1';
case 2, scheme = '2';
case 3, scheme = '3';
case 4, scheme = '4';
end
%If scheme is 'c', call to create new NNs, train them then perform forecast
if(scheme == '1')
nCreate;
end
%If scheme is 'r', call to retrain the existing NNs
if(scheme == '2')
nRetrain;
end
%If scheme is 'f', call to load the existing NN model
if(scheme == '3')
nForecast;
end
%If scheme is 'e', verifies and quit session if 'yes' is selected else continue
if(scheme == '4')
button = questdlg('Quit session?', 'Exit Dialog','Yes','No','No');
switch button
case 'Yes', disp(' ');
disp('Session has ended!!');
disp(' ');
break;
case 'No', quit cancel;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%File name : ncreate.m
%Description : This file prepares the input & output data for the NNs. It creates then trains the
% NNs to mature them.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Clear command screen and set target demand ouput to start at point 2
clc;
targetStartAt = 2;
disp('Program will now CREATE a Neural Network for training and forecasting...');
disp(' ');
disp('To capture the pattern of the signal, the model is ')
disp('set to accept dataPoints x 2 sets of training examples; ');
disp('1 set of demand + 1 sets of price. ');
disp(' ');
disp('The normalised demand data , is to be taken as the ')
disp('output value for the first iteration of training examples. ');
disp(' ');
disp('Press ENTER key to continue...');
pause;
%Clear command screen then prompt for no. of training cycles
%For current program, 1 cycle = user set no. of iterations (ie: dataPoints)
clc;
cycle = input('Input the number of training cycles: ');
numOfTimes = resolution + 1;
%Creating and training the NNs for the respective
%demand and price coefficient signals
for x = 1: numOfTimes
%Clearing variables
clear targetDemand;
clear inputs;
clear output;
clc;
if(x == 1)
neuralNetwork = ['Neural network settings for approximation level ',
num2str(resolution)];
else
neuralNetwork = ['Neural network settings for detail level ', num2str(x - 1)];
end
disp(neuralNetwork);
disp(' ');
%Set no. of input nodes and hidden neurons for the
%respective demand and price coefficient signal
numOfInputs = 2;
inputValue = ['Number of neural network INPUT units is set at ', num2str(numOfInputs)];
disp(inputValue);
disp(' ');
numOfOutput = 1;
outValue = ['Output is set to ', num2str(numOfOutput)];
disp(outValue);
disp(' ');
numOfHiddens = input('Enter the no. of HIDDEN units for the NN hidden : ');
hiddenValue = ['Number of neural network HIDDEN units is set at ',
num2str(numOfHiddens)];
disp(hiddenValue);
disp(' ');
%Setting no. of training examples
trainingLength = dataPoints;
%Set target outputs of the training examples
if(x == 1)
targetDemand = normDemand(targetStartAt: 1 + trainingLength);
else
targetDemand = normDemandDetail(x - 1, targetStartAt: 1 + trainingLength);
end
%Preparing training examples
%Setting training i/ps to be 2 with user defined no. of iterations (dataPoints)
y = 0;
while y < trainingLength
if(x == 1)
inputs(1, y + 1) = normDemand(y + 1);
inputs(2, y + 1) = normPrice(y + 1);
else
inputs(1, y + 1) = normDemandDetail(x - 1, y + 1);
inputs(2, y + 1) = normPriceDetail(x - 1, y + 1);
end
output(y + 1, :) = targetDemand(y + 1);
y = y + 1;
end
inputs = (inputs');
%Setting no. of training cycles
[ni, np] = size(targetDemand); % <== [ni, np] tells the NN how long is 1 cycle;
size(targetDemand)
clear nn;
%Create new neural network for respective coefficient signal
%NET = MLP(NIN, NHIDDEN, NOUT, FUNC)
nn = mlp(numOfInputs, numOfHiddens, numOfOutput, 'linear');
%NN options
options = zeros(1, 18);
options(1) = 1; %Provides display of error values
options(14) = cycle * ni * np;
%Training the neural network
%netopt(net, options, x, t, alg);
nn = netopt(nn, options, inputs, output, 'scg');
%Save the neural network
filename = ['nn', num2str(x)];
save(filename, 'nn');
disp(' ');
msg = ['Neural network successfully CREATED and saved as => ', filename];
disp(msg);
if(x < 3)
disp(' ');
disp('Press ENTER key to continue training the next NN...');
else
disp(' ');
disp('Model is now ready to forecast!!');
disp(' ');
disp('Press ENTER key to continue...');
end
pause;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%File name : nforecast.m
%Description : This file loads the trained NNs for load forcasting and %recombines the predicted
% outputs from the NNs to form the final predicted load series.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Clear command screen and variables
clc;
clear forecastResult;
clear actualDemand;
clear predicted;
clear actualWithPredicted;
disp('Neural networks loaded, performing electricity demand forecast...');
disp(' ');
pointsAhead = input('Enter no. of points to predict (should be < 336): ');
numOfTimes = resolution + 1;
%Predict coefficient signals using respective NNs
for x = 1 : numOfTimes
%Loading NN
filename = ['nn', num2str(x)];
clear nn
load(filename);
clear in;
numOfInputs = nn.nin;
y = 0;
%Preparing details to forecast
while y < pointsAhead
if(x == 1)
propData(1, y + 1) = normDemand(y + 1);
propData(2, y + 1) = normPrice(y + 1);
else
propData(1, y + 1) = normDemandDetail(x - 1, y + 1);
propData(2, y + 1) = normPriceDetail(x - 1, y + 1);
end
y = y + 1;
end
if(x == 1)
forecast = mlpfwd(nn, propData') * maxDemand;
else
forecast = mlpfwd(nn, propData') * maxDemandDetail(x - 1) - 4000;
end
forecastResult(x, :) = forecast';
end
%For debugging
% forecastResult
actualDemand = tDemandArray(2: 1 + pointsAhead);
predicted = sum(forecastResult, 1)';
%Calculating Mean Absolute Error
AbsError = abs(actualDemand - predicted(1: pointsAhead)) ./ actualDemand;
msg = ['Mean Absolute Error = ', num2str(mean(AbsError(1: pointsAhead))), ' !!'];
disp(' ');
disp(msg);
%Plot actual time series against predicted result
figure(3)
actualWithPredicted(:, 1) = actualDemand;
actualWithPredicted(:, 2) = predicted(1: pointsAhead);
plot(actualWithPredicted);
graph = ['Mean Absolute Error = ', num2str(mean(AbsError))];
title(graph);
legend('Actual', 'Forecasted');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%File name : nretrain.m
%Description : This file loads the existing NNs and trains them again.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
%Prompt for the starting point for training
disp('Program will now RETRAIN the Neural Networks ')
disp('with the SAME intial data series again...');
disp(' ');
disp('To capture the pattern of the signal, the model is ')
disp('set to accept dataPoints x 2 sets of training examples; ');
disp('1 set of demand + 1 sets of price. ');
disp(' ');
disp('The normalised demand data , is to be taken as the ')
disp('output value for the first iteration of training examples. ');
disp(' ');
msg = ['Data points to be used for reTraining the NNs is from 1 to ',...
num2str(dataPoints)];
disp(msg);
disp(' ');
disp('Press ENTER key to continue...');
pause;
%Clear command screen
clc;
%Prompt for no. of training cycles
%For current program, 1 cycle = user set no. of iterations (ie: dataPoints)
cycle = input('Input number of cycles to retrain the NNs: ');
numOfTimes = resolution + 1;
%Loading existing NNs for training
for x = 1: numOfTimes
%Re-initialising variables
clear targetDemand;
clear inputs;
clear output;
clc;
%Loading NN for the respective demand and temperature coefficient signals
filename = ['nn', num2str(x)];
clear nn
load(filename);
%Getting the size of NN
numOfInputs = nn.nin;
numOfHiddens = nn.nhidden;
numOfOutput = 1;
%Setting length of reTraining examples and target outputs
reTrainLength = dataPoints;
targetLength = reTrainLength;
targetStartAt = 2;
%Set target outputs of the training examples
if(x == 1)
targetDemand = normDemand(targetStartAt: 1 + targetLength);
else
targetDemand = normDemandDetail(x - 1, targetStartAt: 1 + targetLength);
end
%Preparing training examples
%Setting training i/ps to be 2 with user set no. of iterations (dataPoints)
y = 0;
while y < reTrainLength
if(x == 1)
inputs(1, y + 1) = normDemand(y + 1);
inputs(2, y + 1) = normPrice(y + 1);
else
inputs(1, y + 1) = normDemandDetail(x - 1, y + 1);
inputs(2, y + 1) = normPriceDetail(x - 1, y + 1);
end
output(y + 1, :) = targetDemand(y + 1);
y = y + 1;
end
inputs = (inputs');
%Setting no. of training cycles
[ni, np] = size(targetDemand); % <== [ni, np] tells the NN how long is 1 cycle;
size(targetDemand) %With reference to line 106
%NN options
options = zeros(1, 18);
options(1) = 1; %Provides display of error values
options(14) = cycle * ni * np;
%Training the neural network
%netopt(net, options, x, t, alg);
nn = netopt(nn, options, inputs, output, 'scg');
%Save the neural network
filename = ['nn', num2str(x)];
save(filename, 'nn');
disp(' ');
msg = ['Neural network => ', filename, ' <= successfully RETRAINED and saved!! '];
disp(msg);
if(x < 3)
disp(' ');
disp('Press ENTER key to continue training the next NN...');
else
disp(' ');
disp('Model is now ready to forecast again!!');
disp(' ');
disp('Press ENTER key to continue...');
end
pause;
end