使用matlab在直角坐标下使用牛顿拉夫逊算法计算潮流
程序员文章站
2022-03-07 11:19:18
本文是使用matlab在直角坐标系系统下,运用牛顿-拉夫逊算法进行的计算电力系统潮流,直角坐标系下进行运算的优点是易于理解,但缺点是矩阵比较庞大。...
使用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