%% 2.计及风速和元件故障的风电场出力

% 2.1得到N台机组M状态的风电场出力模型

% 风速单位m/s ,切出功率单位MW，相关出力数据来自威布尔分布预测的风速，相关计算查看百度文库“风速威布尔分布和ARMA模型预测matlab程序”

clc

clear

loadresult_weibull

Generator.Wind=struct('vin',3,'vout',25,'vr',15,'Pr',1.5,'FOR',0.028,'lamda',5,'mu',175.2);

N=10;

M=6;

% 风电功率转换函数

[~,k1,k2]=Power([],Generator);

StateWeibull=zeros(M,2);

fori=1:M

StateWeibull(i,1)=(i-1)/(M-1)*Generator.Wind.Pr;

ifStateWeibull(i,1)==0

% StateWeibull(i,2)=1-exp(-(Vin/c)^k)+exp(-(Vout/c)^k))+

StateWeibull(i,2)=wblcdf((((2*i-1)/(2*M-2))*Generator.Wind.Pr-k2)/k1,c,k)+exp(-(Generator.Wind .vout/c)^k) ;

end

ifStateWeibull(i,1)>0&&StateWeibull(i,1)<Generator.Wind.Pr

StateWeibull(i,2)=wblcdf((((2*i-1)/(2*M-2))*Generator.Wind.Pr-k2)/k1,c,k)-wblcdf((((2*i-3)/(2*M -2))*Generator.Wind.Pr-k2)/k1,c,k);

end

ifStateWeibull(i,1)==Generator.Wind.Pr

StateWeibull(i,2)=wblcdf((Generator.Wind.Pr-k2)/k1,c,k)-wblcdf((((2*i-3)/(2*M-2))*Generator.Wi nd.Pr-k2)/k1,c,k);

end

end

% StateWeibull表示单台机组的6状态出力，StateWeibull6表示N台风机构成的风电场的6状态出力

StateWeibull6=[StateWeibull(:,1)*N StateWeibull(:,2)];

figure(2)

bar(StateWeibull6(:,1),StateWeibull6(:,2))

grid on

title('只计及风速Weibull分布的风电场6状态出力模型')

% 2.2计算N台风电机组运行状态

% StateRun=binornd(N,Generator.Wind.FOR,N,1);

StateFOR=[(0:N)' zeros(N+1,1)];

fori=1:length(StateFOR)

StateFOR(i,2)=binopdf(i-1,N,1-Generator.Wind.FOR);

end

% 2.3计算考虑FOR的风电场M状态出力模型（Weibull）StateFORWeibull6=zeros(M,2);

fori=1:M

StateFORWeibull6(i,1)=(i-1)/(M-1)*Generator.Wind.Pr*N;

end

state=StateFOR(:,1)*StateWeibull(:,1)';

prob=StateFOR(:,2)*StateWeibull(:,2)';

[n,m]=size(state);

fori=1:n

for j=1:m

for k=1:M+1

if abs(state(i,j)-StateFORWeibull6(k,1))<=1/(2*M-2)*Generator.Wind.Pr*N StateFORWeibull6(k,2)=StateFORWeibull6(k,2)+prob(i,j);

break

else

continue

end

end

end

end

figure(3)

bar(StateFORWeibull6(:,1),StateFORWeibull6(:,2));

grid on

title('计及元件随机故障的风电场6状态出力模型（Weibull分布）'); save('result_WindFarmOutput.mat')