含风电场电力系统的可靠性评估Matlab程序

clcclearyear=50%模拟的年限forl=1:year;%%%%%%%%%%%%%%%%%%%%%%%%%第一步：数据导入与预处理%%%%%%%%%%%%%%%%%%%%SW0=load('windspeed.txt');%载入原始风速数据SW0=xlsread('windspeed.xls');SW0=SW0';SW0=SW0/10*3.6;%原始数据的风速单位为0.1m/s，这里转化为km/hN=size(SW0,2);mu=mean(SW0);sigma=var(SW0);sigma=sigma^0.5;%求样本的平均值和标准差y=(SW0-mu)./sigma;%数据预处理%figure(1);%subplot(211);%autocorr(y);%画出自相关图%title('自相关图');%subplot(212);%parcorr(y);%画出偏自相关图%title('偏相关图');%%%%%%%%%%%%%%第二步：根据AIC准则确定ARMA模型的阶数%%%%%%%%%%%%%%%%%%forn=2:7;m=armax(y',[n,n-1]);fai=-m.a;theta=m.c;%把armax函数得到的参数，取出来fori=1:n;y1(i)=y(i);endNoise=m.NoiseVariance^0.5;e=normrnd(0,Noise,1,N);fort=n+1:1:N;y1(t)=0;forj=2:n+1;y1(t)=y1(t)+fai(j)*y1(t-(j-1));endfork=1:n;y1(t)=y1(t)+theta(k)*e(t-(k-1));endend%y1(t)为预测值s(n)=0;fori1=1:N;residual=y1(i1)-y(i1);s(n)=s(n)+residual^2;%求取残差平方和endAIC(n)=N*log(s(n))+2*n-1;%求AICendarma=AIC(1,2:7);%n=2-7时，各AIC的值[AIC,n1]=min(arma);n1=n1+1;%n1为得到的ARMA模型的阶数%%%%%%%%%%%%%%%%%第三步：用ARMA模型预测风速并确定风速分布%%%%%%%%%%%%%%%%%m=armax(y',[n1,n1-1]);fai=-m.a;theta=m.c;fori=1:n1;y2(i)=y(i);endNoise=m.NoiseVariance^0.5;e=normrnd(0,Noise,1,8736);fort=n1+1:1:8736;y2(t)=0;forj=2:n1+1;y2(t)=y2(t)+fai(j)*y2(t-(j-1));endfork=1:n1;y2(t)=y2(t)+theta(k)*e(t-(k-1));endend%y1(t)为预测值y2;SW=y2.*sigma+mu;%figure(2)%subplot(121);%hist(SW0,100);%xlabel('风速');ylabel('频数');title('原始风速的分布');%subplot(122);%hist(SW,100);%xlabel('风速');ylabel('频数');title('预测风速的分布');%得到8736个小时的预测风速：SW1*8736%%%%%%%第四步：风电场的转移过程，确定其一年中三种状态分别存在的时长%%%%%%%%风电机组的三状态模型lambdaRD=5.84;lambdaRF=7.96;lambdaDR=48.3;lambdaFR=58.4;lambdaDF=0;lambdaFD=0;%风力发电机的3个状态的转移率T=[111;lambdaRD-lambdaDR-lambdaDFlambdaFD;lambdaRFlambdaDF-lambdaFR-lambdaFD];T=inv(T);probWTG=T*[1;0;0];%风力发电机分别处于运行、降额以及停运状态的概率TR=8760/(lambdaRD+lambdaRF);TD=8760/(lambdaDR+lambdaDF);TF=8760/(lambdaFD+lambdaFR);WTGnum=25;R1=rand(WTGnum,1000);%R1确定风力发电机所处的状态R2=rand(WTGnum,1000);%R2确定该状态所持续的时间D=zeros(WTGnum,1000);%D记录每个状态所处的时间time=zeros(WTGnum,1001);%time记录每个状态变化的时间节点alpha1=zeros(WTGnum,1000);%alpha1记录每个状态分别是什么N1=8736;%N1为模拟的小时数fori2=1:WTGnum%i2表示有多少台风力发电机fori3=1:1000;ifR1(i2,i3)probWTG(1,1)alpha1(i2,i3)=1;D(i2,i3)=-TR*log(R2(i2,i3));%alpha记录第i2台发电机，在第i3个时刻的状态，D是相应的持续时长elseifR1(i2,i3)probWTG(1,1)+probWTG(2,1);alpha1(i2,i3)=0.8;D(i2,i3)=-TD*log(R2(i2,i3));elsealpha1(i2,i3)=0;D(i2,i3)=-TF*log(R2(i2,i3));endtime(i2,i3+1)=time(i2,i3)+D(i2,i3);iftime(i2,i3+1)=N1;breakendendr(i2)=i3;%r记录了第i2台发电机状态的总变化次数end%记录每次风力发电机所处的状态，以及该状态存在的时长，达到1年后停止循环%a=time(1,1:r(1)+1);%b=[alpha1(1,1:r(1)),ones(1,1)];%a1=time(2,1:r(2)+1);%b1=[alpha1(2,1:r(2)),ones(1,1)];%figure(3)%subplot(211)%stairs(a,b,'LineWidth',2.0);%作出风力发电机输出功率状态图%title('一年中风力发电机的输出状态：正常、降额、停运');%xlabel('时间');ylabel('状态');%axis([0N101.2]);%gridon;%subplot(212)%stairs(a1,b1,'LineWidth',2.0);%xlabel('时间');ylabel('状态');%axis([0N101.2]);%gridon;fori2=1:WTGnum;fori4=1:N1;%N1=8736forj2=1:r(i2)+1;ifi4time(i2,j2+1)alpha2(i2,i4)=alpha1(i2,j2);breakelseendendendend%得到一年中25台风力发电机8736小时的每个小时所处的停运状态：alpha225*8736%%%%%%%%%%%%%%%第五步：考虑尾流效应以及地形对风速的影响，并计算功率%%%%%%%%%%%%beta1=struct('landformfactor',{111.10781.29080.7747},'wakeeffect',{10.92470.92950.94050.9046});%尾流+地形%beta1=struct('landformfactor',{11111},'wakeeffect',{11111});%都不考虑%beta1=struct('landformfactor',{11111},'wakeeffect',{10.92470.92470.92470.9247});%尾流fori6=1:N1;x=alpha2(:,i6);x=x';x=reshape(x,5,5)';flag=ones(5,5);fori7=2:5;forj3=1:5;ifx(i7-1,j3)==0;flag(i7,j3)=beta1(i7).landformfactor*flag(i7-1,j3);elseflag(i7,j3)=beta1(i7).landformfactor*beta1(i7).wakeeffect*flag(i7-1,j3);endendendflag=reshape(flag',1,25);flag=flag';flag1(:,i6)=flag;end%得到尾流效应以及地形对风速的影响：flag125*8736%%%%%%%%%%%%%%%%%第六步：根据预测得到的风速计算风力发电机发出的功率%%%%%%%%%%%%%%%%%Vci=14.4;Vr=45;Vco=90;Pr=1.5;%风力发电机切入、额定、切除风速km/h，额定功率1.5MWalpha=1/((Vci-Vr)^2);beta=((Vci+Vr)/(2*Vr))^3;A=alpha*(Vci*(Vci+Vr)-4*Vci*Vr*beta);B=alpha*(4*(Vci+Vr)*beta-(3*Vci+Vr));C=alpha*(2-4*beta);SW=[SW;SW;SW;SW;SW];SW=[SW;SW;SW;SW;SW];SW1=SW.*flag1;forj4=1:WTGnum;forj1=1:N1;ifSW1(j4,j1)VciP(j4,j1)=0;elseifSW1(j4,j1)VrP(j4,j1)=Pr*(A+B*SW1(j4,j1)+C*SW1(j4,j1)*SW1(j4,j1));elseifSW1(j4,j1)VcoP(j4,j1)=Pr;elseP(j4,j1)=0;endendendP1=P.*alpha2;%得到25台风力发电机8736小时的每小时的输出功率：P125*8736%figure(4)%subplot(121);%hist(P1(1,:),10);%xlabel('功率');%title('考虑尾流、地形的功率分布-1号风机');%subplot(122)%hist(P1(2,:),10);%xlabel('功率');%title('考虑尾流、地形的功率分布-2号风机2');pWTG=zeros(1,N1);fori8=1:N1;fori9=1:WTGnum;pWTG(i8)=pWTG(i8)+P1(i9,i8);endendpWTG=pWTG*15;%得到风力发电机每小时输出的总功率：pWTG1*8736%%%%%%%%%%%%%%%%%%%%第七步：计算常规发电机发出的功率%%%%%%%%%%%%%%%%%FOR=load('forced_outage_rate.txt');FOR=FOR';pG=zeros(1,N1);forj5=1:N1;power=[12121212122020202050505050505076767676100100100155155155155197197197350400400];r1=rand(1,32);forj6=1:32;ifr1(j6)=FOR(j6);power(j6)=0;endpG(j5)=pG(j5)+power(j6);%第i2次模拟中常规机组输出的总功率pG(i2)endend%得到常规发电机每小时的输出总功率：pG1*8736loads=load('load.txt');loads1=loads(:,2)'*2850;%导入负荷：loads1*8736%figure(5)%h=1:N1;%z=pWTG+pG;%plot(h,z,'+');%holdon%plot(h,loads1);%legend('输出总功率','负荷');%画总输出功率与负荷曲线%%%%%%%%%%%%%%%%%%%%第八步：计算LOLP和EENS%%%%%%%%%%%%%%%%%%%%%%%%%%X=0;W=0;fori10=1:N1;ifpWTG(i10)+pG(i10)loads1(i10);x1=1;w(i10)=(load