模拟负荷不确定性——拉丁超立方抽样生成及缩减场景（Matlab全代码）

WW、forever

已于 2023-03-02 22:38:29 修改

阅读量3.2k

点赞数 8

分类专栏：论文复现：模拟风电、光伏和负荷不确定性文章标签： matlab 负荷不确定性负荷预测场景缩减

于 2023-03-02 22:34:51 首次发布

本文链接：https://blog.csdn.net/qq_44246618/article/details/129300796

版权

论文复现：模拟风电、光伏和负荷不确定性专栏收录该内容

3 篇文章 60 订阅

订阅专栏

本文介绍了一种通过拉丁超立方抽样在MATLAB中模拟负荷不确定性并进行场景缩减的方法。代码涵盖场景生成、缩减过程，还提供了如何更换分布函数（如从正态分布改为Weibull分布）的指导。实验结果显示了1000种场景及缩减后的场景概率，适用于高精度（约5%预测误差）的负荷预测研究。

摘要由CSDN通过智能技术生成

与风电、光伏不确定性不同，负荷不确定性的分布函数、预测误差有别于风电、光伏出力的不确定性，本代码通过拉丁超立方抽样和快速前代消除法模拟了负荷的不确定性，并提供了如何修改分布函数/概率密度函数的思路。

据文献表明，负荷预测准确率较高，约为预测负荷的5%。本文提供全代码及修改为其他函数分布的思路。

1、MATLAB代码

1.1 场景生成

clc 
clear
%% 参数设置
Z = 24;                                     % 以一天24h为周期进行场景生成
iterations = 1000;                          % 抽样次数 1000次-1000种场景
S = 10;                                     % 代表性场景数量 5个
load('Load_forecast.mat')
step = 24;                                  % 步长 h（总计一天）
P_test    = Load_forcast;                   % 转置-待生成数据
mu        = 0*ones(Z,1);                    % 正态分布-均值
sigma     = 0.05*P_test';                   % 正态分布-标准差（此时取5%Pforecast）
leng      = 50;
%% 场景生成
% 概率密度函数
y1 = @(x) exp(-(x-mu)^2/(2*sigma^2))/(sigma*(2*pi)^0.5);
% 累计分布函数
CDF = @(x)@(mu)@(sigma) normcdf(x,mu,sigma);
x = leng : 0.1 : leng;
y2 = CDF(x);
% 拉丁超立方抽原始样本
segmentSize = 1 / iterations;                                              % 每层大小
for j = 1:Z
for i = 0 : iterations-1                                                   % 逐层随机抽样
    segmentMin = i * segmentSize;
    segmentMax = (i+1) * segmentSize;
    samplePoint(j,i+1) = segmentMin + rand() * segmentSize;
end
end
for i= 1:Z
L_KN(i,:) = round(rand(1,iterations)*iterations);
end
row   = corr(L_KN');                                                       % 相关系数矩阵
[x,y] = eig(row);                                                          % 求特征值，判断是否是正定矩阵（特征值全部大于0） 
% 若是正定矩阵，才可以后续
D = chol(row,'lower');                                                     % L即为非奇异下三角矩阵
G_KN = D'*L_KN;

row_s = corr(samplePoint');                                                % 相关系数矩阵
D_s = chol(row_s,'lower');
paixu_before = D_s*G_KN;
for i = 1:Z
paixu(i,:) = Order_Vector(paixu_before(i,:));                              % 上式的顺序矩阵，每行从大到小排序
samplePoint(i,:) = samplePoint(i,randperm(iterations));                    % 乱序得到原始样本
% samplePoint(i,:) = samplePoint(i,paixu(i,:));                            % 乱序得到原始样本
end
% 映射得到最终样本
for i=1:Z
    for j=1:iterations
        Value(i,j) = ICDF(samplePoint(i,j),mu(i),sigma(i));
        if Value(i,j)>2*sigma(i)                                           % 抽样生成的比较随机，故进行上下界限制
            Value(i,j)=2*sigma(i);
        end
        if Value(i,j)<-2*sigma(i)
            Value(i,j)=-2*sigma(i);
        end
    end
end

for i=1:Z
    for j=1:iterations
        a_thou(i,j)=P_test(i)+Value(i,j);
%         plot(Pw(i)+Value(i,j));
    end
end
figure
hold on
for j=1:iterations
    plot(a_thou(:,j));
end
plot(P_test,'r-s','linewidth',2)
plot(P_test+2*sigma','b-s','linewidth',1.5)
plot(P_test-2*sigma','b-s','linewidth',1.5)

1.2 场景缩减

%% 场景削减
Ws = Value; 
m = iterations;
Ws_d = Ws;
% 场景削减
p_i = 1/m*ones(1,m); 
x=zeros(m,m); 
for i=1:m
    for j=1:m
            x(i,j)=sum(abs(Ws(:,i)-Ws(:,j)));
    end
end

% 计算每个场景与剩余场景的概率距离之和y
y0=zeros(1,m);
for i=1:m
y0(i)=1/m*sum(x(:,i));
end
k=length(y0);

% 不断削减场景，直到剩余S个场景
while(k>S)
d=find(y0==min(y0));                            % 选定与剩余场景的概率距离之和最小的场景
x_2=x+5000*eye(k);                             % (5000敏感)构造新的x，以便找出风电场景Ws中与场景d几何距离最小的场景r
r=find(x_2(:,d)==min(x_2(:,d)));
p_i(r)=p_i(r)+p_i(d);                           % 将d场景的概率加到r场景上
% 在场景中删除d场景
p_i(d)=[]; 
Ws_d(:,d)=[];
x(d,:)=[];
x(:,d)=[];
y0(d)=[];
k=length(y0);
end

figure
[ss,gg]=meshgrid(1:S,1:Z );
plot3(ss,gg,Ws_d,'-');

grid
xlabel('场景');
ylabel('时刻');
zlabel('出力值');
title('场景削减图');
figure
hold on
for j=1:S
    h(j)=plot(P_test'+Ws_d(:,j));                % 五种场景
end
h = legend( h,'1','2','3','4','5','6','7','8','9','10');
for j=1:S
    a_five(j,:)=P_test'+Ws_d(:,j);               % 五种场景
end

figure
% axes2=axes('position',[0.12,0.12,0.6,0.6]); 
bar(p_i,0.5);
a_pro=p_i;
set(gca,'XTickLabel',{'场景1','场景2','场景3','场景4','场景5','场景6','场景7','场景8','场景9','场景10'});
for i=1:S 
    text(i,p_i(i)+0.01,num2str(p_i(i)),'VerticalAlignment','bottom','HorizontalAlignment','center');% 就是用test加数值，这个0.03看情况定，根据数值大小，再改就好了
end
 xlabel '场景'
 ylabel '概率'
 title '各场景概率结果'

1.3 其他函数

function ICDF_value = ICDF(x,mu,sigma)
%ICDF 此处显示有关此函数的摘要
ICDF_value = norminv(x,mu,sigma);
end

 function OrderVector= Order_Vector(vector)
       array=sort(vector,'ascend');    %  获得排序后的数组
       OrderVector=zeros(length(vector),1);
     for i=1:length(vector)
          temp=find(array==vector(i));   %find()函数，找寻元素在排序后的数组中的位置，有可能不止一个，因为可能有重复的元素
         if   length(temp)~=1               %若有重复的元素
              OrderVector(i)=temp(1);    %重复元素的赋予相同的秩
         else
               OrderVector(i)=temp;
        end
    end
end
% 数据
Load_forcast = [100	234	345	235	67	456	343	235	656	564	765	453	323	345	76	453	231	242	125	731	235	573	231	435];