使用matlab在直角坐标下使用牛顿拉夫逊算法计算潮流

小编这周刚进行了课程设计，根据何仰赞的《电力系统分析》的内容编写程序，写完后觉得分享给大家会更有意义，这是小编第一次发博客，有不妥之处还请大家包涵，同时也欢迎大家纠错。

在这里用一个框图简单给大家解释一下，牛拉法进行潮流运算到底怎么回事：

本人使用的变压器模型如下

算例如下：

主程序

clear %清除变量
clc %清屏
%数据输入（标幺值）
%交流线参数：I 侧母线，J侧母线，阻抗，1/2接地导纳
Line=[1 4 0.12+0.5i 0.01920i;
      4 2 0.08+0.4i 0.01413i;
      2 3 0.10+0.4i 0.01528i];
%变压器参数：I侧母线，J侧母线，阻抗，变比（需要折算到i侧）
Trans=[3 1 0.3i 0.909];
n=4;%四个节点
PQ=2;
PV=1;
%按照标号的节点参数 P,Q,U,;
P=[-0.30 -0.55 -0.5];
Q=[-0.18 -0.13 0];
U=[0 0 1.10 1.05];
%变压器Π型等效导纳电路
Yt=zeros(size(Trans,1),3);
Yt(:,1)=Trans(:,4)./Trans(:,3);%互导纳
Yt(:,2)=(1-Trans(:,4))./Trans(:,3);%i侧对地导纳
Yt(:,3)=Trans(:,4).*(Trans(:,4)-1)./Trans(:,3);%j侧对地导纳
%I侧母线 J侧母线 互导纳 I侧自导纳 J侧自导纳
Trans_pi=[Trans(:,1:2) Yt(:,1) Yt(:,2) Yt(:,3)];
%节点导纳矩阵运算
Y=zeros(n);
%计算导线间导纳
for k=1:size(Line,1)
    i=Line(k,1); j=Line(k,2);
    Y(i,j)=-1/Line(k,3);
    Y(j,i)=Y(i,j);
    Y(i,i)=Y(i,i)+Line(k,4)+1/Line(k,3);
    Y(j,j)=Y(j,j)+Line(k,4)+1/Line(k,3);
end
%计算变压器的导纳
for k=1:size(Trans,1)
    i=Trans(k,1);j=Trans(k,2);
    Y(i,j)=-Trans_pi(k,3);
    Y(j,i)=Y(i,j);
    Y(i,i)=Y(i,i)+Trans_pi(k,4)+Trans_pi(k,3);
    Y(j,j)=Y(j,j)+Trans_pi(k,5)+Trans_pi(k,3);
end
%提取感抗矩阵的实部和虚部
G=real(Y);
B=imag(Y);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算潮流
%设出初值
e=[1 1 1.10 1.05];
f=[0 0 0 0];
for k=1:100
    [dW]=Unbalanced(n,PQ,PV,P,Q,U,e,f,G,B);
    if max(abs(dW))<0.00001
        fprintf('迭代%d次结束\n',k-1);
        break;
    end
    [J]=Jacobi(n,PQ,PV,e,f,G,B);
    dV=(-J)\dW';
    for i=1:n-1
        e(i)=e(i)+dV(2*i-1);
        f(i)=f(i)+dV(2*i);
    end
    %创建cell数组存取迭代过程中全部J,e,f,dW的矩阵,Jall代表所有的J
    Jall{k}=J;
    dWall{k}=dW;
    eall{k}=e;
    fall{k}=f;
    
end

%依次打印每次迭代的各项数据
for i=1:k-1
    Usall{i}=eall{i}+fall{i}*1i;
end

for i=1:k-1
    fprintf('第%d次迭代节点电压\n',i);
    disp(Usall{i})
    fprintf('%c',8)
    fprintf('第%d次迭代节点不平衡量\n',i);
    disp(dWall{i})  
end
for i=1:k-1
    fprintf('第%d次迭代Jacobi矩阵\n',i);
    disp(Jall{i})
