牛顿法——电力系统潮流计算

最新推荐文章于 2024-06-20 08:38:01 发布

小白砰

最新推荐文章于 2024-06-20 08:38:01 发布

阅读量624

点赞数 1

文章标签：算法

本文链接：https://blog.csdn.net/m0_54086968/article/details/131530984

版权

潮流计算

计算算例

算法流程图

如下：

计算过程及代码

1、输入电网及注入节点数据

2、计算导纳阵

%% 输入节点之间阻抗和导纳
Y0=zeros(5,15);
Y0=[0,0,0,complex(0,0.25),1/complex(0.04,0.25),complex(0,0.25),0,1/complex(0.1,0.35),0,0,0,0,0,0,0;
    complex(0,0.25),1/complex(0.04,0.25),complex(0,0.25),0,0,0,complex(0,0.25),1/complex(0.08,0.3),complex(0,0.25),(1-1.05)/complex(0,1.05*1.05*0.015),1/complex(0,1.05*0.015),(1.05-1)/complex(0,1.05*0.015),0,0,0;
    0,1/complex(0.1,0.35),0,complex(0,0.25),1/complex(0.08,0.3),complex(0,0.25),0,0,0,0,0,0,(1-1.05)/complex(0,1.05*1.05*0.03),1/complex(0,1.05*0.03),(1.05-1)/complex(0,1.05*0.03);
    0,0,0,(1.05-1)/complex(0,1.05*0.015),1/complex(0,1.05*0.015),(1-1.05)/complex(0,1.05*1.05*0.015),0,0,0,0,0,0,0,0,0;
    0,0,0,0,0,0,(1.05-1)/complex(0,1.05*0.03),1/complex(0,1.05*0.03),(1-1.05)/complex(0,1.05*1.05*0.03),0,0,0,0,0,0];

%虚数可以直接用2+i*4表示
%% 计算导纳阵Y
N=5;
Y=zeros(N,N);
for i=1:N
    for j=1:N
        if i==j
            Y(i,j)=sum(Y0(i,:))-Y0(i,3)-Y0(i,6)-Y0(i,9)-Y0(i,12)-Y0(i,15);
        else
            Y(i,j)=-Y0(i,j*3-1);
        end
    end
end

3、判断节点类型

%% 判断节点类型
for i=1:N 
    if data(i,2)==0 && data(i,3)==0                                        %如果PG、QG为0，那就是PQ节点
        fprintf('%d节点是PQ节点。\n',i);
        [p_q_derta(2*i-1),p_q_derta(2*i)]=PQ(p_q_matrix(2*i-1),p_q_matrix(2*i),Y,e_f_matrix,i);
    elseif isnan(data(i,2)) && isnan(data(i,3))                            %如果PG、QG为NAN，PL、QL为0，那就是V0节点
        fprintf('%d节点是V0节点。\n',i);
    else
        fprintf('%d节点是PV节点。\n',i);
        [p_q_derta(2*i-1),p_q_derta(2*i)]=PV(p_v_matrix(2*i-1),p_v_matrix(2*i),Y,e_f_matrix,i);
    end
end

4、计算雅可比矩阵并输出

%% 计算雅可比矩阵
for j=1:4
     for i=1:2*N
         syms (['e_f',num2str(i)]);
     end
    
     for i=1:N 
        if data(i,2)==0 && data(i,3)==0                                       
                p_derta_dd=p_q_matrix(2*i-1)-...
                             ['e_f',num2str(2*i-1)]*((real(Y(i,1))*e_f1-imag(Y(i,1))*e_f2)+...
                                                     (real(Y(i,2))*e_f3-imag(Y(i,2))*e_f4)+...
                                                     (real(Y(i,3))*e_f5-imag(Y(i,3))*e_f6)+...
                                                     (real(Y(i,4))*e_f7-imag(Y(i,4))*e_f8)+...
                                                     (real(Y(i,5))*e_f9-imag(Y(i,5))*e_f10))-...
                             ['e_f',num2str(2*i)]*((real(Y(i,1))*e_f2+imag(Y(i,1))*e_f1)+...
                                                   (real(Y(i,2))*e_f4+imag(Y(i,2))*e_f3)+...
                                                   (real(Y(i,3))*e_f6+imag(Y(i,3))*e_f5)+...
                                                   (real(Y(i,4))*e_f8+imag(Y(i,4))*e_f7)+...
                                                   (real(Y(i,5))*e_f10+imag(Y(i,5))*e_f9));
                J1=jacobian(p_derta_dd,[e_f1,e_f2,e_f3,e_f4,e_f5,e_f6,e_f7,e_f8,e_f9,e_f10]);
                
                q_derta_dd=p_q_matrix(2*i)-...
                              ['e_f',num2str(2*i)]*((real(Y(i,1))*e_f1-imag(Y(i,1))*e_f2)+...
                                                    (real(Y(i,2))*e_f3-imag(Y(i,2))*e_f4)+...
                                                    (real(Y(i,3))*e_f5-imag(Y(i,3))*e_f6)+...
                                                    (real(Y(i,4))*e_f7-imag(Y(i,4))*e_f8)+...
                                                    (real(Y(i,5))*e_f9-imag(Y(i,5))*e_f10))+...
                            ['e_f',num2str(2*i-1)]*((real(Y(i,1))*e_f2+imag(Y(i,1))*e_f1)+...
                                                    (real(Y(i,2))*e_f4+imag(Y(i,2))*e_f3)+...
                                                    (real(Y(i,3))*e_f6+imag(Y(i,3))*e_f5)+...
                                                    (real(Y(i,4))*e_f8+imag(Y(i,4))*e_f7)+...
                                                    (real(Y(i,5))*e_f10+imag(Y(i,5))*e_f9));
                J2=jacobian(q_derta_dd,[e_f1,e_f2,e_f3,e_f4,e_f5,e_f6,e_f7,e_f8,e_f9,e_f10]);
               J0=[J1;J2];
               if i==1
                   J0_0=J0;
               else
                   J0_0=[J0_0;J0];
               end
        elseif isnan(data(i,2)) && isnan(data(i,3))                          
            J0=zeros(2,10);
            if i==1
                J0_0=J0;
            else
                J0_0=[J0_0;J0];
            end
        else
                p_derta_dd=p_q_matrix(2*i-1)-...
                             ['e_f',num2str(2*i-1)]*((real(Y(i,1))*e_f1-imag(Y(i,1))*e_f2)+...
                                                     (real(Y(i,2))*e_f3-imag(Y(i,2))*e_f4)+...
                                                     (real(Y(i,3))*e_f5-imag(Y(i,3))*e_f6)+...
                                                     (real(Y(i,4))*e_f7-imag(Y(i,4))*e_f8)+...
                                                     (real(Y(i,5))*e_f9-imag(Y(i,5))*e_f10))-...
                             ['e_f',num2str(2*i)]*((real(Y(i,1))*e_f2+imag(Y(i,1))*e_f1)+...
                                                   (real(Y(i,2))*e_f4+imag(Y(i,2))*e_f3)+...