timeseries.m

小波神经网络的matlab代码

%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 <nCreate> to create new NNs, train them then perform forecast 
      if(scheme == '1') 
           nCreate; 
      end 

      %If scheme is 'r', call <nRetrain> to retrain the existing NNs 
      if(scheme == '2') 
           nRetrain; 
      end 

      %If scheme is 'f', call <nForecast> 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 <point 2>, 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);