end
Is=zeros(n);
Us=e+f*1i;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%求最终每个节点的功率
for i=1:n
    for j=1:n
        Is(i)=Is(i)+conj(Y(i,j))*conj(Us(j));
    end
    %每个节点的功率
    Ps(i)=Us(i)*Is(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%求支路电流
for k=1:size(Line,1)
    i=Line(k,1); j=Line(k,2);
    I(i,j)=conj((abs(Us(i))^2*conj(Line(k,4))+Us(i)*(conj(Us(i))-conj(Us(j)))*conj(1/Line(k,3)))/Us(i));
    I(j,i)=conj((abs(Us(j))^2*conj(Line(k,4))+Us(j)*(conj(Us(j))-conj(Us(i)))*conj(1/Line(k,3)))/Us(j));
end
for k=1:size(Trans,1)
    i=Trans_pi(k,1); j=Trans_pi(k,2);
    I(i,j)=conj((abs(Us(i))^2*conj(Trans_pi(k,4))+Us(i)*(conj(Us(i))-conj(Us(j)))*conj(Trans_pi(k,3)))/Us(i));
    I(j,i)=conj((abs(Us(j))^2*conj(Trans_pi(k,5))+Us(j)*(conj(Us(j))-conj(Us(i)))*conj(Trans_pi(k,3)))/Us(j));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%求线路功率
for k=1:size(Line,1)
    i=Line(k,1); j=Line(k,2);
    S(i,j)=abs(Us(i))^2*conj(Line(k,4))+Us(i)*(conj(Us(i))-conj(Us(j)))*conj(1/Line(k,3));
    S(j,i)=abs(Us(j))^2*conj(Line(k,4))+Us(j)*(conj(Us(j))-conj(Us(i)))*conj(1/Line(k,3));
end
%%%%%%%%%
%以下是两个求变压器线路功率的办法，取一个即可
for k=1:size(Trans_pi,1)
    i=Trans_pi(k,1); j=Trans_pi(k,2);
    S(i,j)=abs(Us(i))^2*conj(Trans_pi(k,4))+Us(i)*(conj(Us(i))-conj(Us(j)))*conj(Trans_pi(k,3));
    S(j,i)=abs(Us(j))^2*conj(Trans_pi(k,5))+Us(j)*(conj(Us(j))-conj(Us(i)))*conj(Trans_pi(k,3));
end

%for k=1:size(Trans,1)
%   i=Trans(k,1);j=Trans(k,2);
%   S(i,j)=Us(i)*(conj(Us(i))-conj(Us(j))*Trans(k,4))*conj(1/Trans(3));
%   S(j,i)=Us(j)*Trans(k,4)*(conj(Us(j))*Trans(k,4)-conj(Us(i)))*conj(1/Trans(3));
%end
%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('节点电压为\n')
for i=1:n
    if(real(Us(i))||imag(Us(i))~=0)
        fprintf('V%d=%fj%f\n',i,real(Us(i)),imag(Us(i)));
    end
end
fprintf('节支路电流为\n')
for i=1:n
    for j=1:n
        if(real(I(i,j))||imag(I(i,j))~=0)
            fprintf('I%d%d=%f+j%f\n',i,j,real(I(i,j)),imag(I(i,j)));
        end
        
    end
end
fprintf('最终每个线路的功率\n')
for i=1:n
    for j=1:n
        if(real(S(i,j))||imag(S(i,j))~=0)
            fprintf('S%d%d=%f+j%f\n',i,j,real(S(i,j)),imag(S(i,j)));
        end
    end
end
fprintf('平衡节点功率\n')
fprintf('P%d+jQ%d=%f+j%f\n',n,n,real(Ps(n)),imag(Ps(n)));
fprintf('网络损耗\n%f+j%f\n',real(sum(Ps)),imag(sum(Ps)))

下面是不平衡量（Unbalanced）程序

%功率电压不平衡量
function [dW]=Unbalanced(n,PQ,PV,P,Q,U,e,f,G,B)
for i=1:n-1
    dP(i)=P(i);
end
for i=1:PQ
    dQ(i)=Q(i);
end
for i=1:n-1
    for j=1:n
        dP(i)=dP(i)-(e(i)*(e(j)*G(i,j)-f(j)*B(i,j))+f(i)*(f(j)*G(i,j)+e(j)*B(i,j)));
    end
end
for i=1:PQ
    for j=1:n
        dQ(i)=dQ(i)-(f(i)*(e(j)*G(i,j)-f(j)*B(i,j))-e(i)*(f(j)*G(i,j)+e(j)*B(i,j)));
    end
end
for i=PQ+1:PQ+PV 
    dU(i)=U(i)^2-e(i)^2-f(i)^2;
end
%构建dw
k=1;
for i=1:PQ+PV
    dW(k)=dP(i);
    k=k+2;
end
k=2;
for i=1:PQ
    dW(k)=dQ(i);
    k=k+2;
end
k=2*PQ+2;
for i=PQ+1:PQ+PV
    dW(k)=dU(i);
end

以下是雅可比矩阵（J）

function[J]=Jacobi(n,PQ,PV,e,f,G,B)
J=zeros(2*n-2);
for i=1:PQ
    for j=1:n-1
        J(2*i-1,2*j-1)=-(G(i,j)*e(i)+B(i,j)*f(i));
        J(2*i-1,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
        J(2*i,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i);
        J(2*i,2*j)=G(i,j)*e(i)+B(i,j)*f(i);
    end
    for k=1:n
        J(2*i-1,2*i-1)=J(2*i-1,2*i-1)-(G(i,k)*e(k)-B(i,k)*f(k));
        J(2*i-1,2*i)=J(2*i-1,2*i)-(G(i,k)*f(k)+B(i,k)*e(k));
        J(2*i,2*i-1)=J(2*i,2*i-1)+G(i,k)*f(k)+B(i,k)*e(k);
        J(2*i,2*i)=J(2*i,2*i)-(G(i,k)*e(k)-B(i,k)*f(k));
    end
end
for i=PQ+1:PQ+PV
    for j=1:n-1
        J(2*i-1,2*j-1)=-(G(i,j)*e(i)+B(i,j)*f(i));
        J(2*i-1,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
        J(2*i,2*j-1)=0;
        J(2*i,2*j)=0;
    end
    for k=1:n
        J(2*i-1,2*i-1)=J(2*i-1,2*i-1)-(G(i,k)*e(k)-B(i,k)*f(k));
        J(2*i-1,2*i)=J(2*i-1,2*i)-(G(i,k)*f(k)+B(i,k)*e(k));
    end
    J(2*i,2*i-1)=-2*e(i);
    J(2*i,2*i)=-2*f(i);
end
end

参考文献

【1】《电力系统分析（第四版）》 何仰赞 温增银，华中科技大学出版社，2017

写在最后

各位同学在大学生活中除了建模以外应该很少用到matlab，所以小编写这篇文章的初衷就是帮助想申优的同学加深思考，而不是被程序束缚，小编在写的时候也是边写边查，现学现卖，matlab在许多领域都有重要应用，希望大家能熟练使用。最后祝大家可以用这些代码写出属于自己的思想！
最后附上程序文件
链接: 文件.提取码0104

本文地址：https://blog.csdn.net/abcdesisjf/article/details/110123485