计及DG的基于半不变量和Gram_charlier级数潮流计算(matlab版+python版)

程序名称##

计及DG的基于半不变量和Gram_charlier级数潮流计算(matlab版+python版)

程序功能(对象)

1. 适用于任意大小的纯交流电网,支持节点和支路的增删；
2. 适用于接入多个服从不同概率分布的风电、光伏等分布式电源；
3. 采用Gram_Charlier级数拟合随机变量的概率分布函数；
4. 构建原点矩计算半不变量函数，以及DG原点矩计算函数；
5. 误差分析：与蒙特卡洛计算结果对比，误差<10^-2；

程序函数说明

matlab版：

1. 主程序：new_mian_newload
2. 输入算例子程序：所有case文件
3. 牛拉法基波子程序：NLPowerFlow文件
4. 节点导纳子程序：createYbus文件
5. 雅克比矩阵子程序：Jacobi文件(n-1+m*n-1+m)
6. 节点功率与电压的灵敏度子程序：dSbus_dV文件
7. 雅克比矩阵子程序：formJ文件(2n*2n)
8. 主程序内子函数：
   8.1 pv_aerfa();光伏原点矩函数
   8.2 wt_aerfa;风电原点矩函数
   8.3 aerfa_gama();原点矩计算半不变量函数
   8.4 pdf();Gram_Charlier级数展开函数
   

python版：

1. 主程序输入输出：BBBLflow_test.py
2. 子函数包：BBBLflow.py
   2.1读表函数：

def read_sheet(data_matrix):

2.2节点导纳：

def createYbus(bus,branch):

2.3切片函数：

def slice_matrix(Matrix,Queue_row,Queue_col):

2.4节点注入功率函数：

def calculate_PiiQii(busNum,Bus_V,Ybus):

2.5雅克比矩阵函数：

def jacobi(V,Y,ac_pq,ac_pv,Pi,Qi):

2.6牛拉法函数：

def NL_flow(ac_baseMVA,ac_bus,ac_branch0,ac_gen,ac_bn,ac_pv,ac_pq,m):

2.7光伏原点矩计算函数：

def pv_aerfa(a1,b1,Ppv_max1):

2.8风电原点矩计算函数：

def wt_aerfa(K,C,K1,K2):

2.9原点矩计算半不变量函数：

def aerfa_gama(aerfa):

2.10Gram_Charlier级数展开：

def c2pdf(x,c,sigma):

2.11灵敏度函数：

def dSbus_dV(Ybus,V):

2.12随机潮流计算主函数：

def main(MC_bus,CM_bus,MC_branch,CM_branch,MC_gen,CM_gen,MC_dg0,CM_dg0):

所有子函数matlab版与python版命名保持一致。

程序计算步骤及流程图

以matlab版程序作为讲解。

1 参数初始化

参数初始化，读取网络参数，在此标幺化，随机潮流是在确定性基波潮流计算基础之上进行的，参数初始化和基波潮流学习(极坐标下牛顿—拉夫逊潮流计算(matlab版+python版))： https://blog.csdn.net/WConstelltion/article/details/123773744.

2 蒙特卡洛计算

2.1负荷抽样：

    p_load(i,:)=normrnd(mu_load_p(i),sigma_load_p(i),1,m);%正态分布
    q_load(i,:)=normrnd(mu_load_q(i),sigma_load_q(i),1,m);

2.2DG抽样：(风速、光强、储能功率)

     v_FD = wblrnd(FD_c,FD_k,m,1); %风速抽样
    cd_GF= betarnd(GF_a,GF_b,1,m);   %光强抽样

2.3DG功率计算：

      for i=1:m                                                            
          if v_FD(i)<vci||v_FD(i)>vco
              p_FD(i)=0;
          elseif v_FD(i)>vci&&v_FD(i)<vr
              p_FD(i)=k1*v_FD(i)+k2;
          elseif v_FD(i)>vr&&v_FD(i)<vco
              p_FD(i)=PnW;
          end
     end
                
     p_GF = cd_GF*A*yita/1000;  

2.4修正节点负荷：

     MC_bus(node_DG(k),3)= MC_bus(node_DG(k),3)-p_DG(k,i);
     MC_bus(node_DG(k),4)= MC_bus(node_DG(k),4)-q_DG(k,i);

