基于交替迭代法的通用型交直流潮流计算方法(matlab版+python版)

程序名称##

基于交替迭代法的通用型交直流潮流计算(matlab版+python版)

程序功能(对象)

1. 适用于任意大小的交流电网,支持节点和支路的增删；
2. 适用于接入多个直流环节，直流环节内允许接入风光储等DG；
3. 换流器采用多种控制方式的VS,且认为换流器自身损耗不可忽略；
4. 通用型分析：1、交流电网拓扑允许随意变动；
   2、直流电网拓扑允许随意变动，直流环节个数及位置允许随意变动；
   3、交流电网与直流环节均可接入风光储等多个分布式电源；
   4、VSC的控制方式与接入位置允许随意变动；
5. 通用型限制：去除直流环节需要一个交流系统，否则需要单独编写专用程序，不再通用(大家也可以尝试探讨修改思路)，直流环节个数需要提前知道；
6. 误差分析：与统一迭代法进行对比分析，误差<10^-2；

程序函数说明

matlab版：

1. 主程序：ACDC_Flow
2. 交流系统参数输入子程序：所有ac_case文件
3. 直流系统参数输入子程序：所有dc_case文件
4. 牛拉法基波子程序：NLPowerFlow文件
5. 节点导纳子程序：createYbus文件
6. 交流雅克比矩阵子程序：Jacobi文件(n-1+m*n-1+m)
7. 直流雅克比矩阵子程序：CreateDCJac文件
   图一为改进的IEEE9节点系统算例程序列表
   图二为改进的IEEE33节点系统算例程序列表
   
   

python版：

1. 主程序输入输出：ACDC_Flow_big_test.py
2. 子函数包：ACDC_Flow.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 createDCJac(Y,Udc):

2.8交直流潮流子函数：

def ACDC_FLOW(sheets,DC_baseMVA,DC_bus,DC_branch,DC_VSC,AC_baseMVA,AC_bus,AC_branch,AC_gen,AC_DG,GBUS,YBUS,base_loss):

程序计算步骤及流程图

以matlab版程序作为讲解。

1 参数初始化

参数初始化，读取网络参数，在此标幺化，随机潮流是在确定性基波潮流计算基础之上进行的，参数初始化和基波潮流学习： https://blog.csdn.net/WConstelltion/article/details/123751611.
包括：交流节点负荷，支路阻抗，发电机以及分布式电源参数初始化；直流节点负荷，支路电阻以及分布式电源参数初始；VSC的控制方式以及接入容量和位置初始化，7种控制方式如下图所示。

交替迭代法的计算流程图如下图所示，后续程序依据该逻辑进行讲解

2功率站潮流计算

%目的：计算功率站损耗、功率站直流注入功率Pcdc、直流节点注入功率Pdci
%方式1,2的初值Ps,Qs,Us=1.0,δs=0；方式3的初值Ps,Qs=0,Us,δs=0；
2.1功率站直流侧功率及损耗计算：(Ic为换流器内部流过交流的电流)

    Uc_p = Uf_p-Ic.*Zc;                                         %功率站等效电压源电压
    Sc_p = Uc_p.*Ic; 
Pcloss = a.*abs(Ic).*abs(Ic)+b.*abs(Ic)+c; 

2.2节点注入直流网络功率计算：

    Pdc = -dc_bus(:,3);                                           %直流负荷注入功率
    Pdc(p_vsc) = Pcdc(p_vsc)+Pdc(p_vsc);               %%含VSC注入功率和负荷注入功率，也可增加DG

