% QRLSTM
% 数据集 列为特征，行为样本数目
clc
clear
load Train.mat
% 训练集
Train(1,:) =[];
y = Train.demand(1:2000,:);
x = Train{1:2000,3:end};
clearvars -except x y
% 数据归一化
[xnorm,xopt] = mapminmax(x',0,1);
xnorm = mat2cell(xnorm,size(xnorm,1),ones(1,size(xnorm,2)));
[ynorm,yopt] = mapminmax(y',0,1);
ynorm = ynorm';
%%
% LSTM 层设置，参数设置
inputSize = size(xnorm{1},1); %数据输入x的特征维度
outputSize = size(ynorm,2); %数据输出y的维度
numhidden_units1=50; %网络单元
numhidden_units2= 20;
numhidden_units3=100;
%%
% lstm
for i = 0.1:0.1:0.9
layers = [ ...
sequenceInputLayer(inputSize,'name','input') %输入层设置
lstmLayer(numhidden_units1,'Outputmode','sequence','name','hidden1')
dropoutLayer(0.3,'name','dropout_1')
lstmLayer(numhidden_units2,'Outputmode','last','name','hidden2')
dropoutLayer(0.3,'name','drdiopout_2')
fullyConnectedLayer(outputSize,'name','fullconnect') % 全连接层设置（影响输出维度）（cell层出来的输出层） %
quanRegressionLayer('out',i)];
% 参数设定
opts = trainingOptions('adam', ...
'MaxEpochs',10, ...
'GradientThreshold',1,...
'ExecutionEnvironment','cpu',...
'InitialLearnRate',0.001, ...
'LearnRateSchedule','piecewise', ...
'LearnRateDropPeriod',2, ... %2个epoch后学习率更新
'LearnRateDropFactor',0.5, ...
'Shuffle','once',... % 时间序列长度
'SequenceLength',1,...
'MiniBatchSize',24,...
'Verbose',0);
%
% 网络训练
tic
net1 = trainNetwork(xnorm,ynorm,layers,opts);
net2(floor(i*10)) = trainNetwork(xnorm,ynorm,layers,opts);
end
% 测试集
load Test.mat
Test(1,:) =[];
clear x y
y = Test.demand;
x = Test{:,3:end};
% 归一化
[xnorm,xopt] = mapminmax(x',0,1);
xnorm = mat2cell(xnorm,size(xnorm,1),ones(1,size(xnorm,2)));
[ynorm,yopt] = mapminmax(y',0,1);
ynorm = ynorm';
color=[111,168,86;128,199,252;112,138,248;184,84,246]/255;
% 迭代计算
% 点预测
for i = 1:length(xnorm)
YpredictNorm1= net1.predict(xnorm(i));
Ypredict1(i) = mapminmax('reverse',YpredictNorm1,yopt);
end
% 区间预测
for i = 1:length(net2)
YpredictNorm2(i,:) = net2(i).predict(xnorm);
Ypredict2(i,:) = mapminmax('reverse',YpredictNorm2(i,:),yopt);
end
%% 可视化
x =[1:168];
figure
fill([x,x(end:-1:1)],[smooth(Ypredict2(1,:),1)',smooth(Ypredict2(end,end:-1:1),1)'],'r','FaceColor',[1 0.8 0.8],'EdgeColor','none')
hold on
plot(Ypredict1,'-','Color',color(3,:),'LineWidth',1)
hold on
plot(y','-','Color',color(4,:),'LineWidth',1)
legend('90%置信区间','预测值','实际值')
title('区间预测')
xlabel('Time')
ylabel('load(KW)')
title('预测周负荷概率分布')
time =datetime('2022-10-1'):datetime('2022-10-7');
xticks([0:24:24*7]);
xticklabels(string(time));
xtickangle(45)
legend('90%置信区间','预测值','实际值')
picp = PICP(Ypredict2,y);
pimwp = PIMWP(Ypredict2,y);
% 区间覆盖率
function picp = PICP(Ypredict,Yreal)
% 维度保持一致（样本个数）
if size(Ypredict,1) ~= size(Yreal,1)
Ypredict = Ypredict';
end
RangeForm = [ Ypredict(:,1),Ypredict(:,end)];
Num = 0;
for i = 1:length(Yreal)
Num = Num + (Yreal(i) >= RangeForm(i,1) && Yreal(i) <= RangeForm(i,2));
end
picp = Num / length(Yreal);
end
% 区间平均宽度百分比
function pimwp = PIMWP(Ypredict,Yreal)
% 维度保持一致（样本个数）
if size(Ypredict,1) ~= size(Yreal,1)
Ypredict = Ypredict';
end
pimwp = 1/length(Yreal) * sum( (Ypredict(:,end) - Ypredict(:,1))...
./ Yreal) * 100;
end
