《量化投资:以MATLAB为工具》

MATLAB技术论坛

 找回密码
 注册账号
查看: 4275|回复: 5
收起左侧

[其它] 关于蒙特卡罗法计算电力系统可靠性指标

[复制链接]
发表于 2012-4-16 13:04:24 | 显示全部楼层 |阅读模式
文献求助
标题: 求关于蒙特卡罗法计算电力系统可靠性指标程序的详细注释
赏金: 1 贝壳
链接: -
数据库:
关键字:
作者:
刊名:
时间:
备注: 求程序详细注释 做毕业设计需要
function [MVAbase, bus, gen, branch, success, et] =runpf
[baseMVA, bus, gen, branch] = loadcase('caseRTS79');
[i2e, bus, gen, branch] = ext2int(bus, gen, branch);
[probline,probgen]=failprob;
[A,lpr,equ,Pgmax,goalA,busPg]=loadpro;

limB=zeros(1,48);             %limB是1x48的全0矩阵
ranbr=size(branch,1);         %ranbr=矩阵branch的行数
lineB=zeros(ranbr,ranbr);     %lineB是ranbr x ranbr的全0矩阵
for i=1:ranbr                 %i从0到ranbr
    lineB(i,i)=1/branch(i,4); %方阵lineB的对角元素分别是1除以branch第4列的相应行数
end
Pload=bus(:,3);               %Pload是取矩阵bus的第3列的所有元素
Pload(13,:)=[];               %删除Pload的第13行的所有元素
sumload=0;                    %定义sumload=0
for i=1:size(bus,1)           %i从1到矩阵bus的行数
    sumload=sumload+bus(i,3);
end                           %sumload=矩阵bus第3列所有元素之和
sumpg=0;                      %定义sumpg=0
for i=1:length(busPg)         %i从1到矩阵busPg的长度
    sumpg=sumpg+busPg(i,1);
end                           %sumpg=busPg第1列所有元素之和
refPg=591-sumload+sumpg;      
Pmax=branch(:,8);             %Pmax是矩阵branch第8列的所有元素
lolp=0;                       %定义电力不足概率LOLP=0
edns=0;                       %定义缺供期望电力EDNS=0
vari=0;                       %
sumcut=0;                     %定义sumcut=0
sumsqcut=0;                   %定义sumsqcut=0
B=[];
state=[];
for stct=1:50000
    stvari=mc(probline,probgen);
    lengthst=length(stvari);
    numstate=length(state);
    lolp=lolp*(stct-1)/stct;
    edns=edns*(stct-1)/stct;
         ednsarray(1,stct)=edns;
     lolparray(1,stct)=lolp;

    if ~lengthst
          vari=sumsqcut-2*sumcut*edns+stct*edns^2;
       vari=vari/stct^2;
       vindex(1,stct)=sqrt(vari)/edns;
       ednsarray(1,stct)=edns;
       lolparray(1,stct)=lolp;
       continue;
    else
        flag=0;
        for k=1:length(state)
            if lengthst==length(state(1,k).st);
                if stvari==state(1,k).st
                    state(1,k).num=state(1,k).num+1;
                    flag=1;
                    break;
                end
            end
        end
        if ~flag
            state(1,numstate+1).st=stvari;
            state(1,numstate+1).num=1;
        end
    end
    if flag
        if state(1,k).cutload
             sumcut=sumcut+state(1,k).cutload;
            sumsqcut=sumsqcut+state(1,k).cutload^2;
            lolp=lolp+1/stct;
            edns=edns+state(1,k).cutload/stct;
                        vari=sumsqcut-2*sumcut*edns+stct*edns^2;
       vari=vari/stct^2;
                        ednsarray(1,stct)=edns;
            lolparray(1,stct)=lolp;
        end
        vindex(1,stct)=sqrt(vari)/edns;
        continue;
    end
    clear stvari;

    ischange=0;
    sPgmax=Pgmax;
    sbusPg=busPg;
    srefPg=refPg;
    outbr=0;
    outgen=0;
    for lenct=1:length(state(1,length(state)).st)
        if state(1,length(state)).st(1,lenct)<39
            outbr=outbr+1;
            branch(state(1,length(state)).st(1,lenct),11)=0;
            memobr(1,outbr).loc=state(1,length(state)).st(1,lenct);
            memobr(1,outbr).b=lineB(state(1,length(state)).st(1,lenct),4);
            lineB(state(1,length(state)).st(1,lenct),4)=0;
            ischange=1;
            clear B;
           
        else
            gavri=state(1,length(state)).st(1,lenct)-38;
            gen(gavri,8)=0;
            srefPg=srefPg-gen(gavri,2);
            outgen=outgen+1;
            memogen(1,outgen)=gavri;
            if gen(gavri,1)<13
                sPgmax(1,gen(gavri,1))=sPgmax(1,gen(gavri,1))-gen(gavri,9);
                sbusPg(gen(gavri,1),1)=sbusPg(gen(gavri,1),1)-gen(gavri,2);
            end
            if gen(gavri,1)==13
                srefPg=-1;
                sPgmax(1,24)=Pgmax(1,24)-gen(gavri,9);
            end
            if gen(gavri,1)>13
                sPgmax(1,gen(gavri,1)-1)=sPgmax(1,gen(gavri,1)-1)-gen(gavri,9);
                sbusPg(gen(gavri,1)-1,1)=sbusPg(gen(gavri,1)-1,1)-gen(gavri,2);
            end
        end
    end