2.5确定性潮流计算：

    [~,  ~, ~,~,mc_bus,mc_branch, Ymc]=NLPowerFlow(MC_baseMVA, MC_bus ,MC_branch, MC_gen);

3 半不变量计算

本文算法流程引用前辈的计算思路并基于改进：https://download.csdn.net/download/destiny9613/9793309?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522164886605016782246437671%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fall.%2522%257D&request_id=164886605016782246437671&biz_id=1&utm_medium=distribute.pc_search_result.none-task-download-2_allfirst_rank_ecpm_v1~rank_v31_ecpm_bkp-5-9793309.142^v5pc_search_quality_down,157^v4control&utm_term=%E5%8D%8A%E4%B8%8D%E5%8F%98%E9%87%8F&spm=1018.2226.3001.4187.5
3.1负荷半不变量计算：(一阶,二阶)

   gama_p(:,2)=sigma_P_L.^2;   
   gama_q(:,2)=sigma_Q_L.^2;

3.2风电半不变量计算

      aerfa_PW =  wt_aerfa(FD_k, FD_c, k1, k2);
      gama_pw_p = aerfa_gama( aerfa_PW );

3.3光伏半不变量计算

      aerfa_PV = pv_aerfa(GF_a ,GF_b ,Ppv_max); 
      gama_pv_p = aerfa_gama(aerfa_pv_p); 

3.4计及分布式电源接入后的半不变量

      gama_p(node_DG(k),2)=abs(gama_p(node_DG(k),2)-gama_dg_p(k,2));
      CM_bus(node_DG(k),3)=CM_bus(node_DG(k),3)-gama_dg_p(k,1)*MC_baseMVA;

3.5确定性潮流计算

[H0, N0, M0, L0, cm_bus, cm_branch, Ybus]=NLPowerFlow(CM_baseMVA, CM_bus,CM_branch,CM_gen); %基波潮流

3.6节点电压半不变量

s1=Jacco\eye(size(Jacco));
gama_v(:,i)=(s1.^i)*gama_pq(:,i);

3.7支路潮流灵敏度矩阵

G=[G11 G12;G21 G22];%查阅参看文献
t1=G*s1;
gama_xianlu(:,i)=(t1.^i)*gama_pq(:,i);

3.8Gram-charlie展开系数
以节点电压为例分析

g_vm(:,i)=gama_vm(:,i)./(gama_vm(:,2).^(i/2));
c_vm=zeros(size(g_vm));
c_vm(:,3)=-g_vm(:,3);
c_vm(:,4)=g_vm(:,4);
c_vm(:,5)=-g_vm(:,5);
c_vm(:,6)=g_vm(:,6)+10*g_vm(:,3).^2;

3.9概率密度求解
以节点电压为例分析

for i=1:npq
    x=0.95:0.0001:1.04;
    x1=(x-mu_vm(i))/sigma_vm(i);
    pdf_vm(i,:)=c2pdf(x1,c_vm(i,:),sigma_vm(i));
end

8. 计算流程图如下图所示
   

计及DG的半不变量和Gram_chaarlier级数随机潮流计算流程图

程序误差分析

1. 两种方法接入DG前节点电压计算结果对比图，曲线几乎吻合，误差小于2%，这是方法本身的特性造成的
   
   两种方法的部分计算数值对比
   

2. matlab版18节点接入光伏的计算结果：(Bus_V为标幺值)
   

3. python版接入DG之后计算结果：(与MATLAB版结果一致)
   

程序链接

前辈程序(仅考虑负荷波动): https://download.csdn.net/download/destiny9613/9793309?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522164886605016782246437671%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fall.%2522%257D&request_id=164886605016782246437671&biz_id=1&utm_medium=distribute.pc_search_result.none-task-download-2_allfirst_rank_ecpm_v1~rank_v31_ecpm_bkp-5-9793309.142^v5pc_search_quality_down,157^v4control&utm_term=%E5%8D%8A%E4%B8%8D%E5%8F%98%E9%87%8F&spm=1018.2226.3001.4187.5
如需博主计及DG的完整版matlab版/python版(恒流源或诺顿等效)私信，也可以V博主一下：lemonyoungman，欢迎批评指正。