《牛拉法潮流计算(共10页).docx》由会员分享，可在线阅读，更多相关《牛拉法潮流计算(共10页).docx（10页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上%本程序的功能是用牛拉法进行潮流计算%原理介绍详见鞠平著电气工程%默认数据为鞠平著电气工程例8.4所示数据%B1是支路参数矩阵%第一列和第二列是节点编号。节点编号由小到大编写%对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点编号%第三列为支路的串列阻抗参数，含变压器支路此值为变压器短路电抗%第四列为支路的对地导纳参数，含变压器支路此值不代入计算%第五烈为含变压器支路的变压器的变比，变压器非标准电压比%第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，“0”为不含有变压器%B2为节点参数矩阵%第一列为节点注入发电功率参数%第二列为节点负荷功
2、率参数%第三列为节点电压参数%第四列%第五列%第六列为节点类型参数，“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数%X为节点号和对地参数矩阵%第一列为节点编号%第二列为节点对地参数clear;clc;num=input(是否采用默认数据？(1-默认数据;2-手动输入);if num=1 n=4; n1=4; isb=4; pr=0.00001; B1=1 2 0.1667i 0 0.8864 1;1 3 0.1302+0.2479i 0.0258i 1 0;1 4 0.1736+0.3306i 0.0344i 1 0;3 4 0.2603+0.4959i 0.0518i 1 0; B
3、2=0 0 1 0 0 2;0 -0.5-0.3i 1 0 0 2;0.2 0 1.05 0 0 3;0 -0.15-0.1i 1.05 0 0 1; X=1 0;2 0.05i;3 0;4 0;else n=input(请输入节点数:n=); n1=input(请输入支路数:n1=); isb=input(请输入平衡节点号:isb=); pr=input(请输入误差精度:pr=); B1=input(请输入支路参数:B1=); B2=input(请输入节点参数:B2=); X=input(节点号和对地参数:X=);endTimes=1; %迭代次数%创建节点导纳矩阵Y=zeros(n);fo
4、r i=1:n1 if B1(i,6)=0 %不含变压器的支路 p=B1(i,1); q=B1(i,2); Y(p,q)=Y(p,q)-1/B1(i,3); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4); Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4); else %含有变压器的支路 p=B1(i,1); q=B1(i,2); Y(p,q)=Y(p,q)-B1(i,5)/B1(i,3); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+B1(i,5)/B1(i,3)+(1-B1(i,5)/B1(i,3);
5、 Y(q,q)=Y(q,q)+B1(i,5)/B1(i,3)+(B1(i,5)*(B1(i,5)-1)/B1(i,3); endendfor i=1:n1 Y(i,i)=Y(i,i)+X(i,2); %计及补偿电容电纳enddisp(导纳矩阵为：);disp(Y); %显示导纳矩阵 %初始化OrgS、DetaSOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1); %创建OrgS，用于存储初始功率参数h=0;j=0;for i=1:n %对PQ节点的处理 if i=isb&B2(i,6)=2 %不是平衡点&是PQ点 h=h+1; for j=1:n %公式8-74
6、%Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej) %Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej) OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1) =OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(
7、j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end endendfor i=1:n %对PV节点的处理，注意这时不可再将h初始化为0 if i=isb&B2(i,6)=3 %不是平衡点&是PV点 h=h+1; for j=1:n %公式8-75-a %Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej) %Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej) OrgS(2*h-1,1)=OrgS
8、(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end ende
9、nd%创建PVU 用于存储PV节点的初始电压PVU=zeros(n-h-1,1);t=0;for i=1:n if B2(i,6)=3 t=t+1; PVU(t,1)=B2(i,3); endend%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:n %对PQ节点的处理 if i=isb&B2(i,6)=2 h=h+1; DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); %delPi DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1); %delQi endendt=0;for i=1:n %对P
10、V节点的处理，注意这时不可再将h初始化为0 if i=isb&B2(i,6)=3 h=h+1; t=t+1; DetaS(2*h-1,1)=real(B2(i,1)-OrgS(2*h-1,1); %delPi DetaS(2*h,1)=real(PVU(t,1)2+imag(PVU(t,1)2-real(B2(i,3)2-imag(B2(i,3)2; %delUi endend% DetaS%创建I，用于存储节点电流参数i=zeros(n-1,1);h=0;for i=1:n if i=isb h=h+1; I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/
11、conj(B2(i,3);%conj求共轭 endend%创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对PQ节点的处理 if B2(i,6)=2 h=h+1; for j=1:n if j=isb k=k+1; if i=j %对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3
12、)+real(I(h,1); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1); else %非对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*
13、k); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1); end if k=(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行 k=0; end end end endendk=0;for i=1:n %对PV节点的处理 if B2(i,6)=3 h=h+1; for j=1:n if j=isb k=k+1; if i=j %对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j
14、)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=2*imag(B2(i,3); Jacbi(2*h,2*k)=2*real(B2(i,3); else %非对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=0; Jacbi(2*h
15、,2*k)=0; end if k=(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行 k=0; end end end endenddisp(初始雅可比矩阵为：);disp(Jacbi); %求解修正方程，获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS; %inv矩阵求逆% DetaU%修正节点电压j=0;for i=1:n %对PQ节点处理 if B2(i,6)=2 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendfor i=1:
16、n %对PV节点的处理 if B2(i,6)=3 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endend% B2%开始循环*while abs(max(DetaU)pr OrgS=zeros(2*n-2,1); h=0; j=0; for i=1:n if i=isb&B2(i,6)=2 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B
17、2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end end end for i=1:n if i=isb&B2(i,6)=3 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+re
18、al(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end end end% OrgS %创
19、建DetaS h=0; for i=1:n if i=isb&B2(i,6)=2 h=h+1; DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1); end end t=0; for i=1:n if i=isb&B2(i,6)=3 h=h+1; t=t+1;% DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h-1,1)=real(B2(i,1)-OrgS(2*h-1,1); DetaS(2*h,1)=real(PVU(t,1)2
20、+imag(PVU(t,1)2-real(B2(i,3)2-imag(B2(i,3)2; end end% DetaS %创建I i=zeros(n-1,1); h=0; for i=1:n if i=isb h=h+1; I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3); end end% I %创建Jacbi Jacbi=zeros(2*n-2); h=0; k=0; for i=1:n if B2(i,6)=2 h=h+1; for j=1:n if j=isb k=k+1; if i=j Jacbi(2*h-1,2*k-1)
21、=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1); else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y
22、(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1); end if k=(n-1) k=0; end end end end end k=0; for i=1:n if B2(i,6)=3 h=h+1; for j=1:n if j=isb k=k+1; if i=j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*r
23、eal(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=2*imag(B2(i,3); Jacbi(2*h,2*k)=2*real(B2(i,3); else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*r
24、eal(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=0; Jacbi(2*h,2*k)=0; end if k=(n-1) k=0; end end end end end% Jacbi DetaU=zeros(2*n-2,1); DetaU=inv(Jacbi)*DetaS;% DetaU %修正节点电压 j=0; for i=1:n if B2(i,6)=2 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); end end for i=1:n if B2(i,6
25、)=3 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); end end% B2 Times=Times+1; %迭代次数加1enddisp(迭代次数为：);disp(Times);disp(收敛时电压修正量为：);disp(DetaU);for k=1:n E(k)=B2(k,3); e(k)=real(E(k); f(k)=imag(E(k); V(k)=sqrt(e(k)2+f(k)2); sida(k)=atan(f(k)./e(k)*180./pi;end%= 计算各输出量 =disp(各节点的实际电压标幺值E为
26、(节点号从小到大排列)：);disp(E); %显示各节点的实际电压标幺值E用复数表示disp(-);disp(各节点的电压大小V为(节点号从小到大排列)：);disp(V); %显示各节点的电压大小V的模值disp(-);disp(各节点的电压相角sida为(节点号从小到大排列)：);disp(sida); %显示各节点的电压相角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); %计算各节点的功率 S = 电压 X 注入电流的共轭值enddisp(各
27、节点的功率S为(节点号从小到大排列)：);disp(S); %显示各节点的注入功率Sline=zeros(n1,5);disp(-);disp(各条支路的首端功率Si为(顺序同您输入B1时一致)：);for i=1:n1 p=B1(i,1); q=B1(i,2); Sline(i,1)=B1(i,1); Sline(i,2)=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)*conj(1./(B1(i,3)*B1(i,5); Siz(i)=Si(p,q); el
28、se Si(p,q)=E(p)*(conj(E(p)*(1-B1(i,5)/B1(i,3)+(conj(E(p)-conj(E(q)*(B1(i,5)/B1(i,3); Siz(i)=Si(p,q); end SSi(p,q)=Si(p,q); Sline(i,3)=Siz(i); ZF=S(,num2str(p),num2str(q),)=,num2str(SSi(p,q); disp(ZF); enddisp(-);disp(各条支路的末端功率Sj为(顺序同您输入B1时一致)：);for i=1:n1 p=B1(i,1);q=B1(i,2); if B1(i,6)=0 Sj(q,p)=E(
29、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)*(B1(i,5)*(B1(i,5)-1)/B1(i,3)+(conj(E(q)-conj(E(p)*(B1(i,5)/B1(i,3); Sjy(i)=Sj(q,p); end SSj(q,p)=Sj(q,p); Sline(i,4)=Sjy(i); ZF=S(,num2str(q),num2str(p),)=,num2str(SSj(q,p); disp(ZF);enddisp(-);disp(各条支路的功率损耗DS为(顺序同您输入B1时一致)：);for i=1:n1 p=B1(i,1); q=B1(i,2); DS(i)=Si(p,q)+Sj(q,p); DDS(i)=DS(i); Sline(i,5)=DS(i); ZF=DS(,num2str(p),num2str(q),)=,num2str(DDS(i); disp(ZF);enddisp(-);disp(各支路首端编号 末端编号 首端功率 末端功率 线路损耗);disp(Sline);专心-专注-专业