2.3直流潮流计算：

        Pdc = Udcs.*(Gbus*Udcs);                            %计算Pdci
        dPdc = Pdcs - Pdc;                                      %功率偏差(不是两次迭代之差么？）
        if max(abs(dPdc(p))) <= 1e-9                       %非平衡节点收敛判据
            break
        end
        %直流网络的雅克比矩阵，对电压求导获得Gij*Ui+Gii*Ui;
        J = Gbus.*Udcs + sparse(1:dc_nb,1:dc_nb,Gbus*Udcs,dc_nb,dc_nb); %Sparse(i,j,s,m,n)构造满足S(i(k),j(k)) = s(k)的m*n维稀疏矩阵,省略非零源(算法可靠！)
       % J1 = CreateDCJac(Gbus,Udcs);
        J = J(2:dc_nb,2:dc_nb);                                 %不含平衡节点雅克比
        dU = J\dPdc(p);
        Udcs(p) = Udcs(p) + dU;                              %平衡点以外电压修正
        kdc = kdc+1;

2.4功率站电压越限校验：

         p1 = find(Udcs(p_vsc)<0.95);                       %%判断功率站节点是否越限

3 电压站潮流计算

%目的：内部等效节点参数Uc，δc，Uf，δf以及Ss
3.1通过换流器内部有功和无功平衡方程计算内部电压电压Uc，进而求解交流侧容量Sc

        Ic = (Uc-Uf)/Zc;
        It = (Uf-Us)/Zt;
        If = -Uf*(1j*Bf(v_vsc)); 
        %牛顿法求解电压站等效电路参数
        q1 = Pcdc(v_vsc)+real(Uc.*Ic)+a.*abs(Ic).*abs(Ic)+b.*abs(Ic)+c==0;  %交直侧有功平衡，不是Pdcs而是Pcdc???
        q2 = real(Uf.*Ic)-real(Uf.*It)==0;                  %滤波器节点f有功平衡
        q3 = Qs-imag(Uf.*It-Zt.*It.*It)==0;               %变压器支路无功平衡
        q4 = imag(Uf.*Ic)-imag(Uf.*It)+imag(Uf.*If)==0;%f节点无功平衡，滤波器的无功imag(Uf.*If)，电压站的Bf为零此处直接省略滤波之路
    
        [c_e, c_f,f_e,f_f] = vpasolve(q1,q2,q3,q4,c_e, c_f,f_e,f_f );%%求解非线性方程组
        Uf = f_e+f_f*1i;
        Uc = c_e+c_f*1i;
    Ss =Uc.*(Uc-Uf)/Zc-Zc.*((Uc-Uf)/Zc).*((Uc-Uf)/Zc)+Uf.*If-Zt.*It.*It;  %电压站的初次无功为零

2.2交流系统参数修正:(节点类型，换流器交流侧功率)

    ac_bus(ac_pq,[3]) = ac_bus(ac_pq,[3])+Ps(vsc_pq);
    ac_bus(ac_pq,[4]) = ac_bus(ac_pq,[4])+Qs(vsc_pq);%%修正VSC交流侧PQ节点母线功率  

    ac_bus(ac_pv,[2]) = 2;                                       %VSC的电压站控制方式存在导致原PQ节点变为PV节点 
    ac_bus(ac_pv,[8]) = abs(Us(vsc_pv));                  %变换后pv节点电压修正

2.3交流潮流计算：

    [Pac,Qac,Bus_V,kac]=NLPowerFlow(ac_baseMVA, ac_bus,ac_branch,ac_gen,ac_dg);

2.4定电压Us换流站Qs越限校验：

           p2 = find(Qs(vsc_pv)>0.48);

2.5全局收敛判据：

     if (max(abs(dUac_m)) <= 1e-9)&&(max(abs(dUac_a)) <= 1e-9)&&(max(abs(dUdcs)) <= 1e-9)&&(max(abs(dPdc)) <= 1e-9)&&(max(abs(dQs)) <= 1e-9)
        break
     else

2.6 变量迭代修正：

         VSC(:,8) = Uac_m(ac_node);                         %VSC的交流侧Usi
         VSC(:,9) = Uac_a(ac_node);                          %更新δsi        
         VSC(vsc_pv ,7) = Qs(vsc_pv);                        %3VSC换流站交流侧等效PV节点;

程序误差分析

1. IEEE9节点算例拓扑计算结果潮流分布图如下图所示
   
   统一迭代法和交替迭代法误差对比，误差小于3%
   
2. python版算结果与MATLAB版结果一致，且较matlab版更加通用，可以允许不知道直流环节个数就可计算。

程序“配送”

免费版(未计及DG): https://download.csdn.net/download/lisha1002/12148929?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522164835259216780366543218%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fdownload.%2522%257D&request_id=164835259216780366543218&biz_id=1&utm_medium=distribute.pc_search_result.none-task-download-2_downloadfirst_rank_ecpm_v1~rank_v31_ecpm-5-12148929.pc_v2_rank_dl_default&utm_term=%E5%8D%8A%E4%B8%8D%E5%8F%98%E9%87%8F&spm=1018.2226.3001.4451.5.
完整版matlab版/python版私信V：lemonyoungman