《牛拉法潮流计算》由会员分享，可在线阅读，更多相关《牛拉法潮流计算（12页珍藏版）》请在金锄头文库上搜索。
1、自动化 07-1 班 段佳 07051101function nl; %- %= %=牛顿拉夫逊法= %=潮流计算= %= %- % % %使用说明部分 display(% %本程序的功能是用牛顿拉夫逊法进行潮流计算); display(% %本程序要求用户按照一定的格式将电力系统的参数制成 excel 表格，系统运行 时将从 excel 中加载这些参数，随后后即可进行潮流计算); display(% %为了方便运算，用户再给系统节点进行编号时，请按照先 PQ 节点，再 PV 节点， 最后平衡节点的顺序从小到大编号); display(% %电力系统潮流计算 excel 格式支路参数 ：1、支
2、路首端号；2、末端号； 3、支路阻抗；4、支路对地电纳;5、支路的变比 K:1；6、支路首端处于 K 侧为 1，1 侧为 0);display(% %电力系统潮流计算 excel 格式节点参数 ：1、节点号；2、电压大小；3、 相位角；4、发电机有功；5、发电机无功；6、负载有功；7、负载无功；8、节点类型);%= %=数据准备= %= % %-电力系统数据加载部分- clear x=0; Branch=0;%支路参数 Note=0;%节点参数filename, pathname = uigetfile(*.xls, please choose the excel file with your
3、 powersystem parameters );%从外部 excel 导入电力系统潮流计算相关参数tryif filename = 0 x=xlsread(pathname,filename,sheet1, A3:F3);Branch=xlsread(pathname,filename,sheet1, A5:G10);%读支路参数Note=xlsread(pathname,filename,sheet1, A15:H19);%读节点参数end catch%进行出错处理errmsg = lasterr;errordlg(errmsg,Save as Error);rethrow(laster
4、ror);end % %-支路参数初始化部分- SB=100;UB=220;n=1;m=1;pr=0.0001; SB=x(5);%功率基准值 UB=x(6);%电压基准值 n=x(1);%节点数 nl=x(2);%支路数 m=x(3);%PQ 节点的个数 pr=x(4);%误差精度 B1(:,1)=Branch(:,1);%1、支路首端号 B1(:,2)=Branch(:,2);%2、末端号 B1(:,3)=Branch(:,3)+Branch(:,4)*i;%3、支路阻抗 B1(:,4)=Branch(:,5)*i;%4、支路对地电纳 B1(:,5)=Branch(:,6);%5、支路的变
5、比 K:1； B1(:,6)=Branch(:,7);%6、支路首端处于 K 侧为 1，1 侧为 0% % % %-节点参数初始化部分- U=ones(n,1); a=zeros(n,1); Ps=zeros(n,1); Qs=zeros(n,1); P=zeros(n,1); Q=zeros(n,1); detp=zeros(n-1,1); detq=zeros(m,1); deta=zeros(n-1,1); detu=zeros(m,1); k=0;%迭代次数 U=Note(:,2);%各节点电压初始值(标幺值) a=Note(:,3);%各节点电压相位初始值（弧度） Gp=Note(:
6、,4);%各节点发电机有功功率初始值(标幺值) Gq=Note(:,5);%各节点发电机无功功率初始值(标幺值) Lp=Note(:,6);%各节点负载有功功率初始值(标幺值) Lq=Note(:,7);%各节点负载无功功率初始值(标幺值) type=Note(:,8);%节点类型，PQ 节点=1 ，PV 节点=2 ，平衡节点=3for h=1:n Ps(h)=Gp(h)-Lp(h);%各节点注入的有功功率 Qs(h)=Gq(h)-Lq(h);%各节点注入的无功功率end % % %-导纳矩阵计算部分- Y=zeros(n); for h=1:nl%支路数 if B1(h,6)=0%左节点处于
7、低压侧 （6、支路首端处于 K 侧为 1，1 侧为 0） p=B1(h,1);q=B1(h,2); %1、支路首端号；2、末端号；Y(p,q)=Y(p,q)-1./B1(h,3); %非对角元 3、支路阻抗；4、支路对地电纳； 5、支路的变比；Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1./B1(h,3)+B1(h,4);Y(q,q)=Y(q,q)+1./B1(h,3)+B1(h,4);else p=B1(h,1);q=B1(h,2); %1、支路首端号；2、末端号；Y(p,q)=Y(p,q)-1./(B1(h,3)*B1(h,5);%非对角元 3、支路阻抗；4、支路对地电 纳
8、；5、支路的变比；Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1./B1(h,3)+B1(h,4);Y(q,q)=Y(q,q)+1./(B1(h,3)*B1(h,5)2)+B1(h,4);end end %导纳矩阵显示 disp(导纳矩阵 Y=); disp(Y) % % %OK,至此潮流计算所需的数据已经准备好了 %= %=潮流计算= %= %u(i)=e(i)+jf(i);Y(ij)=G(ij)+jB(ij); G=real(Y);B=imag(Y);%分解出导纳阵的实部和虚部 %=计算失配功率初始值detpdetq= for h=1:n-1s=0;for j=1:ns=s+
9、U(j)*(G(h,j)*cos(a(h)-a(j)+B(h,j)*sin(a(h)-a(j);endP(h)=U(h)*s; end for h=1:n-1s=0;for j=1:ns=s+U(j)*(G(h,j)*sin(a(h)-a(j)-B(h,j)*cos(a(h)-a(j);endQ(h)=U(h)*s; end for h=1:n-1detp(h)=Ps(h)-P(h); end for h=1:mdetq(h)=Qs(h)-Q(h); end %=不满足精度要求则进入循环= while(max(abs(detp)=pr|max(abs(detq)=pr)%不满足精度要求则循环%
10、=求取 Jacobi 矩阵=H=zeros(n-1,n-1);N=zeros(n-1,m);K=zeros(m,n-1);L=zeros(m,m);for h=1:n-1for j=1:n-1if h=jH(h,j)=U(h)2*B(h,j)+Q(h);elseH(h,j)=-U(h)*U(j)*(G(h,j)*sin(a(h)-a(j)-B(h,j)*cos(a(h)- a(j);endendendfor h=1:n-1for j=1:mif h=jN(h,j)=-U(h)2*G(h,j)-P(h);elseN(h,j)=-U(h)*U(j)*(G(h,j)*cos(a(h)-a(j)+B(
11、h,j)*sin(a(h)-a(j);endendendfor h=1:mfor j=1:n-1if h=jK(h,j)=U(h)2*G(h,j)-P(h);elseK(h,j)=U(h)*U(j)*(G(h,j)*cos(a(h)-a(j)+B(h,j)*sin(a(h)-a(j);endendendfor h=1:mfor j=1:mif h=jL(h,j)=U(h)2*B(h,j)-Q(h);elseL(h,j)=-U(h)*U(j)*(G(h,j)*sin(a(h)-a(j)-B(h,j)*cos(a(h)- a(j);endendend %=解修正方程，得到修正量detu,deta=Jacobi=H N;K L;display(Jacobi);dets=detp;detq;solutions=-inv(Jacobi)*dets;deta=solutions(1:n-1,:);detu=solutions(n:n-1+m,:);%=迭代过程中的电压=for h=1:n-1a(h)=a(h)+deta(h)
