电力系统分析大作业-

1、电力系统分析大作业电力系统分析程序设计 专业班级 姓 名 学 号 日 期 要求计算所给系统的潮流，设发电机G1的端电压为1 p.u.，发出的有功、无功可调；发电机G2的端电压为1 p.u.，按指定的有功P=0.5 p.u.发电，取=10-4。SB=100MVA，UB=Uav。一牛顿拉夫逊法（1） 程序代码 clear; clc; %节点编号，5号为平衡节点，4号为PV节点，1、2、3号为PQ节点 %输入原始数据，求节点导纳矩阵 y=0; y(1,1) = 1/(0.03+0.1i)+1/(0.025+0.08i)+0.09i+0.07i; y(1,2) = -1/(0.025+0.08i);

2、y(1,3) = -1/(0.03+0.1i); y(2,1) = -1/(0.025+0.08i); y(2,2) = 1/(0.02+0.06i)+1/(0.025+0.08i)+0.05i+0.07i+1/1.05222*1/(0.1905i); y(2,3) = -1/(0.02+0.06i); y(2,4) = 4.9889i; y(3,1) = -1/(0.03+0.1i); y(3,2) = -1/(0.02+0.06i); y(3,3) = 7.7523-28.7757i; y(3,5) = 4.9889i; y(4,2) = 4.9889i; y(4,4) = -5.2493

3、i; y(5,3) = 4.9889i; y(5,5) = -5.2493i; Y=y; display('导纳矩阵为: '); Y %Y 为导纳矩阵 G=real(Y); B=imag(Y); %原始节点功率 Pg=zeros(1,5); Qg=zeros(1,5); Pd=zeros(1,5); Qd=zeros(1,5); S(1)=-0.8055-0.5320i; S(2)=-0.18-0.12i; S(3)=0; S(4)=0.5; S(5)=0; P=real(S); Q=imag(S); %赋初值 U = ones(1,5); a = zeros(1,5); dx

