【物理应用】非序贯蒙特卡洛法评估风电系统【含matlab源码 766期】.zip

clc clear %% 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 % load data % mu=mean(speed);%原始数据的统计参数 % sigma=sqrt(var(speed)); % % % 2.1数据标准化 % Stdspeed=(speed-mu)/sigma; % % Stdspeed=iddata(Stdspeed); % AIC=zeros(10,1); % % for n=2:10 % for n=2:7 % sys=armax(Stdspeed,[n,n-1]); % NoiseStd=sqrt(sys.NoiseVariance); % e=normrnd(0,NoiseStd,8760,1);%产生白噪声序列[时序标号 对应时序的值] % y=zeros(8760,1);%预测后的标准风速[时序标号 对应时序的值] % % % % 2.2利用arma模型产生预测值,1-8760项 % y(1)=e(1);%第一项等于第一个噪声 % for i=2:n %2-n项通项公式一致 % y(i)=-sys.A(2:i)*y(i-1:-1:1)+sys.C(1:i)*e(i:-1:1); % end % % for i=n+1:8760% n+1至8760项通项公式一致 % y(i)=-sys.A(2:n+1)*y(i-1:-1:i-n)+sys.C(1:n)*e(i:-1:i-n+1); % end % % % 2.3计算AIC值 % % 计算残差 % s=0; % for i=1:8760 % s=s+(Stdspeed(i)-y(i))^2; % end % % AIC(n,1)=8760*log(s)+2*n; % end % [minAIC,optn]=min(AIC(2:7));%找到AIC最小的阶数 % % % 2.3得到最优阶数，代入arma模型预测风速 % sys=armax(Stdspeed,[optn,optn-1]); % Noise.Std=sqrt(sys.NoiseVariance); % e=normrnd(0,Noise.Std,8760,1);%产生白噪声序列[时序标号 对应时序的值] % y=zeros(8760,1);%预测后的标准风速[时序标号 对应时序的值] % % % % y(1)=e(1);%第一项等于第一个噪声 % for i=2:optn %2-n项通项公式一致 % y(i)=-sys.A(2:i)*y(i-1:-1:1)+sys.C(1:i)*e(i:-1:1); % end % % for i=optn+1:8760% n+1至8760项通项公式一致 % y(i)=-sys.A(2:optn+1)*y(i-1:-1:i-optn)+sys.C(1:optn)*e(i:-1:i-optn+1); % end % % 得到预测值y后反变换为实际风速 % Simspeed=y*sigma+mu; % % % 计算arma预测风速的概率分布 % [count,x]=hist(Simspeed,ceil(max(Simspeed)/0.5));%x是区间中心数,组距0.5 % prob3=count/8760/0.5;%计算arma预测数据概率密度 ，频数除以数据种数，除以组距 % % % 时序作图比较 % figure(4) % plot(1:400,speed(1:400),'r-.',1:400,Simspeed(1:400),'b-'); % legend('实际风速','ARMA模拟风速'); % % 概率分布作图比较 % figure(5) % bar(x,prob3,1) % title('ARMA预测概率分布'); % % save('result_arma.mat') %% 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） % 3.1 求出常规机组的出力模型，按类构成多状态模型 % RBTS发电系统中共有6类常规机组，%11台常规机组数据 % %2台5MW水电机组%% %1台10MW热电机组%% %4台20MW水电机组%% %1台20MW热电机组%% %1台40MW水电机组%% %2台40MW热电机组% % Generator.Norm(数量 单机容量 FOR) % Generator.Norm=[2 5 0.01 % 1 10 0.02 % 4 20 0.015 % 1 20 0.025 % 1 40 0.02 % 2 40 0.03]; % % %StateNorm中每一个元胞表示一个类机组的多状态 % StateNorm=cell(1,length(Generator.Norm)); % State=[]; % for i=1:length(Generator.Norm) % for j=1:Generator.Norm(i,1)+1 % % State[出力 概率]指某类机组多状态的二项分布概率 % State(j,:)=[(j-1)*Generator.Norm(i,2) binopdf(j-1,Generator.Norm(i,1),1-Generator.Norm(i,3))]; % end % StateNorm{i}=State; % end % Gener