风电场电力系统可靠性评估matlab程序

风电场电力系统可靠性评估的matlab程序，运用蒙特卡洛方法做的！

clc
clear

year=50   %模拟的年限

for l=1:year;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%  第一步：数据导入与预处理  %%%%%%%%%%%%%%%%%%%
    
    %SW0=load（‘windspeed.txt‘）;    %载入原始风速数据
    SW0=xlsread（‘windspeed.xls‘）;
    SW0=SW0‘;
    SW0=SW0/10*3.6;    %原始数据的风速单位为0.1m/s，这里转化为km/h
    N=size（SW02）;
    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模型的阶数  %%%%%%%%%%%%%%%%%%
    
    for n=2:7;
        m=armax（y‘[nn-1]）;
        fai=-m.a;
        theta=m.c;    %把armax函数得到的参数，取出来
        
        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）为预测值
        
        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时，各AIC的值
    [AICn1]=min（arma）;
    n1=n1+1;    %n1为得到的ARMA模型的阶数
    
    %%%%%%%%%%%%%%%%%  第三步：用ARMA模型预测风速并确定风速分布  %%%%%%%%%%%%%%%%%
    
    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）为预测值
    y2;
    SW=y2.*sigma+mu;
    
    %     figure（2）
    %     subplot（121）;
    %     hist（SW0100）;
    %     xlabel（‘风速‘）;ylabel（‘频数‘）;title（‘原始风速的分布‘）;
    %     subplot（122）;
    %     hist（SW100）;
    %     xlabel（‘风速‘）;ylabel（‘频数‘）;title（‘预测风速的分布‘）;
    
    
    %得到8736个小时的预测风速：SW 1*8736
    
    
    %%%%%%%  第四步：风电场的转移过程，确定其一年中三种状态分别存在的时长%%%%%%%
    
    %风电机组的三状态模型
    lambdaRD=5.84;lambdaRF=7.96;lambdaDR=48.3;lambdaFR=58.4;lambdaDF=0;lambdaFD=0;    %风力发电机的3个状态的转移率
    T=[1 1 1;lambdaRD -lambdaDR-lambdaDF lambdaFD;lambdaRF lambdaDF -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（WTGnum1000）;    %R1确定风力发电机所处的状态
    R2=rand（WTGnum1000）;    %R2确定该状态所持续的时间
    D=zeros（WTGnum1000）;    %D记录每个状态所处的时间
    time=zeros（WTGnum1001）;    %time记录每个状态变化的时间节点
    alpha1=zeros（WTGnum1000）;    %alpha1记录每个状态分别是什么
    
    N1

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

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

----------- ---------  ---------- -----  ----

                 9037                    1