MATLAB潮流计算仿真设计14892.pdf

《MATLAB潮流计算仿真设计14892.pdf》由会员分享，可在线阅读，更多相关《MATLAB潮流计算仿真设计14892.pdf（8页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、.附录A MATLAB程序%本程序的功能是用牛顿拉夫逊法进行潮流计算%B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳%5、支路的变比；6、支路首端处于K侧为1，1侧为0%B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值%4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量%6、节点分类标号 clear;n=10;%input(请输入节点数：n=);nl=11;%input(请输入支路数：nl=);isb=1;%input(请输入平衡母线节点号：isb=);pr=0.00001;%input(请输入误差精度：pr=);B1=1 2 1.755e

2、-2+4.155e-2i 0.26i 1 0;1 4 3.159e-2+7.479e-2i 0.1215i 1 0;1 6 3.159e-2+7.479e-2i 0.1215i 1 0;2 3 3.68e-3+0.11135i 0 0.909 1;4 5 3.68e-3+0.11135i 0 0.909 1;4 6 2.808e-2+6.648e-2i 0.108i 1 0;6 7 3.0865e-3+0.0833i 0 0.909 1;6 8 3.159e-2+7.479e-2i 0.1215i 1 0;6 10 2.457e-2+5.817e-2i 0.0945i 1 0;8 9 3.08

3、65e-3+0.0833i 0 0.909 1;8 10 2.808e-2+6.648e-2i 0.108i 1 0;%input(请输入由支路参数形成的矩阵：B1=);B2=0 0 1.05 1.05 0 1;0 0 1 0 0 2;0 0.6+0.3718i 1 0 0 2;0 0 1 0 0 2;0 0.4+0.247i 1 0 0 2;0 0 1 0 0 2;0 0.35+0.2169i 1 0 0 2;0 0 1 0 0 2;0 0.5+0.3099i 1 0 0 2;0.8 0 1.05 1.05 0 3;%input(请输入各节点参数形成的矩阵：B2=);Y=zeros(n);e

4、=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);%修改部分 ym=0;SB=100;UB=220;%ym=input(您输入的参数是标么值？（若不是则输入一个不为零的数值）);if ym=0%SB=input(请输入功率基准值:SB=);.%UB=input(请输入电压基准值:UB=);YB=SB./UB./UB;BB1=B1;BB2=B2;for i=1:nl B1(i,3)=B1(i,3)*YB;B1(i,4)=B1(i,4)./YB;end disp(B1矩阵B1=);disp(B1)for i=1:n

5、 B2(i,1)=B2(i,1)./SB;B2(i,2)=B2(i,2)./SB;B2(i,3)=B2(i,3)./UB;B2(i,4)=B2(i,4)./UB;B2(i,5)=B2(i,5)./SB;end disp(B2矩阵B2=);disp(B2)end%-for i=1:nl%支路数 if B1(i,6)=0%左节点处于低压侧 p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);end Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5);%非对角元 Y(q,p)=Y(p,q);Y(q,q)=Y(q,q)+1./(B1(i,3)*B1

6、(i,5)2)+B1(i,4)./2;%对角元K侧 Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;%对角元1侧 end%求导纳矩阵 disp(导纳矩阵 Y=);disp(Y)%-G=real(Y);B=imag(Y);%分解出导纳阵的实部和虚部 for i=1:n%给定各节点初始电压的实部和虚部 e(i)=real(B2(i,3);f(i)=imag(B2(i,3);V(i)=B2(i,4);%PV节点电压给定模值 end for i=1:n%给定各节点注入功率 S(i)=B2(i,1)-B2(i,2);%i节点注入功率SG-SL .B(i,i)=B(i,i)+B2(

7、i,5);%i节点无功补偿量 end%=P=real(S);Q=imag(S);ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;while IT2=0 IT2=0;a=a+1;for i=1:n if i=isb%非平衡节点 C(i)=0;D(i)=0;for j1=1:n C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%(Gij*ej-Bij*fj)D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%(Gij*fj+Bij*ej)end P1=C(i)*e(i)+f(i)*D(i);%节 点 功 率 P 计 算 ei (G

8、ij*ej-Bij*fj)+fi(Gij*fj+Bij*ej)Q1=C(i)*f(i)-e(i)*D(i);%节 点 功 率 Q 计 算 fi (Gij*ej-Bij*fj)-ei(Gij*fj+Bij*ej)%求P,Q V2=e(i)2+f(i)2;%电压模平方%=以下针对非PV节点来求取功率差及Jacobi矩阵元素=if B2(i,6)=3%非PV节点 DP=P(i)-P1;%节点有功功率差 DQ=Q(i)-Q1;%节点无功功率差%=以上为除平衡节点外其它节点的功率计算=%=求取Jacobi矩阵=for j1=1:n if j1=isb&j1=i%非平衡节点&非对角元 X1=-G(i,j1

9、)*e(i)-B(i,j1)*f(i);%dP/de=-dQ/df X2=B(i,j1)*e(i)-G(i,j1)*f(i);%dP/df=dQ/de X3=X2;%X2=dp/df X3=dQ/de X4=-X1;%X1=dP/de X4=dQ/df p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;elseif j1=i&j1=isb%非平衡节点&对角元 X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%dP/de X2=-D(i)+B(i,i)

10、*e(i)-G(i,i)*f(i);%dP/df X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);%dQ/de X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);%dQ/df p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列Q m=p+1;J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列P.J(m,q)=X2;end end else%=下面是针对PV节点来求取Jacobi矩阵的元素=DP=P(i)-P1;%PV节点有功误差 DV=V(i)2-V2;%PV节点电压误差 for j1=1:n if

