序贯蒙特卡洛可靠性评估.rar

对含储能和风电的电力系统进行了可靠性评估，利用序贯蒙特卡洛法，把风储系统接入IEEE-RBTS系统来仿真，探讨了风电场、储能系统、储能容量和储能最大充放电功率对系统可靠性的具体影响，可以运行。

clear;
clc;
tic

year=2
%模拟的年限

for l=1:year;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%  第一步：数据导入与预处理  %%%%%%%%%%%%%%%%%%%
    
    SW0=load（‘windspeed.txt‘）;    %载入原始风速数据
    SW0=SW0‘;                     %将风速数据从横列转置成竖列
    SW0=SW0/10*36;                %原始数据的风速单位为m/s，这里转化为km/h
    N=size（SW02）;                %把SW0的列数赋值给N
    mu=mean（SW0）;                 %求SW0的平均值
    sigma=var（SW0）;
    sigma=sigma^0.5;              %求SW0的标准差
    y=（SW0-mu）./sigma;            %数据预处理
    
    N1=8736;                      %N1为模拟的小时数
    WTGnum=25;                    %WTGnum为风力发电机数量
    
    %%%%%%%%%%%%%%  第二步：根据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=8736/（lambdaRD+lambdaRF）;TD=8736/（lambdaDR+lambdaDF）;TF=8736/（lambdaFD+lambdaFR）;
    
    
    R1=rand（WTGnum1000）;         %R1确定风力发电机所处的状态
    R2=rand（WTGnum1000）;         %R2确定该状态所持续的时间
    D=zeros（WTGnum1000）;         %D记录每个状态所处的时间
    time=zeros（WTGnu

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

     文件        251  2017-04-09 20:59  windspeed.txt

     文件      69317  2013-01-19 19:43  load.txt

     文件      11176  2017-04-17 17:03  wind.m

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

                80744                    3