风电等值容量以及风电可信容量的matlab仿真,包含仿真操作录像,代码注释

clc; clear; close all; warning off; addpath(genpath(pwd)); rng('default') %% 1.计算风速weibull分布 % 数据处理 tic load data; mu=mean(speed);%原始数据的统计参数 sigma=sqrt(var(speed)); % 计算威布尔分布参数 parmhat=wblfit(speed); k=parmhat(2); c=parmhat(1); % k=(sigma/mu)^-1.086; % c=mu/gamma(1+1/k); % 威布尔分布拟合 [y,x]=hist(speed,ceil(max(speed)/0.5));%x是区间中心数,组距-1.5 prob1=y/8760/0.5;%计算原始数据概率密度 ，频数除以数据种数，除以组距 prob2=(k/c)*(x/c).^(k-1).*exp(-(x/c).^k);%威布尔分布 figure(1) title('Weibull分布拟合图'); bar(x,prob1,1) hold on plot(x,prob2,'r') legend('历史数据','Weibull拟合结果') % legend('Weibull拟合结果') hold off save('result_weibull.mat') toc % c=cumsum(prob2) % plot(x,c) %% 2.ARMA模型预测风速 clc clear tic load data y=speed(1:300); Data=y; %共300个数据 SourceData=Data(1:250,1); %前250个训练集 step=50; %后50个测试 TempData=SourceData; TempData=detrend(TempData);%去趋势线 TrendData=SourceData-TempData;%趋势函数 %--------差分，平稳化时间序列--------- H=adftest(TempData); difftime=0; SaveDiffData=[]; while ~H SaveDiffData=[SaveDiffData,TempData(1,1)]; TempData=diff(TempData);%差分，平稳化时间序列 difftime=difftime+1;%差分次数 H=adftest(TempData);%adf检验，判断时间序列是否平稳化 end %---------模型定阶或识别-------------- u = iddata(TempData); test = []; for p=1:5 %自回归对应PACF,给定滞后长度上限p和q，一般取为T/10、ln(T)或T^(1/2),这里取T/10=12 for q=1:5 %移动平均对应ACF m = armax(u,[p q]); AIC = aic(m); %armax(p,q),计算AIC test = [test;p q AIC]; end end for k=1:size(test,1) if test(k,3) == min(test(:,3)) %选择AIC值最小的模型 p_test = test(k,1); q_test = test(k,2); break; end end %------1阶预测----------------- TempData=[TempData;zeros(step,1)]; n=iddata(TempData); %m = armax(u(1:ls),[p_test q_test]); %armax(p,q),[p_test q_test]对应AIC值最小，自动回归滑动平均模型 m = armax(u,[p_test q_test]); % ------------------------------------------- P1=predict(m,n,1); PreR=P1.OutputData; PreR=PreR'; Noise.std=sqrt(m.NoiseVariance); e=normrnd(0,Noise.std,1,300); for i=251:300 PreR(i)=-m.A(2:p_test+1)*PreR(i-1:-1:i-p_test)'+m.C(1:q_test+1)*e(i:-1:i-q_test)'; end % ------------------------------------------- %----------还原差分----------------- if size(SaveDiffData,2)~=0 for index=size(SaveDiffData,2):-1:1 PreR=cumsum([SaveDiffData(index),PreR]); end end %-------------------预测趋势并返回结果---------------- mp1=polyfit([1:size(TrendData',2)],TrendData',1); xt=[]; for j=1:step xt=[xt,size(TrendData',2)+j]; end TrendResult=polyval(mp1,xt); PreData=TrendResult+PreR(size(SourceData',2)+1:size(PreR,2)); tempx=[TrendData',TrendResult]+PreR; % tempx为预测结果 plot(tempx,'r-.'); hold on plot(Data,'b'); legend('ARMA拟合时序曲线','实际时序风速'); save('resultarma.mat'); toc %% 2.计及风速和元件故障的风电场出力 % 2.1得到N台机组M状态的风电场出力模型 % 风速单位m/s ,切出功率单位MW % clc clear load result_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); for i=1:M StateWeibull(i,1)=(i-1)/(M-1)*Generator.Wind.Pr; if StateWeibull(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 if StateWeibull(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 if StateWeibull(i,1)==Generator.Wind.Pr StateWeibull(i,2)=wblcdf((Generator.Wind.Pr-k2)/k1,c,k)-wblcdf((((2*i-3)/(2*M-2))*Generator.Wind.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)]; for i=1:length(StateFOR) StateFOR(i,2)=binopdf(i-1,N,1-Generator.Wind.FOR); end % 2.3计算考虑FOR的风电场M状态出力模型（Weibull） StateFORWeibull6=zeros(M,2); for i=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); for i=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') %% 3.计算含风电场的发电系统可靠性指标（非序贯MC） clc clear load result_WindFarmOutput % 3.1 求出常规机组的出力模型，按类构成多状态模型 % RBTS发电系统中共有6类常规机组，%11台常规机组数据 % %2台5MW水电机组%% %1台10MW热电机组%% %4台20MW水电机组%% %1台20MW热电机组%% %1台40MW水电机组%% %2台40MW热电机组% Generator.Norm=[5 0.01 5 0.01 10 0.02 20 0.015 20 0.015 20 0.015 20 0.015 20 0.025 40 0.02 40 0.03 40 0.03]; save('process.mat'); % 3.2MC抽样机组确定机组状态 % 3.2.1计算含风电场的RBTS可靠性 % 共有7类机组，常规机组状态在StateNorm【出力 概率】元胞数组中，风电状态在StateFORWeibull6【出力 概率】 I=0 %I用来记录发生却负荷的次数 sumDNS=0; DNS=zeros(200000,1); K=rand(200000,12);%1-11常规 12风电 pwind=zeros(200000,1); for k=1:200000 Pout=zeros(12,1); %得到一次抽样常规机组状态 for i=1:11 if K(k,i)>Generator.Norm(i,2) Pout(i)=Generator.Norm(i,1); else Pout(i)=0; end end %得到一次抽样风电场的出力 t=0; for i=1:M t=t+StateFORWeibull6(i,2); if K(k,12)-t>=0 continue else Pout(12)=StateFORWeibull6(i,1); pwind(k)= Pout(12); break end end sumPout=sum(Pout);%发电系统出力 % 统计系统状态参数 % 持续峰荷185MW RBTSload=185; DNS(k)=max([0;RBTSload-sumPout]); if DNS(k)>0 I=I+1; %I用来记录发生却负荷的次数 sumDNS=sumDNS+DNS(k); end end LOLP1=I/200000; EENS1=sumDNS/200000*8760; LOLE1=I/200000*8760; disp('LOLP1=');disp(LOLP1); disp('EENS1=');disp(EENS1); disp('LOLE1=');disp(LOLE1); save('resultwind.mat','LOLP1','EENS1','LOLE1'); %% 3.2.2计算常规RBTS系统可靠性 clc clear load process I=0; sumDNS=0; DNS=zeros(100000,1); K=rand(100000,11);%1-11常规 7风电 for k=1:100000 Pout=zeros(11,1); %得到一次抽样常规机组状态 for i=1:11 if K(k,i)>Generator.Norm(i,2) Pout(i)=Generator.Norm(i,1); else Pout(i)=0; end end sumPout=sum(Pout);%发电系统出力 % 统计系统状态参数 % 持续峰荷185MW RBTSload=185; DNS(k)=max([0;RBTSload-sumPout]); if DNS(k)>0 I=I+1; %I用来记录发生却负荷的次数 sumDNS=sumDNS+DNS(k); end end % 计算LOLP EENS LOLE LOLP2=I/100000; EENS2=sumDNS/100000*8760; LOLE2=I/100000*8760; disp('LOLP2=');disp(LOLP2); disp('EENS2=');disp(EENS2); disp('LOLE2=');disp(LOLE2); save('result0.mat','LOLP2','EENS2','LOLE2'); clear load process load resultwind load result0 % 计算EFC定义下风电等值容量 P0=0;%设定初值Pe下限为0 R0=EENS2; P=[];%设定变量P、R存储每次迭代后Pe和对应的可靠性指标 R=[]; % 计算初值P(1)的可靠性指标,P(1)对应等值容量上限，和风电容量相同 P(1)