4、=ones(7,1); fx=ones(7,1); count=0; %计算迭代次数 while (max(fx)>1e-4) %收敛判据 for i=1:4 for j=1:4 H(i,j)=0; dP(i)=0; end end for i=1:3 for j=1:3 L(i,j)=0; dQ(i)=0; end end for i=1:4 for j=1:3 N(i,j)=0; end end for i=1:3 for j=1:4 K(i,j)=0; end end %P %失配功率表达式 for i=1:4 x=0; for j=1:5 %dP(i)=dP(i)-U(i)*U(

5、j)*(G(i,j)*cos(a(i)-a(j) + B(i,j)*sin(a(i)-a(j) ); x=x+U(i)*U(j)*(G(i,j)*cos(a(i)-a(j)+B(i,j)*sin(a(i)-a(j); end dP(i)=P(i)-x; end for i=1:3 y=0; for j=1:5 y=y+(U(i)*U(j)*(G(i,j)*sin(a(i)-a(j)-B(i,j)*cos(a(i)-a(j); end dQ(i)=Q(i)-y; end fx=dP,dQ' %求雅克比矩阵 %当i=j时候求H,N,K,L 如下： for i=1:4 for j=1:5 i

6、f i=j H(i,i)=H(i,i)+U(i)*U(j)*(G(i,j)*sin(a(i)-a(j)-B(i,j)*cos(a(i)-a(j) ); end end end for i=1:4 for j=1:4 if i=j H(i,j)=-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j)-B(i,j)*cos(a(i)-a(j); end end end for i=1:3 x=0; for j=1:5 if i=j x=x + (- (U(j)*( G(i,j)*cos(a(i)+a(j) + B(i,j)*sin(a(i)-a(j) ) ); end end N(i,i

7、)=-2*U(i)*G(i,i)+x; end for i=1:4 for j=1:3 if i=j N(i,j)=-U(i)*( G(i, j)*cos(a(i)-a(j) + B(i,j)*sin(a(i)-a(j) ); end end end for i=1:3 for j=1:5 if i=j K(i,i)=K(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j) + B(i,j)*sin(a(i)-a(j) ); end end end for i=1:3 for j=1:4 if i=j K(i,j)=U(i)*U(j)*( G(i,j)*cos(a(i)-a

8、(j) + B(i,j)*sin(a(i)-a(j) ); end end end for i=1:3 x=0; for j=1:5 if i=j x=x+(-( U(j)*(G(i,j)*sin(a(i)-a(j) - B(i,j)*cos(a(i)-a(j); end end L(i,i)=2*U(i)*B(i,i)+x; end for i=1:3 for j=1:3 if i=j L(i,j)=-U(i)*( G(i,j)*sin(a(i)-a(j) - B(i,j)*cos(a(i)-a(j) ); end end end J=H,N;K,L; %J 为雅克比矩阵 dx=-(inv(

9、J)*fx); for i=1:4 da(i)=dx(i); end for i=1:3 dU(i)=dx(i+4); end for i=1:4 a(i)=a(i)+da(i); end for i=1:3 U(i)=U(i)+dU(i); end count=count+1; end %求节点注入的净功率 i=5; for j=1:5 P(i)=P(i)+U(i)*U(j)*(G(i,j)*cos(a(i)-a(j)+B(i,j)*sin(a(i)-a(j); Q(i)=Q(i)+U(i)*U(j)*(G(i,j)*sin(a(i)-a(j)-B(i,j)*cos(a(i)-a(j); e

10、ndS(5)=P(5)+Q(5)*sqrt(-1); if count=3 display('最终结果为: '); H,N,K,L display('雅可比矩阵如下: '); Jdisp('/*/') a./pi*180 disp('/*/') fx,dx,U,adisplay('迭代次数为: ');countdisplay('节点注入净功率如下: '); S%求节点注入电流display('节电注入电流如下: '); I=Ydisplay('迭代次数为: ');co

11、untend 运算结果如下：导纳矩阵为: 最终结果为: H = -20.6021 11.3894 9.2127 0 11.5923 -32.2398 15.5947 5.0528 9.3958 15.6419 -30.1181 0 0 5.0528 0 -5.0528N = -5.6127 3.8473 3.0267 3.2930 -8.7264 5.1580 2.5065 5.0456 -8.0992 0 -0.4914 0K = 7.0121 -3.9156 -3.0965 0 -3.2661 9.0432 -5.2770 -0.5001 -2.4861 -5.1352 8.1121 0L

12、 = -19.6995 11.1907 9.0050 11.6878 -31.4424 15.2431 9.4732 15.3690 -29.4399雅可比矩阵如下: /*/ans = -7.4732 -5.8532 -5.5850 -0.2006 0/*/fx = -0.0017 -0.0021 -0.0020 -0.0001 -0.0003 -0.0004 -0.0004dx = -0.0013 -0.0013 -0.0012 -0.0013 -0.0003 -0.0002 -0.0002U = 0.9916 1.0175 1.0229 1.0000 1.0000a = -0.1304 -

13、0.1022 -0.0975 -0.0035 0迭代次数为: count = 3节点注入净功率如下: S = -0.8055-0.5320i -0.1800-0.1200i 0 0.5000 0.4966 + 0.1706i节电注入电流如下:迭代次数为: count = 3二PQ解耦迭代程序代码n=5;m=3;% 各参数初始化 Y=zeros(n,n); %导纳矩阵U=ones(n,1); %电压矢量a=zeros(n,1); %相角矢量Ps=zeros(n,1);Qs=zeros(n,1);P=zeros(n,1);Q=zeros(n,1);p=zeros(n-1,1);q=zeros(m,

14、1);aa=zeros(n-1,1);u=zeros(m,1);k=0;%节点导纳矩阵的生成I=sqrt(-1); %由于后面用到ij作为循环变量，故重定义虚数单位Y=6.311-I*20.3722 -2.7523+I*9.1743 -3.5587+I*11.3879 0 0; -3.5587+I*11.3879 8.5587-I*30.9993 -5+I*15 I*4.9889 0; -2.7523+I*9.1743 -5+I*15 7.7523-I*28.7557 0 I*4.9889; 0 I*4.9889 0 -I*5.2493 0; 0 0 I*4.9889 0 -I*5.2493;

15、%初值的设定U=1;1;1;1;1;a=0;0;0;0;0;Ps=-0.8055;-0.18;0;0.5;0;Qs=-0.532;-0.12;0;0;0;%求失配功率P(k)和Q(k) for i=1:n-1 s=0; for j=1:n s=s+U(j)*(real(Y(i,j)*cos(a(i)-a(j)+imag(Y(i,j)*sin(a(i)-a(j); end P(i)=U(i)*s;endfor i=1:n-1 s=0; for j=1:n s=s+U(j)*(real(Y(i,j)*sin(a(i)-a(j)-imag(Y(i,j)*cos(a(i)-a(j); end Q(i)

16、=U(i)*s;endfor i=1:n-1 detp(i)=Ps(i)-P(i);endfor i=1:m detq(i)=Qs(i)-Q(i);end %生成B'和B"矩阵 Bp=zeros(n-1,n-1); Bpp=zeros(m,m); for i=1:n-1; for j=1:n-1; Bp(i,j)=imag(Y(i,j); end; end; for i=1:m; for j=1:m; Bpp(i,j)=imag(Y(i,j); end; end;k=0;while (max(detp)>=0.0001 | max(detq)>=0.0001) %

17、设定精确度 aa=(-1)*inv(Bp)*detp(1)/U(1);detp(2)/U(2);detp(3)/U(3);detp(4)/U(4); u=(-1)*inv(Bpp)*detq(1)/U(1);detq(2)/U(2);detq(3)/U(3); %求解修正方程，得到修正量 aa,u for i=1:n-1 %计算 a(节点电压相角),U(节点电压大小) a(i)=a(i)+aa(i); end for i=1:m U(i)=U(i)+u(i); end% detp'detq' %迭代过程中失配功率变化情况 k=k+1;%(记录迭代次数) for i=1:n-1

18、%计算功率误差 p,q s=0; for j=1:n s=s+U(j)*(real(Y(i,j)*cos(a(i)-a(j)+imag(Y(i,j)*sin(a(i)-a(j); end P(i)=U(i)*s; end for i=1:n-1 s=0; for j=1:n s=s+U(j)*(real(Y(i,j)*sin(a(i)-a(j)-imag(Y(i,j)*cos(a(i)-a(j); end Q(i)=U(i)*s; end for i=1:n-1 detp(i)=Ps(i)-P(i); end for i=1:m detq(i)=Qs(i)-Q(i); endendA=a*180/pi; %对相角进行转换（由弧度制装换为角度）display('最终的计算结果如下: ');AUdisplay('生成的节点导纳矩阵为：');Ydisplay('算得矩阵B为：');Bpdisplay('算得矩阵B”为：');Bppdisplay('各节点的电压为:');for i=1:n fprintf('节点%d: %d(%d)°nn',i,U(i),A(i);enddisplay('失配功率为: &