首页 > 编程语言> 文章详细

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

2020-11-25 17:00:29 阅读：265 来源： 互联网

标签：PQ end 导纳直角坐标拉夫 matlab PV Line Trans

使用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
    [dP,dQ,dU,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(Ps(1)+Ps(2)+Ps(3)+Ps(4)),imag(Ps(1)+Ps(2)+Ps(3)+Ps(4)))

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

%功率电压不平衡量
function [dP,dQ,dU,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

标签：PQ,end,导纳,直角坐标,拉夫,matlab,PV,Line,Trans
来源： https://blog.csdn.net/abcdesisjf/article/details/110123485

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

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

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

算例如下：

主程序

参考文献