风电出力的不确定性主要源于预测误差，而研究表明预测误差（e）服从正态分布且大概为预测出力的10%。本代码采用拉丁超立方抽样实现场景生成[1,2]、基于概率距离的快速前代消除法实现场景缩减[3]，以此模拟了风电出力的不确定性。

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 = 5;                                      % 代表性场景数量 5个
load('Pw')
step = 24;                                  % 步长 h（总计一天）
P_test    = Pw';                            % 转置-待生成数据
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
    for j=1:iterations
        Value(i,j) = ICDF……

2.2 场景缩减

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

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+10000*eye(k);……

完整代码请关注作者并私聊获取

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.