【直流潮流】基于直流潮流的电力系统停电分布及自组织临界性分析

1.软件版本

matlab2013b

2.系统概述

整个算法的流程图，OPA模型的流程图如下：

 详细的讲，我们的这个算法的安如下的步骤进行：

步骤1：k=k+1，通过慢动态过程中的几个公式，对Pmax，Fmax进行更新；

步骤2：考虑随机因素进行线路的断开，以一个随机概率来随机断开一条支路；

步骤3：根据慢动态计算得到的参数开始进行慢动态仿真；

步骤4：在慢动态仿真中，使用粒子群算法来计算论文公式5的最小值；

步骤5：进一步计算注入功率和潮流；

步骤6：根据计算得到的潮流和Fmax，根据F/Fmaz > alpha来确定线路是否接近容量极限；

步骤7：如果满足F/Fmaz > alpha，则以一个概率beta进行断开。

步骤8：如果发生断开时间，则返回步骤1，如果没有发生断开事件，则结束当前的动态过程。

3.部分源码

clc;
clear;
close all;
warning off;
addpath 'FUNC\'

%%
MOD      = 4;%选择3种不同的参数
if MOD < 4
%参数初始化
k        = 1;
kmax     = 100000;%点数少的话，就没法得到平滑的仿真图，    
K2       = 2;
sel      = 1;      %1：选择IEEE30模型，2：选择IEEE118模型；
alpha    = 0.7;    %Fjk/Fjk,max
if MOD == 1
beta     = 0.3;   %开断概率
end
if MOD == 2
beta     = 0.7;   %开断概率
end
if MOD == 3
beta     = 1;   %开断概率
end

lemda    = 1.0005; %Pi,(k+1),max/pi,k,max
mu       = 1.003;  %Fi,(k+1),max/Fi,k,max
gamas    = [1,1];  %Pik/Pik,max
gama     = gamas(1) + (gamas(2) - gamas(1))*rand(1,1);


W        = 10;%公式五中的w值
iter_max = 2000;%粒子群优化算法得到最小的值（公式5）
MTKL     = 1;%利用蒙特卡洛的设计思想，对一个结果进行多次仿真，然偶进行求平均
%初始化最大值的放大比例
KKK1     = 11.5;
KKK2     = 1;

%节点数的初始化
if sel == 1%IEEE30初始化参数，直流潮流计算采用工具箱进行计算
   mpc                            = case30;
   BUS                            = mpc.bus;
   GEN                            = mpc.gen;
   LINE                           = mpc.branch;
  [r,c]                           = size(mpc.bus);%r为节点数量
   Bus_Num                        = r;
  [r,c]                           = size(mpc.branch);%r为节点数量
   F_Num                          = r;    
end

if sel == 2%IEEE118初始化参数，直流潮流计算采用工具箱进行计算
   mpc                            = case118;
   BUS                            = mpc.bus;
   GEN                            = mpc.gen;
   LINE                           = mpc.branch;   
  [r,c]                           = size(mpc.bus);%r为节点数量  
   Bus_Num                        = r;
  [r,c]                           = size(mpc.branch);%r为节点数量
   F_Num                          = r;      
end


%变量初始化
p      = zeros(Bus_Num,kmax);  %实际负荷需求
F      = zeros(F_Num,kmax);    %潮流
Pmax   = zeros(Bus_Num,kmax);  %Pmax的最大值
Fmax   = zeros(F_Num,kmax);    %Fmax的最大值
CUT    = zeros(F_Num,kmax);    %判断某一天某个支路是否被断开
c      = zeros(Bus_Num,kmax);  %cij

%定义两类集合：发电机节点G和负荷节点L
%定义两类集合：发电机节点G和负荷节点L
G      = zeros(Bus_Num,1);     %发电机的集合
L      = zeros(Bus_Num,1);     %负荷点的集合
PG     = zeros(Bus_Num,kmax);  %
PL     = zeros(Bus_Num,kmax);  %
P      = zeros(Bus_Num,kmax);  %
[G,L]  = func_findGL(Bus_Num,mpc.gen,mpc.bus);


%%
%直流潮流计算
%通过初始的IEEEXX电网的类型，计算出其对应的节点负荷和线路潮流,即计算k = 1的情况
%即计算K = 1的时候，系统初始的F和P
%根据公式5来计算
%根据IEEEXX标准数据得到节点的初始负荷
p(:,1) = mpc.bus(:,3);

%计算导纳矩阵A
Ak = func_Admittance_matrix(mpc.bus,mpc.branch);

%初始的注入功率Pk(发电机出力和负荷)
PG(G,1) =  mpc.gen(:,2);   %发电机节点
c(G,1)  =  1;
PL(L,1) =  mpc.bus(L,3); %负荷节点

%然后根据Fk = Ak * Pk计算潮流
P(:,1) = PG(:,1) + PL(:,1);%这里做相加，其实就是数组的合并，即出力节点和负荷节点的合并变为Pik
F(:,1) = Ak * P(:,1);

%初始化网络结构
BUS    = mpc.bus;
GEN    = mpc.gen;
LINE   = mpc.branch;   

