基于蒙特卡洛法 分析了风电和光伏的出力情况，之后以IEEE33节点为例进行概率潮流计算（matlab程序实现）

clear; clc; times=5000; %蒙特卡洛仿真次数 baseMVA = 100000;%SB % 光伏有功服从Beta分布 Ppv_samp=zeros(1,times); % Beta分布的两个形状参数 a_pv=0.6869; b_pv=2.1320; % 光伏发电相关参数：组件总面积S_pv、光电转换率prey_pv、最大光强rmax(MW/m2) S_pv=10; prey_pv=0.14; rmax=700; % 光伏有功出力样本 pv_samp(1,:)=betarnd(a_pv,b_pv,1,times); Ppv_samp(1,:)=pv_samp(1,:)*rmax*S_pv*prey_pv/baseMVA; % weibull分布的两个形状参数 k_wt=1.637; c_wt=5.218; wt_samp =wblrnd(c_wt,k_wt,1,times); %产生服从weibull分布的样本 PN_wt=1000; vci=3 ; vN=13; vco=25 ; for i=1:times %得到风电出力样本 if wt_samp(i)<vci Pwt_samp(i)=0; end if wt_samp(i)>vci&&wt_samp(i)<vN Pwt_samp(i)=(wt_samp(i)-vci)/(vN-vci)*PN_wt; if Pwt_samp(i)>PN_wt Pwt_samp(i)=PN_wt; end end if wt_samp(i)>vN&&wt_samp(i)< vco Pwt_samp(i)=PN_wt; end if wt_samp(i)> vco Pwt_samp(i)=0; end end for i=1:times Pwt_samp(i)= Pwt_samp(i)/baseMVA ; end mpc=case33bw; %得到33节点的数据 P=mpc.bus(:,3); %提取33节点的P Q数据 Q=mpc.bus(:,4); ld_ero=0.05; %正态分布的方差 for k=1:33 Pld_samp(k,:)=normrnd(P(k),P(k)*ld_ero,1,times); %%生成服从正态分布的随机数 Qld_samp(k,:)=normrnd(Q(k),P(k)*ld_ero,1,times); end for i=1:times mpc.bus(:,3)= Pld_samp(:,i)'; %更新该次潮流计算的负荷数据 mpc.bus(:,4) =Qld_samp(:,i)'; mpc.bus(14,3)=mpc.bus(14,3)-Ppv_samp(1,i); %更新接入的光伏功率 mpc.bus(14,4)=mpc.bus(14,4)- Ppv_samp(1,i)*0.484; mpc.bus(30,3)=mpc.bus(30,3)- Pwt_samp(1,i); %更新接入的光伏功率 mpc.bus(30,4)=mpc.bus(30,4)-Pwt_samp(1,i)*0.484; result=runpf( mpc); line_flow(i,1:32)=result.branch(1:32,14)'; %得到线路有功功率 Vbus(i,:)=result.bus(:,8)'; %得到电压结果 end