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

光伏出力的不确定性主要源于预测误差，而研究表明预测误差（e）服从正态分布且大概为预测出力的10%。本代码采用拉丁超立方抽样实现场景生成[1,2]、基于概率距离的快速前代消除法实现场景缩减[3]，以此模拟了光伏出力的不确定性。与风电不确定性模拟不同之处在于——光伏存在0出力点，在不确定模拟中容易出现超出索引范围的情况，本代码解决了该问题。
复现论文
[1] SCI-Clustering based unit commitment with wind power uncertainty. Energy Conversion and Management, 2016, 111:89-102.
[2] 中文核心-含风光水的虚拟电厂与配电公司协调调度模型[J].电力系统自动化,2015,39(09):75-81+207.
关注专栏可查看Matlab全部代码

1、Matlab代码

1.1 场景生成

% 不确定性模拟之场景生成 + 场景缩减
clc 
clear
%% 参数设置
Z = 24;                                     % 以一天24h为周期进行场景生成
iterations = 1000;                          % 抽样次数 1000次-1000种场景
S = 10;                                      % 代表性场景数量 5个
load('PV_forecast.mat')
step = 24;                                  % 步长 h（总计一天）
P_test    = PV_forecast+0.001;               % 转置-待生成数据
mu        = 0*ones(Z,1);                    % 正态分布-均值
sigma     = 0.1*P_test';                    % 正态分布-标准差（此时取10%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 '各场景概率结果'

其他函数:

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

PV_forecast = [0	0	0	0	0	0	0	499	1538	3101	5467	5840	5669	5910	5690	5138	3354	128	0	0	0	0	0	0];

2、结果输出

2.1 1000种生成的场景集

2.2 缩减后的10种场景

2.3 10种场景对应概率