%       if (stct==1)|ischange
        B = makeBdc(baseMVA, bus, branch);
        subB=full(B);
        subB(13,:)=[];
        subB(:,13)=[];
        swp=lineB*A*inv(subB);
        swp1=swp*Pload;
        maxArray=Pmax+swp1;
        minArray=swp1-Pmax;
        maxArray=[maxArray;-minArray];
        lprA=swp*lpr;
        lprA=[lprA;-lprA];
        clear minArray
        clear B
        clear subB
%       end
   
    state(1,length(state)).cutload=0.0;
    if srefPg>0
        brflow=swp*(sbusPg-Pload);
        cutload=0;
        for ctbranch=1:38
            if abs(brflow(ctbranch,1))>branch(ctbranch,8)
                limA=[Pload',bus(13,3),sPgmax];
                [m,cutload]=linprog(goalA,lprA,maxArray,equ,sumload,limB,limA);
                if cutload>1
                    state(1,length(state)).cutload=cutload;
                end
                break;
            end
        end
    else
        limA=[Pload',bus(13,3),sPgmax];
        [m,cutload]=linprog(goalA,lprA,maxArray,equ,sumload,limB,limA);
        if cutload>1
             state(1,length(state)).cutload=cutload;
        end
    end
    if state(1,length(state)).cutload
                    sumcut=sumcut+state(1,length(state)).cutload;
            sumsqcut=sumsqcut+state(1,length(state)).cutload^2;
        lolp=lolp+1/stct;
        edns=edns+state(1,length(state)).cutload/stct;
         vari=sumsqcut-2*sumcut*edns+stct*edns^2;
        vari=vari/stct^2;
        ednsarray(1,stct)=edns;
        lolparray(1,stct)=lolp;
    end
    vindex(1,stct)=sqrt(vari)/edns;
    success = 1;
    for i=1:outbr
        branch(memobr(1,i).loc,11)=1;
        lineB(memobr(1,i).loc,4)=memobr(1,i).b;
    end
    for i=1:outgen
        gen(memogen(1,i),8)=1;
    end
    clear memobr;
    clear memogen;
%     if (stct>10)&(vindex(1,stct)<0.017)
%         break
%     end
end
layer=zeros(1,15);
for i=1:length(state)
    layer(1,length(state(1,i).st))=layer(1,length(state(1,i).st))+state(1,i).num*state(1,i).cutload/stct;
end

lolp
edns
dlmwrite('E:\study\edns1.txt', ednsarray);
dlmwrite('E:\study\lolp1.txt', lolparray);
dlmwrite('E:\study\var1.txt', vindex);
dlmwrite('E:\study\layer1.txt', layer);
plot(vindex);
hold on
plot(layer)
return;
rudeMC.rar (18.16 KB, 下载次数: 100)
发表于 2012-8-24 11:05:01 | 显示全部楼层
本帖最后由 liyandump 于 2012-8-24 11:06 编辑

好东西啊!很多人在找这些源代码都找不到,我老师前几天刚找到发给我,就在这里看到你贴出来。我也正在研究中
发表于 2012-12-21 15:22:32 | 显示全部楼层
好东西谢谢分享
发表于 2013-7-30 10:29:00 | 显示全部楼层
谢谢楼主,找了很久,终于找到了有用的资料。
回复 支持 反对

使用道具 举报

发表于 2015-3-25 11:26:17 | 显示全部楼层
谢谢楼主,好东西下来学习研究
回复 支持 反对

使用道具 举报

发表于 2017-5-16 19:26:31 | 显示全部楼层
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 注册账号

本版积分规则

QQ|网站地图|MATLAB技术论坛|Simulink仿真论坛 ( 蜀ICP备19014457号 

GMT+8, 2020-10-1 08:14 , Processed in 0.091047 second(s), 17 queries , Gzip On, MemCached On.

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表