%初始化最大值的参数
%初试情况下，我们将最大值定义所有节点的最大值
%注意，这里最大值的初始约束条件，你也可以自己设定
Pmax(:,1)    = KKK1*P(:,1); %节点的发电能力和负荷
for i = 1:F_Num
    Fmax(i,1)= KKK2*F(i,1);
end

%保存停电次数变量
CUT_NUM = zeros(kmax,1);
%保存停电规模变量
CUT_ARE = zeros(kmax,1);
S1      = 0;

%OPA模型的循环
for k = 2:kmax%慢循环过程
    LINE   = mpc.branch;
    flag   = 0;
    %步骤一
    %步骤一
    %由慢动态确定最后的最大值
    %对于第k天，确定发电机的最大出力和负荷需求后Fi,k,max  Pi,k,max      
    %更新负荷，发电机和部分线路容量
    %更新P
    Pmax(:,k) = lemda*Pmax(:,k-1);
 
    %更新F
    for i = 1:F_Num
        if CUT(i,k-1) == 1%判断有没有开断的线路
           Fmax(i,k) = mu*Fmax(i,k-1);
        else
           Fmax(i,k) =    Fmax(i,k-1); 
        end
    end

    %考虑随机因素进行线路的断开
    %这里以一个0~1的随机数和0.5做比较来判断是否收到随机因素的干扰
    %断开的一个随机位置的线路
    PP = rand(1,1);
    
    if PP >= 0.5%断开
       NUMS = floor(F_Num*rand(1,1));
       if NUMS == 0
          NUMS = 1;
       end
       LINE(NUMS,3:end) = 0;  
%        CUT(NUMS,k)      = 1;
    else%不断开
 
    end
    
    while(flag == 0) %慢动态循环中下进行快动态循环
        
        %负荷节点的浮动
        p(:,k)   = gama*Pmax(:,k-1);%实际的负荷波动 
        
        disp('迭代次数：');
        k
        
        %步骤二
        %步骤二
        %根据公式5来计算出P和F
        %这里利用粒子群的思想，进行函数最小值的求解
        %这里编写了func_pso_calculate_min粒子群优化函数来计算公式5的最小值
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        [V_score2,PP] = func_pso_calculate_min(W,Bus_Num,iter_max,Pmax(:,k),c(:,1),p(:,k)); 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Ak     = func_Admittance_matrix(BUS,LINE);
        F(:,k) = Ak*PP;

        %步骤三
        %步骤三
        %随机开断线路
        %产生一个随机的数作为概率
        TMP = 0;
        for i = 1:F_Num           
            %如果满足如下的条件，则以一个固定的概率进行断开
            if abs(F(i,k)/Fmax(i,k)) >= alpha
               Pp = rand(1,1);  
               if Pp >= 1 - beta%开断概率
                  CUT(i,k) = 1;
                  flag     = 1;
               else
                  CUT(i,k) = 0;
               end
            end 
        end
        
        if flag == 0;
           break;
        end
        
        %步骤四
        %步骤四
        %根据CUT断开矩阵，重新更新网络结构
        %当检测到某个CUT为1的时候，说明要对该线路进行断开
        IND = [];
        cnt = 0;
        for i = 1:F_Num
            if CUT(i,k) == 1;%断开相关的线路
               cnt = cnt + 1;
               LINE(i,3:end) = 0; 
               IND(cnt) = i;
            end
        end
    end
    %这里只是为了对比，所以实际的规模计算见图2,3,4的计算方法，这里简化
    CUT_NUM(k) = sum(CUT(:,k));
    I(k)       = (sum(p(L,k))/(abs(sum(Fmax(G,k)))));
end

%%
%画出最后的仿真结果图，仿真出文献“基于直流潮流的电力系统停电分布及自组织临界性分析”的所有图
%每个仿真文件仿真对应如下的一个图
%(O)Fig.2：改变beta得到不同线路可靠性下的停电分布
%(X)Fig.3：改变mu得到不同线路容量增加方式下的停电分布
%(X)Fig.4：改变gama得到不同负荷波动下的停电分布
%(X)Fig.5：趋于临界状态的过程
%(X)Fig.6：停电概率分布

if MOD == 1
   [n,m] = hist(CUT_NUM,20);
   m = m/max(m);
   save Result\fig2\data1.mat m n
end
if MOD == 2
   [n,m] = hist(CUT_NUM,20);
   m = m/max(m);
   save Result\fig2\data2.mat m n
end
if MOD == 3
   [n,m] = hist(CUT_NUM,20);
   m = m/max(m);
   save Result\fig2\data3.mat m n
end
else
    figure;
    load Result\fig2\data1.mat
    loglog(m,n,'r','LineWidth',2);
    hold on
    load Result\fig2\data2.mat
    loglog(m,n,'b','LineWidth',2);
    hold on   
    load Result\fig2\data3.mat
    loglog(m,n,'m','LineWidth',2);
    hold on        
    xlabel('停电规模');
    ylabel('停电次数');
    axis([0.045,1,1e1,3e3]);
    grid on;
    legend('b=0.3','b=0.7','b=1');
end

4.仿真结论

 

 

 

 

5.参考文献

[1]佚名. 基于直流潮流的电力系统停电分布及自组织临界性分析[J]. 电网技术, 2006, 30(14):6.A02-09