丁采样算法

潮流计算中的拉丁采样算法的源程序，可运用其进行概率潮流计算

%C1z，C2z，C3z的求解
C1=[1.0 0.5 0.8 0.7 0.6;
    0.5 1.0 0.7 0.6 0.8;
    0.8 0.7 1.0 0.5 0.6;
    0.7 0.6 0.5 1.0 0.9;
    0.6 0.8 0.6 0.9 1.0];
C2=[1.0 0.5 0.8 0.7 0.6 0.4;
    0.5 1.0 0.7 0.6 0.8 0.5;
    0.8 0.7 1.0 0.5 0.6 0.5;
    0.7 0.6 0.5 1.0 0.9 0.6;
    0.6 0.8 0.6 0.9 1.0 0.6;
    0.4 0.5 0.5 0.6 0.6 1.0];
C3=[1.0 0.5;
    0.5 1.0];
k=2.15;
namada=9.0;
Ew=namada*gamma（1+1/k）; %威布尔分布的期望
sigma=sqrt（namada^2*gamma（1+2/k）-Ew^2）;
T=ones（2）;
for i=1:2
    for j=1:2
        if i~=j
         T（ij）=1.063-0.004*C3（ij）-0.2*（sigma/Ew）-0.001*C3（ij）+0.337*（sigma^2/Ew^2）+0.007*C3（ij）*（sigma/Ew）-0.007*（sigma^2/Ew^2）;
        else
            T（ij）=1;
        end
    end
end
C3z=ones（2）;
for i=1:2
    for j=1:2
        C3z（ij）=T（ij）*C3（ij）;
    end
end
C1z=C1;C2z=C2;

%chol分解求B1B2B3;
BB1=（chol（C1z））‘;
BB2=（chol（C2z））‘;
BB3=（chol（C3z））‘;

%利用LHS方法对n个相互独立的标准正态随机分布随机变量进行采样，采样方法见于晗博士论文
%2-6节点负荷采样
N=input（‘请输入采样规模:‘）;
X1=zeros（5N）;   
for i=1:5
    x1=rand（1N）;
    for n=1:N
        X1（in）=（x1（1n）+n-1）/N;
    end
end
W1=norminv（X101）;   %LHS采样
%排序
L1=zeros（5N）;  %随机数据
for i=1:5
    L1（i:）=randperm（N）;
end
Pl1=corrcoef（L1‘）;   %相关系数矩阵

D1=（chol（Pl1））‘;    %相关系数矩阵分解出的下三角矩阵
D1=inv（D1）*L1;      %顺序矩阵参考数据
[ignoreLd1]=sort（D12）;
Ld1;       %顺序矩阵

A=sort（W12）;
W11=zeros（5N）;
for i=1:5
    for j=1:N
        W11（iLd1（ij））=A（ij）;
    end
end
W11;

%对9-14节点采样
X2=zeros（6N）;
for i=1:6
    x2=rand（1N）;
    for n=1:N
        X2（in）=（x2（1n）+n-1）/N;
    end
end
W2=norminv（X201）;
%排序
L2=zeros（6N）;  %随机数据
for i=1:6
    L2（i:）=randperm（N）;
end
Pl2=corrcoef（L2‘）;   %相关系数矩阵

D2=（chol（Pl2））‘;    %相关系数矩阵分解出的下三角矩阵
D2=inv（D2）*L2;      %顺序矩阵参考数据
[ignoreLd2]=sort（D22）;
Ld2;       %顺序矩阵

A=sort（W22）;
W22=zeros（6N）;
for i=1:6
    for j=1:N
        W22（iLd2（ij））=A（ij）;
    end
end
W22;

%对风机负荷采样
X3=zeros（2N）;
for i=1:2
    x3=rand（1N）;
    for n=1:N
        X3（in）=（x3（1n）+n-1）/N;
    end
end
W3=norminv（X301）;
%排序
L3=zeros（2N）;  %随机数据
for i=1:2
    L3（i:）=randperm（N）;
end
Pl3=corrcoef（L3‘）;   %相关系数矩阵

D3=（chol（Pl3））‘;    %相关系数矩阵分解出的下三角矩阵
D3=inv（D3）*L3;      %顺序矩阵参考数据
[ignoreLd3]=sort（D32）;
Ld3;       %顺序矩阵

A=sort（W32）;
W33=zeros（2N）;
for i=1:2
    for j=1:N
        W33（iLd3（ij））=A（ij）;
    end
end
W33;


%求Z1Z2Z3
Z1=BB1*W11;
Z2=BB2*W22;
Z3=BB3*W33;


%对输入随机变量X进行LHS采样，并按顺序矩阵Ls（上面求得的Z1，Z2，Z3）进行排序，得到最终的样本矩阵S11S22S33
%求2-6节点的最终样本矩阵（不带风机）
Orig1=[0.217 0.06 0.863
      0.942 0.3  0.98
      0.478 0.11 0.31
      0.076 0.03 0.979
      0.112 0.03 0.831];
X1=zeros（1N）;
S1=zeros（5N）;
for i=1:5
  X1=rand（1N）;
    for n=1:N 
        X1（1n）=（X1（1n）+n-1）/N;
    end
    S1（i:）=norminv（X1Orig1（i1）Orig1（i2））;
end

[ignoreLs1]=sort（Z12）;
Ls1;       %顺序矩阵

A1=sort（S12）;
S11=zeros（5N）;
for i=1:5;
    for j=1:N;
        S11（iLs1（ij））=A1（ij）;
    end
end

%求9-14节点的最终样本矩阵

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        4509  2013-05-01 10:25  BandC.m