引言
这里介绍一篇在二维的基础上复现三维Copula的应用。主要内容如下：
新能源的随机性、波动性及间歇性为电力系统规划带来困扰，对风、光出力和负荷的变化规律进行合理刻画，生成典型出力场景是新能源规划的常用方法。针对具有相关性的风、光和负荷出力典型场景难以生成的问题，本文首先得到风、光和负荷的最优边缘分布估计表达式，然后建立多种基于 Copula 函数的风、光和负荷电场出力联合分布模型，判断各个模型的拟合优度，选取最优 Copula 函数作为风电、光伏和负荷联合概率分布，最后采用最优 Copula 联合概率分布生成出力场景。算例分析表明，所得的出力场景符合其相关性，在反映某地区风光实际出力时有更高的准确性，可为电力系统可靠性分析和电网规划提供参考。
选取最优 Copula 函数作为风电、光伏和负荷联合概率分布，最后采用最优 Copula 联合概率分布生成风、光和负荷联合概率下的出力场景。
代码
clear all
clc
close all
%function fdm
%% 参数初始值
n_scenario = 200; % 场景数
n_reduction = 10; % 剩余场景数
data=xlsread('G:\数据1.xlsx','8月','I2:K290');
wf1=data(:,1);
wf2=data(:,2);
wf2( wf2 == 0 ) = 1e-4;
[D_U1, PD_U1]=marginfitdist(wf1);
EP1 = cdf(PD_U1{1},wf1);
IEP1 = icdf(PD_U1{1},EP1);
[D_U2, PD_U2]=marginfitdist(wf2);
EP2 = cdf(PD_U2{1},wf2);
IEP2 = icdf(PD_U2{1},EP2);
[D_U3, PD_U3]=marginfitdist(data(:,3));
EP3 = cdf(PD_U3{1},data(:,3));
IEP3 = icdf(PD_U3{1},EP3);
datacdf=[EP1 EP2 EP3];
datacdf( datacdf == 0 ) = 1e-4;
[A,fam,theta] = mucopulaselect(datacdf);
P_wt = ones(24, n_scenario);%风生成
P_pv = ones(24, n_scenario);%光生成
P_sum = ones(24, n_scenario);%负荷生成
% %%
% %生成200个场景
for t = 1:200
UVW = simrvine(24,A,fam,theta);
P_wt(:,t) = icdf(PD_U1{1},UVW(:,1));
P_pv(:,t) = icdf(PD_U2{1},UVW(:,2));
P_sum(:,t) = icdf(PD_U3{1},UVW(:,3));
end
figure(1)
% %这里以下代表三维的概率copulas结果出图，可以根据自己的需要计算%%
len=50;
example=zeros(len,len,len);
Pro=1/len:1/len:1;%表示边缘分布概率区间为0.1到1
for i=1:len
for j=1:len
for k=1:len
newdata=[Pro(i);Pro(j);Pro(k)]';
cdfvalue = vinecdf(newdata,A,fam,theta);
example(i,j,k)=cdfvalue;
end
end
end
save example
%
X = linspace(0.1,1,len); %X轴范围
Y = linspace(0.1,1,len); %Y轴范围
Z = linspace(0.1,1,len); %Z轴范围
%
[X1,Y1,Z1] = meshgrid(X,Y,Z);
xs = 0.2:0.2:0.8; %x轴切线
ys = 0.2:0.2:0.8; %y轴切线
zs =0.2:0.2:0.8; %z轴切线
S=example(1:len,1:len,1:len);
h = slice(X1,Y1,Z1,S,xs,ys,zs); %S为该坐标的值
shading flat
%shading interp
colormap('jet')
colorbar
set(gca,'LineWidth',2,'fontname','Times New Roman','FontSize',20,'FontWeight','bold');
disp('ok')
%% 场景生成图
figure()
subplot(2,3,1)
plot(P_wt, '--')
hold on
l2 = xlabel('t/h');
set(l2, 'Fontname', 'Times New Roman', 'FontSize', 50)
l3 = ylabel('P/kW');
set(l3, 'Fontname', 'Times New Roman', 'FontSize', 50)
set(gca, 'FontName', 'Times New Roman', 'FontSize', 50)
subplot(2,3,4)
plot(P_pv, '-')
hold on
l2 = xlabel('t/h');
set(l2, 'Fontname', 'Times New Roman', 'FontSize', 50)
l3 = ylabel('P/kW');
set(l3, 'Fontname', 'Times New Roman', 'FontSize', 50)
set(gca, 'FontName', 'Times New Roman', 'FontSize', 50)
%%
%原理：确定初始场景集合的一个子集，并给其重新分配场景概率，使保留场景的概率分布Q与初始场景集合的概率P之间的某种概率距离最短（即，P与Q相近），
%从而削减概率小的概率，将其加到与其场景的概率距离最近的场景上。
%%
%计算各个场景之间的概率距离
zdistance = zeros(n_scenario, n_scenario);
for i = 1:n_scenario
for j = 1:n_scenario
if i == j
distance(i, j) = 0; % 距离
else
distance(i, j) = sqrt(sum((P_sum(:, i) - P_sum(:, j)).^2));
end
end
end
p = ones(n_scenario, 1)*(1/n_scenario); % 各场景初始概率
%% 寻找最小概率距离场景
b2 = [];
distance(distance == 0) = inf;
for n = 1:(n_scenario - n_reduction) % 削减1990次，保留10个概率最高场景
[mink, index] = min(distance, [], 2);%index每行最小坐标列 %mink 每行最小数值 % min(k1, [], 2) 求取每行的最小值； min(k1, [], 1)求取每列的最小值
%删去index2 行 %%min(mink.*p) 概率最低。。。被淘汰
[mink11, index2] = min(mink.*p);
b = index2;
%减少一个场景
distance(b, :) = [];
distance(:, b) = [];
b2 = [b2;b];
%新概率生成
a = index(index2);%与被削减场景的概率距离最近的场景a
%新场景概率a=原来对应场景概率a+概率重新分配系数*与此情景概率距离最近场景index2
p(a) = p(index2) + p(a);
%一次循环后新的概率和场景
p(b) = [];
P_sum(:, b) = [];
P_wt(:, b) = [];
P_pv(:, b) = [];
%%
end%%%%一轮循环结束，场景削减1个。
%% 削减后的场景
subplot(2,3,2)
plot(wf1, '*')
hold on
plot(P_wt);
grid on
l2 = xlabel('t/h');
set(l2, 'Fontname', 'Times New Roman', 'FontSize', 50)
l3 = ylabel('P/kW');
set(l3, 'Fontname', 'Times New Roman', 'FontSize', 50)
set(gca, 'FontName', 'Times New Roman', 'FontSize', 50)
subplot(2,3,5)
plot(wf2, '*')
hold on
plot(P_pv);
grid on
l2 = xlabel('t/h');
set(l2, 'Fontname', 'Times New Roman', 'FontSize', 50)
l3 = ylabel('P/kW');
set(l3, 'Fontname', 'Times New Roman', 'FontSize', 50)
set(gca, 'FontName', 'Times New Roman', 'FontSize', 50)
%% 削减后的场景三维图
subplot(2,3,3)
bar3(P_wt')
l2 = xlabel('t/h');
set(l2, 'Fontname', 'Times New Roman', 'FontSize', 10)
ylabel('场景编号');
zlabel('风电出力');
subplot(2,3,6)
bar3(P_pv')
l2 = xlabel('t/h');
set(l2, 'Fontname', 'Times New Roman', 'FontSize', 10)
ylabel('场景编号');
zlabel('光伏出力');
%% 各个场景的概率
figure(3)
bar(p)
ylim([0, 0.30]);
xlabel('场景编号');
ylabel('概率');
set(gca, 'FontSize', 20)
no = 1;
figure(4)
[ss,gg]=meshgrid(1:n_reduction,1:24);
plot3(ss,gg,P_wt, 'linewidth', 2);
title(['考虑相关性生成的风电出力', num2str(n_reduction), '个场景'])
xlabel('场景'); ylabel('时刻');zlabel('风电出力值');
set(gca, 'FontSize', 20)
set(gca,'LineWidth',2);
figure(5)
[ss,gg]=meshgrid(1:n_reduction,1:24);
plot3(ss,gg,P_sum, 'linewidth', 2);
title(['考虑相关性生成的负荷', num2str(n_reduction), '个场景'])
xlabel('场景'); ylabel('时刻');zlabel('负荷值');
set(gca,'FontSize', 20)
set(gca,'LineWidth',2);
figure(6)
subplot(1,3,1)
for i=1:10
hold on
h=cdfplot(P_wt(:,i));
set(h,'LineStyle', '-', 'LineWidth',2)
end
title('');
grid off;
%h1=legend({'1','2','3','4','5','6','7','8','9','10'},'FontSize',16);
%xlabel('Value','FontSize',16,'fontname','Times New Roman');
%legend('Location', 'Best');
%set(h1, 'Box', 'off');
ylabel('Cumulative probability','FontSize',16,'fontname','Times New Roman');
set(gca,'fontname','Times New Roman','FontWeight','bold','FontSize',20);
set(gca,'LineWidth',1.5);
subplot(1,3,2)
for i=1:n_reduction
hold on
h=cdfplot(P_pv(:,i));
set(h,'LineStyle', '-', 'LineWidth',2)
end
title('');
grid off;
%h1=legend({'1','2','3','4','5','6','7','8','9','10'},'FontSize',16);
%xlabel('Value','FontSize',16,'fontname','Times New Roman');
%legend('Location', 'Best');
%set(h1, 'Box', 'off');
%ylabel('Cumulative probability','FontSize',16,'fontname','Times New Roman');
set(gca,'fontname','Times New Roman','FontWeight','bold','FontSize',20);
set(gca,'LineWidth',1.5);
subplot(1,3,3)
for i=1:10
hold on
h=cdfplot(P_sum(:,i));
set(h,'LineStyle', '-', 'LineWidth',2)
end
title('');
grid off;
h1=legend({'1','2','3','4','5','6','7','8','9','10'},'FontSize',16);
%xlabel('Value','FontSize',16,'fontname','Times New Roman');
legend('Location', 'Best');
set(h1, 'Box', 'off');
%ylabel('Cumulative probability','FontSize',16,'fontname','Times New Roman');
set(gca,'fontname','Times New Roman','FontWeight','bold','FontSize',20);
set(gca,'LineWidth',1.5);
效果
如果需要更多关于Copula函数应用的联合重现期、同现重现期和条件重现期等的代码（MATLAB和R语言）（代码可直接运行，替换数据即可运行出结果）如需请关注公众号回复 copula
欢迎关注公众号，获取更多科研代码和前沿论文资讯等相关内容
参考文献
[1]宋宇,李涵.基于核密度估计和Copula函数的风、光出力场景生成[J].电气技术,2022,23(01):56-63.
[2]林顺富,刘持涛,李东东,等.考虑电能交互的冷热电区域多微网系统双层多场景协同优化配置[J].中国电机工程学报,2020,40(05):1409-1421.DOI:10.13334/j.0258-8013.pcsee.190275.
[3]白凯峰,顾洁,彭虹桥,等.融合风光出力场景生成的多能互补微网系统优化配置[J].电力系统自动化,2018,42(15):133-141.
576
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
