关于蒙特卡罗法计算电力系统可靠性指标程序的详细注释
function =runpf
= loadcase('caseRTS79');
= ext2int(bus, gen, branch);
=failprob;
=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=;
lprA=swp*lpr;
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=;
=linprog(goalA,lprA,maxArray,equ,sumload,limB,limA);
if cutload>1
state(1,length(state)).cutload=cutload;
end
break;
end
end
else
limA=;
=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;
好复杂的样子
我也用过蒙特卡洛，可以交流下
蒙特卡罗算法在MATLAB中怎么实现呀，还有随机数怎么生成？跪求帮助！
蒙特卡罗算法在MATLAB中怎么实现呀，还有随机数怎么生成？跪求帮助！
好好学习一下
刚开始学习
慢慢来，希望能提高自己的能力
页:
[1]