蒙特卡洛算法计算电力系统可靠度matlab程序

蒙特卡洛模拟法是一种模拟法，由电子计算机来模拟一个过程的实现，重复 一定次数，然后计算系统的风险指标。它是用随机数来模拟数学与物理问题以求其近似解的一种通用方法。它可以解决带有随机性的问题和确定性问题；处理具有指数型分布和非指数型分布的问题；以解决易于建立数学模型的问题和不能建立数学模型的问题或建立了数学模型而难以求得数值解的问题。

clear;clc;
%% 系统数据
%  编号  类型  幅值  相角  发电有功  发电无功  负荷有功  负荷无功
NODE=[ 1  3  1.04  0  0        0  0    0
       2  2  1.025 0  1.63     0  0    0
       3  2  1.025 0  0.85     0  0    0
       4  0  1     0  0        0  0    0
       5  0  1     0  0        0  1.25 0.50
       6  0  1     0  0        0  0.9  0.30
       7  0  1     0  0        0  0    0
       8  0  1     0  0        0  1    0.35
       9  0  1     0  0        0  0    0
    ];
%  编号1 编号2  类型  R  X  B/2  K
BRANCH=[ 2  7  1  0      0.0625  0      1
         3  9  1  0      0.0586  0      1
         1  4  1  0      0.0576  0      1
         4  5  0  0.010  0.085   0.088  1 
         4  6  0  0.017  0.092   0.079  1
         5  7  0  0.032  0.161   0.153  1
         6  9  0  0.039  0.170   0.179  1
         8  9  0  0.0119 0.1008  0.1045 1
         7  8  0  0.0085 0.072   0.0745 1
    ];
%发电机参数：  发电机编号    H       xd‘       xd         xq       D 
GEN=[          1       23.64   0.0608    0.146      0.0969     0
               2        6.4    0.1198    0.8958     0.8645     0
               3        3.01   0.1813    1.3125     1.2578     0         ];
%

%% 节点导纳矩阵形成
a=size（NODE）; num_node=a（1）;
b=size（BRANCH）; num_branch=b（1）;
G=zeros（num_node）; B=zeros（num_node）;
for ii=1:num_branch
    ST  = BRANCH（ii1）;   
    ED  = BRANCH（ii2）;   
    R   = BRANCH（ii4）;
    X   = BRANCH（ii5）;
    BB  = BRANCH（ii6）;
    K   = BRANCH（ii7）;   
    RX2 = R^2+X^2;
    G（STED） = -K^2*R/RX2;
    G（EDST） = -K^2*R/RX2;
    B（STED） =  K^2*X/RX2;
    B（EDST） =  K^2*X/RX2;
    G（STST） = G（STST） + K^2*R/RX2;
    B（STST） = B（STST） - K^2*X/RX2;
    B（STST） = B（STST） + BB;
    G（EDED） = G（EDED） + K^2*R/RX2;
    B（EDED） = B（EDED） - K^2*X/RX2;
    B（EDED） = B（EDED） + BB;
end
clear ii R X BB RX2 ST ED;

%% 稳态计算
time=0;
while（time==0||flag>1e-5）
u=NODE（:3）; th=NODE（:4）; pg=NODE（:5）-NODE（:7）;qg=NODE（:6）-NODE（:8）;
kp=1;kq=1;sum_pi=0;sum_qi=0;sum_qi1=0;dp=zeros（81）;dq=zeros（61）;pi=zeros（91）;qi=zeros（91）;
%功率平衡方程
for ii=2:9
    for jj=1:9
        sum_pi=sum_pi+u（jj）*（G（iijj）*cos（th（ii）-th（jj））+B（iijj）*sin（th（ii）-th（jj）））;
        if ii>3
            sum_qi=sum_qi+u（jj）*（G（iijj）*sin（th（ii）-th（jj））-B（iijj）*cos（th（ii）-th（jj）））;
        end
    end
    dp（kp1）=pg（ii）-u（ii）*sum_pi;
    kp=kp+1;
    if ii>3
        dq（kq1）=qg（ii）-u（ii）*sum_qi;
        kq=kq+1;
    end
    sum_pi=0;sum_qi=0;
end
clear ii jj kp kq sum_pi sum_qi;

%雅克比矩阵
sum_pi1=0;sum_qi1=0;
for ii=1:9
    for jj=1:9
        sum_pi1=sum_pi1+u（jj）*（G（iijj）*cos（th（ii）-th（jj））+B（iijj）*sin（th（ii）-th（jj）））;
        sum_qi1=sum_qi1+u（jj）*（G（iijj）*sin（th（ii）-th（jj））-B（iijj）*cos（th（ii）-th（jj）））;
    end
    pi（ii1）=u（ii）*sum_pi1;
    qi（ii1）=u（ii）*sum_qi1;
    sum_pi1=0;sum_qi1=0;
end
clear ii jj sum_pi1 sum_qi1;

H=zeros（8）;N=zeros（86）;K=zeros（68）;L=zeros（6）;
for ii=2:9
    for jj=2:9
        if ii==jj  H（ii-1ii-1）=u（ii）^2*B（iiii）+qi（i

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        6800  2018-11-06 15:49  WSCC.m