newton-raphson power flow calculation.m

牛顿拉夫逊算法进行潮流计算

clear all
clc
p10=0.2;p20=-0.2;p30=-0.3;q10=0.2;q20=-0.1;q30=-0.4;count=0;
f0=[p10;p20;p30;q10;q20;q30]
syms a1 a2 a3 v1 v2 v3 ;

v= [v1;v2;v3];a=[a1;a2;a3];

p(1)= -5*v(1)*cos(a(1))+15*v(1)*sin(a(1))+13*v(1)^2-3.3333*v(1)*v(2)*cos(a(1)-a(2))+10*v(1)*v(2)*sin(a(1)-a(2))-5*v(1)*v(3)*cos(a(1)-a(3))+15*v(1)*v(3)*sin(a(1)-a(3));
p(2)=-5*v(2)*cos(a(2))+10*v(2)*sin(a(2))-3.3333*v(2)*v(1)*cos(a(2)-a(1))+10*v(2)*v(1)*sin(a(2)-a(1))+8.3333*v(2)^2;
p(3)=-10*v(3)*cos(a(3))+30*v(3)*sin(a(3))-5*v(3)*v(1)*cos(a(3)-a(1))+15*v(3)*v(1)*sin(a(3)-a(1))+15*v(3)^2;
q(1)=-5*v(1)*sin(a(1))-15*v(1)*cos(a(1))+39.84*v(1)^2-3.3333*v(1)*v(2)*sin(a(1)-a(2))-10*v(1)*v(2)*cos(a(1)-a(2))-5*v(1)*v(3)*sin(a(1)-a(3))-15*v(1)*v(3)*cos(a(1)-a(3));
q(2)=-5*v(2)*sin(a(2))-10*v(2)*cos(a(2))-3.3333*v(2)*v(1)*sin(a(2)-a(1))-10*v(2)*v(1)*cos(a(2)-a(1))+19.89*v(2)^2;
q(3)=-10*v(3)*sin(a(3))-30*v(3)*cos(a(3))-5*v(3)*v(1)*sin(a(3)-a(1))-15*v(3)*v(1)*cos(a(3)-a(1))+44.89*v(3)^2;
x=[v;a];
f=[p,q];
Jcb = jacobian(f,x);
df=f0-f';
 a1=0;a2=0;a3=0;
 v1=1;v2=1;v3=1;
 while(count<3)
 F=subs(f)
 dF=subs(df)
 J=subs(Jcb)
 invJ=inv(subs(Jcb))
 dx=invJ*dF
 a1=a1+dx(1);
 a2=a2+dx(2);
 a3=a3+dx(3);
 v1=v1+dx(4);
 v2=v2+dx(5);
 v3=v3+dx(6);
 xx=subs(x)
 count=count+1;
end


%     
%     count=count+1;
% 	dd=-inv(subs(Jcb))*subs(w);
%     dx=dd(1);
%     dy=dd(2);
%     x=x+dx;
%     y=y+dy;
%     
% 
% 
% x = x
% y = y

% Ans are as follows:
% x = -8.5613e-023
% y = 6.3909e-023