11、j1=isb&j1=i%非平衡节点&非对角元 X1=-G(i,j1)*e(i)-B(i,j1)*f(i);%dP/de X2=B(i,j1)*e(i)-G(i,j1)*f(i);%dP/df X5=0;X6=0;p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;elseif j1=i&j1=isb%非平衡节点&对角元 X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%dP/de X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%dP

12、/df X5=-2*e(i);X6=-2*f(i);p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;end end end end end%=以上为求雅可比矩阵的各个元素=for k=3:N0%N0=2*n（从第三行开始，第一、二行是平衡节点）k1=k+1;N1=N;%N=N0+1 即 N=2*n+1扩展列P、Q for k2=k1:N1%扩展列P、Q J(k,k2)=J(k,k2)./J(k,k);%非对角元规格化 end J(k,k)=1;%对角元规格化 if

13、k=3%不是第三行%=k4=k-1;for k3=3:k4%用k3行从第三行开始到当前行前的k4行消去 for k2=k1:N1%k3行后各行下三角元素 J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算.end J(k3,k)=0;end if k=N0 break;end%=for k3=k1:N0 for k2=k1:N1 J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算 end J(k3,k)=0;end else for k3=k1:N0 for k2=k1:N1 J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k

14、2);%消去运算 end J(k3,k)=0;end end end%=上面是用线性变换方式将Jacobi矩阵化成单位矩阵=for k=3:2:N0-1 L=(k+1)./2;e(L)=e(L)-J(k,N);%修改节点电压实部 k1=k+1;f(L)=f(L)-J(k1,N);%修改节点电压虚部 end%-修改节点电压-for k=3:N0 DET=abs(J(k,N);if DET=pr%电压偏差量是否满足要求 IT2=IT2+1;%不满足要求的节点数加1 end end ICT2(a)=IT2;ICT1=ICT1+1;end%用高斯消去法解w=-J*V disp(迭代次数：);disp(

15、ICT1);disp(没有达到精度要求的个数：);disp(ICT2);.for k=1:n V(k)=sqrt(e(k)2+f(k)2);sida(k)=atan(f(k)./e(k)*180./pi;E(k)=e(k)+f(k)*j;end%=计算各输出量=disp(各节点的实际电压标幺值E为(节点号从小到大排列)：);disp(E);EE=E*UB;disp(EE);disp(-);disp(各节点的电压大小V为(节点号从小到大排列)：);disp(V);VV=V*UB;disp(VV);disp(-);disp(各节点的电压相角sida为(节点号从小到大排列)：);disp(sida)

16、;for p=1:n C(p)=0;for q=1:n C(p)=C(p)+conj(Y(p,q)*conj(E(q);end S(p)=E(p)*C(p);end disp(各节点的功率S为(节点号从小到大排列)：);disp(S);disp();SS=S*SB;disp(SS);disp(-);disp(各条支路的首端功率Si为(顺序同您输入B1时一致)：);for i=1:nl p=B1(i,1);q=B1(i,2);if B1(i,6)=0 Si(p,q)=E(p)*(conj(E(p)*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5)-conj(E(q)*con

17、j(1./(B1(i,3)*B1(i,5);Siz(i)=Si(p,q);else Si(p,q)=E(p)*(conj(E(p)*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5)-conj(E(q)*conj(1./(B1(i,3)*B1(i,5);Siz(i)=Si(p,q);.end disp(Si(p,q);SSi(p,q)=Si(p,q)*SB;ZF=S(,num2str(p),num2str(q),)=,num2str(SSi(p,q);disp(ZF);%disp(SSi(p,q);disp(-);end disp(各条支路的末端功率Sj为(顺序同您输入B

18、1时一致)：);for i=1:nl p=B1(i,1);q=B1(i,2);if B1(i,6)=0 Sj(q,p)=E(q)*(conj(E(q)*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5)-conj(E(p)*conj(1./(B1(i,3)*B1(i,5);Sjy(i)=Sj(q,p);else Sj(q,p)=E(q)*(conj(E(q)*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5)-conj(E(p)*conj(1./(B1(i,3)*B1(i,5);Sjy(i)=Sj(q,p);end disp(Sj(q,p);SSj(q

19、,p)=Sj(q,p)*SB;ZF=S(,num2str(q),num2str(p),)=,num2str(SSj(q,p);disp(ZF);%disp(SSj(q,p);disp(-);end disp(各条支路的功率损耗DS为(顺序同您输入B1时一致)：);for i=1:nl p=B1(i,1);q=B1(i,2);DS(i)=Si(p,q)+Sj(q,p);disp(DS(i);DDS(i)=DS(i)*SB;ZF=DS(,num2str(p),num2str(q),)=,num2str(DDS(i);disp(ZF);%disp(DDS(i);disp(-);end figure(1);subplot(2,2,1);plot(V);xlabel(节点号);ylabel(电压标幺值);.grid on;subplot(2,2,2);plot(sida);xlabel(节点号);ylabel(电压角度);grid on;subplot(2,2,3);bar(real(S);ylabel(节点注入有功);grid on;subplot(2,2,4);bar(Siz);ylabel(支路首端无功);grid on;