clc
clear

year=50???%模擬的年限

for?l=1:year;
????
????%%%%%%%%%%%%%%%%%%%%%%%%%??第一步：數(shù)據(jù)導(dǎo)入與預(yù)處理??%%%%%%%%%%%%%%%%%%%
????
????%SW0=load（‘windspeed.txt‘）;????%載入原始風(fēng)速數(shù)據(jù)
????SW0=xlsread（‘windspeed.xls‘）;
????SW0=SW0‘;
????SW0=SW0/10*3.6;????%原始數(shù)據(jù)的風(fēng)速單位為0.1m/s，這里轉(zhuǎn)化為km/h
????N=size（SW02）;
????mu=mean（SW0）;
????sigma=var（SW0）;
????sigma=sigma^0.5;????%求樣本的平均值和標(biāo)準(zhǔn)差
????y=（SW0-mu）./sigma;????%數(shù)據(jù)預(yù)處理
????
????%?figure（1）;
????%?subplot（211）;
????%?autocorr（y）;?????%畫出自相關(guān)圖
????%?title（‘自相關(guān)圖‘）;
????%?subplot（212）;
????%?parcorr（y）;????%畫出偏自相關(guān)圖
????%?title（‘偏相關(guān)圖‘）;
????
????
????%%%%%%%%%%%%%%??第二步：根據(jù)AIC準(zhǔn)則確定ARMA模型的階數(shù)??%%%%%%%%%%%%%%%%%%
????
????for?n=2:7;
????????m=armax（y‘[nn-1]）;
????????fai=-m.a;
????????theta=m.c;????%把a(bǔ)rmax函數(shù)得到的參數(shù)，取出來(lái)
????????
????????for?i=1:n;
????????????y1（i）=y（i）;
????????end
????????
????????Noise=m.NoiseVariance^0.5;
????????e=normrnd（0Noise1N）;
????????for?t=n+1:1:N;
????????????y1（t）=0;
????????????for?j=2:n+1;
????????????????y1（t）=y1（t）+fai（j）*y1（t-（j-1））;
????????????end
????????????for?k=1:n;
????????????????y1（t）=y1（t）+theta（k）*e（t-（k-1））;
????????????end
????????end????%y1（t）為預(yù)測(cè)值
????????
????????s（n）=0;
????????for?i1=1:N;
????????????residual=y1（i1）-y（i1）;
????????????s（n）=s（n）+residual^2;????%求取殘差平方和
????????end
????????
????????AIC（n）=N*log（s（n））+2*n-1;????%求AIC
????end
????
????arma=AIC（12:7）;????%n=2-7時(shí)，各AIC的值
????[AICn1]=min（arma）;
????n1=n1+1;????%n1為得到的ARMA模型的階數(shù)
????
????%%%%%%%%%%%%%%%%%??第三步：用ARMA模型預(yù)測(cè)風(fēng)速并確定風(fēng)速分布??%%%%%%%%%%%%%%%%%
????
????m=armax（y‘[n1n1-1]）;
????fai=-m.a;
????theta=m.c;
????
????for?i=1:n1;
????????y2（i）=y（i）;
????end
????
????Noise=m.NoiseVariance^0.5;
????e=normrnd（0Noise18736）;
????for?t=n1+1:1:8736;
????????y2（t）=0;
????????for?j=2:n1+1;
????????????y2（t）=y2（t）+fai（j）*y2（t-（j-1））;
????????end
????????for?k=1:n1;
????????????y2（t）=y2（t）+theta（k）*e（t-（k-1））;
????????end
????end????%y1（t）為預(yù)測(cè)值
????y2;
????SW=y2.*sigma+mu;
????
????%?????figure（2）
????%?????subplot（121）;
????%?????hist（SW0100）;
????%?????xlabel（‘風(fēng)速‘）;ylabel（‘頻數(shù)‘）;title（‘原始風(fēng)速的分布‘）;
????%?????subplot（122）;
????%?????hist（SW100）;
????%?????xlabel（‘風(fēng)速‘）;ylabel（‘頻數(shù)‘）;title（‘預(yù)測(cè)風(fēng)速的分布‘）;
????
????
????%得到8736個(gè)小時(shí)的預(yù)測(cè)風(fēng)速：SW?1*8736
????
????
????%%%%%%%??第四步：風(fēng)電場(chǎng)的轉(zhuǎn)移過(guò)程，確定其一年中三種狀態(tài)分別存在的時(shí)長(zhǎng)%%%%%%%
????
????%風(fēng)電機(jī)組的三狀態(tài)模型
????lambdaRD=5.84;lambdaRF=7.96;lambdaDR=48.3;lambdaFR=58.4;lambdaDF=0;lambdaFD=0;????%風(fēng)力發(fā)電機(jī)的3個(gè)狀態(tài)的轉(zhuǎn)移率
????T=[1?1?1;lambdaRD?-lambdaDR-lambdaDF?lambdaFD;lambdaRF?lambdaDF?-lambdaFR-lambdaFD];
????T=inv（T）;
????probWTG=T*[1;0;0];????%風(fēng)力發(fā)電機(jī)分別處于運(yùn)行、降額以及停運(yùn)狀態(tài)的概率
????TR=8760/（lambdaRD+lambdaRF）;TD=8760/（lambdaDR+lambdaDF）;TF=8760/（lambdaFD+lambdaFR）;
????
????WTGnum=25;
????
????R1=rand（WTGnum1000）;????%R1確定風(fēng)力發(fā)電機(jī)所處的狀態(tài)
????R2=rand（WTGnum1000）;????%R2確定該狀態(tài)所持續(xù)的時(shí)間
????D=zeros（WTGnum1000）;????%D記錄每個(gè)狀態(tài)所處的時(shí)間
????time=zeros（WTGnum1001）;????%time記錄每個(gè)狀態(tài)變化的時(shí)間節(jié)點(diǎn)
????alpha1=zeros（WTGnum1000）;????%alpha1記錄每個(gè)狀態(tài)分別是什么
????
????N1

?屬性????????????大小?????日期????時(shí)間???名稱
-----------?---------??----------?-----??----

?????文件???????9037??2019-01-01?22:05??fengdianchang.m

-----------?---------??----------?-----??----

?????????????????9037????????????????????1