论文复现：模拟风电不确定性——拉丁超立方抽样生成及缩减场景（Matlab全代码）

WW、forever

已于 2023-04-01 09:58:02 修改

阅读量3.8k

点赞数 10

分类专栏：论文复现：模拟风电、光伏和负荷不确定性文章标签： matlab 风电、光伏出力不确定性

于 2023-02-20 21:16:26 首次发布

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

版权

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

3 篇文章 60 订阅

订阅专栏

风电出力的不确定性主要源于预测误差，而研究表明预测误差（e）服从正态分布且大概为预测出力的10%。本代码采用拉丁超立方抽样实现场景生成[1,2]、基于概率距离的快速前代消除法实现场景缩减[3]，以此模拟了风电出力的不确定性。
复现论文
[1] SCI-Clustering based unit commitment with wind power uncertainty. Energy Conversion and Management, 2016, 111:89-102.
[2] 中文核心-基于拉丁超立方采样的含风电电力系统的概率可靠性评估[J].电工技术学报,2016,31(10):193-206.
[3] 中文核心-含风光水的虚拟电厂与配电公司协调调度模型[J].电力系统自动化,2015,39(09):75-81+207.
关注专栏可查看Matlab全部代码
在这里插入图片描述

1 风电不确定性模拟理论

1.1 不确定性模拟

不确定性模拟大致可以分为随机优化和场景分析两类。

场景分析（scenario analysis）是一种通过构建确定性场景来分析电力系统不确定性问题的方式，它是解决含可再生能源的电力系统优化规划运行问题的一种有效途径。
基于抽样的场景生成方法：通过对概率分布进行统计学抽样并将输出样本作为生成的场景来得到离散场景，包括蒙特卡罗（Monte Carlo, MC）方法、**拉丁超立方体采样（Latin Hypercube Sampling, LHS）**方法和马尔可夫链蒙特卡洛（Markov chain Monte Carlo, MCMC）等等。具体抽样方法根据初始数据相关性、需求不同而采用不同方法。
蒙特卡罗与拉丁超立方的区别：蒙特卡罗模拟方法由于采用简单随机抽样，在不确定性模拟中需要高计算时间和大量计算机存储。而拉丁超立方抽样是一种随机分层抽样方法，不需要大规模储存和长时间计算，模拟样本与传统蒙特卡罗方法相比，能更好地反映变量的分布范围。

1.2 风电不确定性模拟

将风电的不确定性转换为风电预测误差的不确定性。
在这里插入图片描述
e（预测误差）服从标准正态分布，u=0，核心在sigma。统计学家已发现sigma=5~10%*Pforecast 。

2 Matlab代码

2.1 场景生成

拉丁超立方抽样法是一种分层抽样方法，已广泛应用于模拟风电预测误差，基于拉丁超立方抽样在95%置信区间范围内生成风电的1000组场景：
在这里插入图片描述

Matlab代码：

% 不确定性模拟之场景生成 + 场景缩减
clc 
clear
%% 参数设置
Z = 24;                                     % 以一天24h为周期进行场景生成
iterations = 1000;                          % 抽样次数 1000次-1000种场景
S = 10;                                      % 代表性场景数量 5个
load('Pw_forecast.mat')
% 风电预测数据
Pw_forecast=[110 147	151	130	112	99	103	111	102	256	479	561	573	566	528	455	311	240	200	169	136	102	73	52];
step = 24;                                  % 步长 h（总计一天）
P_test    = Pw_forcast;                     % 转置-待生成数据
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)

图片输出：
在这里插入图片描述

2.2 场景缩减

由于抽样所产生的场景较多，为提高计算效率，需要减少并挑选代表性场景。基于概率距离的快速前代消除法将1000组场景消减为具有代表性的10组：

Matlab代码：

%% 场景削减
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 '各场景概率结果'
Outputw1_1000=a_thou;
Outputw1_S=a_five;
Outputw1_p=a_pro;

图片输出：
在这里插入图片描述

函数

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