数学建模社区-数学中国

标题: 关于蒙特卡罗法计算电力系统可靠性指标程序的详细注释 [打印本页]

作者: ゞ_轻描丶幸福的 时间: 2015-11-30 10:03
标题: 关于蒙特卡罗法计算电力系统可靠性指标程序的详细注释

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, 下载次数: 8, 售价: 2 点体力)

作者: 吃苹果的梨 时间: 2015-11-30 11:36
好复杂的样子

作者: 851240780 时间: 2015-12-3 16:19
我也用过蒙特卡洛，可以交流下

作者: 2867512731 时间: 2015-12-7 20:40
蒙特卡罗算法在MATLAB中怎么实现呀，还有随机数怎么生成？跪求帮助！

作者: 2867512731 时间: 2015-12-7 20:41
蒙特卡罗算法在MATLAB中怎么实现呀，还有随机数怎么生成？跪求帮助！

作者: FabAcK 时间: 2017-5-9 23:10
好好学习一下

作者: FabAcK 时间: 2017-5-9 23:11
刚开始学习

作者: FabAcK 时间: 2017-5-9 23:12
慢慢来，希望能提高自己的能力

欢迎光临数学建模社区-数学中国 (http://www.madio.net/)