全是自己写的，欢迎交流。

clear;
clc;
%蒙特卡洛模拟次数
times=5000;
%光伏有功服从Beta分布，其中相关的两个参数α与β
a_pv=39.40;
b_pv=26.21;
%光伏发电的相关参数：组件总面积S_pv（m^2）、光电转化效率prey_pv、最大光强rmax（MW/m^2）
S_pv=5000;
prey_pv=0.2;
rmax=0.0015;
%光伏出力场景样本
pv_samp(1,:)=betarnd(a_pv,b_pv,1,times);
Ppv_samp(1,:)=pv_samp(1,:)*rmax*S_pv*prey_pv;
%风速威布尔分布的两个形状参数
k_wt=1.637;
c_wt=5.218;
%风力发电机出力场景样本
Pwt_samp=zeros(1,times);
wt_samp=wblrnd(c_wt,k_wt,1,times);
PN_wt=1.5;
vci=3;
vN=14;
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)=(PN_wt*(wt_samp(i)-vci))/(vN-vci);
    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
%获取33节点的运行数据
mpc=case33bw;
%提取33节点的有功无功注入功率
P=mpc.bus(:,3);
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),Q(k)*ld_ero,1,times);
end
%构建未发生N-1故障时33节点原始的关联矩阵
Gra=zeros(33);
for Guani=1:32
    Gra(mpc.branch(Guani,1),mpc.branch(Guani,2))=1;
    Gra(mpc.branch(Guani,2),mpc.branch(Guani,1))=1;
end
%相应场景的潮流计算
%更新潮流计算的负荷数据
line_flow=zeros(5000,32);
Vbus=zeros(5000,33);
Minsum=zeros(32,5000);
%用于判断配电网联通性的单位矩阵
Danwei=eye(33);
for jk=1:times
    mpc.bus(:,3)=Pld_samp(:,jk);
    mpc.bus(:,4)=Qld_samp(:,jk);
    %在33节点模型中的16节点接入光伏，更新接入光伏的功率
    mpc.bus(16,3)=mpc.bus(16,3)-Ppv_samp(1,jk);
    %在33节点模型中的26节点接入风电，更新接入风电的功率
    mpc.bus(26,3)=mpc.bus(26,3)-Pwt_samp(1,jk);
    mpc.bus(26,4)=mpc.bus(26,4)-Pwt_samp(1,jk)*0.484;
    %接入分布式新能源的配电网潮流计算
    result=runpf(mpc);
    Suml1=0;
    for Sumli=2:33
        Suml1=Suml1+mpc.bus(Sumli,3);
    end
    Minsum(1,jk)=Suml1;
    %得出不同运行场景流过线路的功率（不同行表示不同运行场景）
    line_flow(jk,1:32)=result.branch(1:32,14);
    %得到不同运行场景各节点电压（不同行表示不同运行场景）
    Vbus(jk,:)=result.bus(:,8);
    %33节点配电网发生N-1故障遍历配电网除第一条支路以外的所有支路
    for Gi=2:32
        %断开相应的支路
        mpc.branch(Gi,11)=0;
        %更改相应的邻接矩阵
        Gra(mpc.branch(Gi,1),mpc.branch(Gi,2))=0;
        Gra(mpc.branch(Gi,2),mpc.branch(Gi,1))=0;
        %采用5条不同的联络开关进行转供配电网的潮流计算
        %Jishui是用来统计系统发生N-1故障后不同联络开关系统为联通的次数
        Jishui=0;
        %将不同转供方法的损失量先置零
        Sunzong=[5 5 5 5 5];
        for Li=33:37
            %闭合不同的联络开关
            mpc.branch(Li,11)=1;
            %生成配电网发生N-1+1故障后配电网的关联矩阵
            Gra(mpc.branch(Li,1),mpc.branch(Li,2))=1;
            Gra(mpc.branch(Li,2),mpc.branch(Li,1))=1;
            %判断配电网进行N-1+1过程后系统的联通性
            Keda=Danwei|Gra;
            for j=2:32
                Keda=Keda|Gra^j;
            end
            Zong=sum(Keda(:));
            %如果配电网我满足联通性则系统进行潮流计算
            if Zong==1089
                Jishui=Jishui+1;
                result=runpf(mpc);
                %根据节点电压越限程度计算负荷损失量
                Sunshi=0;
                for Jie=2:33
                    if result.bus(Jie,8)<0.95
                        Ki=(0.95-result.bus(Jie,8))/0.95;
                        if Ki<=0.368
                            Kig=Ki/0.368;
                        else
                            Kig=1;
                        end
                    elseif result.bus(Jie,8)>1.05
                        Ki=(result.bus(Jie,8)-1.05)/1.05;
                        if Ki<=0.143
                            Kig=Ki/0.143;
                        else
                            Kig=1;
                        end
                    else
                        Kig=0;
                    end
                    Sunshi=Sunshi+Kig*result.bus(Jie,3);
                end
                Sunzong(1,Jishui)=Sunshi;
            end
            %将配电网N-1+1运行方案恢复，一并恢复原有邻接矩阵
            mpc.branch(Li,11)=0;
            Gra(mpc.branch(Li,1),mpc.branch(Li,2))=0;
            Gra(mpc.branch(Li,2),mpc.branch(Li,1))=0;
        end
        %选择负荷损失最小作为负荷损失量
        Minsum(Gi,jk)=min(Sunzong,[],'all');
        %将配电网N-1故障恢复，并恢复原有邻接矩阵
        mpc.branch(Gi,11)=1;
        Gra(mpc.branch(Gi,1),mpc.branch(Gi,2))=1;
        Gra(mpc.branch(Gi,2),mpc.branch(Gi,1))=1;
    end
end