牛拉法潮流计算

自动化07-1班段佳 07051101function nl;%------------------------------------------------------------------------%===================================================================%======================牛顿——拉夫逊法==============================%===========================潮流计算=================================%===================================================================%-----------------------------------------------------------------------% % %－－－－－－－－－－－－－－－使用说明部分－－－－－－－－－－－－－－－－－－－－－－－－－－－display('% %本程序的功能是用牛顿——拉夫逊法进行潮流计算');display('% %本程序要求用户按照一定的格式将电力系统的参数制成excel表格，系统运行时将从excel中加载这些参数，随后后即可进行潮流计算');display('% %为了方便运算，用户再给系统节点进行编号时，请按照先PQ节点，再PV节点，最后平衡节点的顺序从小到大编号');display('% %电力系统潮流计算excel格式——支路参数：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳;5、支路的变比K:1；6、支路首端处于K侧为1，1侧为0'); display('% %电力系统潮流计算excel格式——节点参数：1、节点号；2、电压大小；3、相位角；4、发电机有功；5、发电机无功；6、负载有功；7、负载无功；8、节点类型'); %===================================================================%==============================数据准备==============================%===================================================================% %---------------------电力系统数据加载部分-----------------------------------------------clearx=0;Branch=0;%支路参数Note=0;%节点参数[filename, pathname] = uigetfile('*.xls', 'please choose the excel file with your powersystem parameters ');%从外部excel导入电力系统潮流计算相关参数tryif filename ~= 0x=xlsread([pathname,filename],'sheet1', 'A3:F3');Branch=xlsread([pathname,filename],'sheet1', 'A5:G10');%读支路参数Note=xlsread([pathname,filename],'sheet1', 'A15:H19');%读节点参数endcatch%进行出错处理errmsg = lasterr;errordlg(errmsg,'Save as Error');rethrow(lasterror);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、支路的变比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(:,4);%各节点发电机有功功率初始值(标幺值)Gq=Note(:,5);%各节点发电机无功功率初始值(标幺值)Lp=Note(:,6);%各节点负载有功功率初始值(标幺值)Lq=Note(:,7);%各节点负载无功功率初始值(标幺值)type=Note(:,8);%节点类型，PQ节点=1 ，PV节点=2 ，平衡节点=3for h=1:nPs(h)=Gp(h)-Lp(h);%各节点注入的有功功率Qs(h)=Gq(h)-Lq(h);%各节点注入的无功功率end% % %---------------------导纳矩阵计算部分-----------------------------------------------------Y=zeros(n);for h=1:nl %支路数if B1(h,6)==0 %左节点处于低压侧（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);elsep=B1(h,1);q=B1(h,2); %1、支路首端号；2、末端号；Y(p,q)=Y(p,q)-1./(B1(h,3)*B1(h,5));%非对角元 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,5)^2)+B1(h,4);endend%导纳矩阵显示disp('导纳矩阵 Y=');disp(Y)% % %－－－－－－－－－－－－－OK,至此潮流计算所需的数据已经准备好了－－－－－－－－－－－－－－－－%===================================================================%==============================潮流计算==============================%===================================================================%u(i)=e(i)+jf(i);Y(ij)=G(ij)+jB(ij);G=real(Y);B=imag(Y);%分解出导纳阵的实部和虚部%============================计算失配功率初始值detp\detq==========================for h=1:n-1s=0;for j=1:ns=s+U(j)*(G(h,j)*cos(a(h)-a(j))+B(h,j)*sin(a(h)-a(j)));endP(h)=U(h)*s;endfor 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;endfor h=1:n-1detp(h)=Ps(h)-P(h);endfor h=1:mdetq(h)=Qs(h)-Q(h);end%============================不满足精度要求则进入循环========================== while(max(abs(detp))>=pr|max(abs(detq))>=pr)%%不满足精度要求则循环%=================================求取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(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);endfor h=1:mU(h)=U(h)+detu(h);endk=k+1;fprintf('迭代次数k=%d\n',k);disp('节点电压大小（标幺值）');disp(U);disp('节点电压相位角（弧度）');disp(a);%===========================迭代过程中的失配功率detp\detq===========================for h=1:n-1s=0;for j=1:ns=s+U(j)*(G(h,j)*cos(a(h)-a(j))+B(h,j)*sin(a(h)-a(j)));endP(h)=U(h)*s;endfor 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;endfor h=1:n-1detp(h)=Ps(h)-P(h);endfor h=1:mdetq(h)=Qs(h)-Q(h);enddisp('迭代过程中的有功失配功率（标幺值）');disp(detp);disp('迭代过程中的无功失配功率（标幺值）');disp(detq);end% % %－－－－－－－－－－－－－OK,至此潮流计算已经完成了－－－－－－－－－－－－－－－－%===================================================================%==============================计算结果输出到工作区========================== %===================================================================%=================================迭代次数、各节点电压和视在功率==============================disp('计算结果');fprintf('总的迭代次数k=%d\n',k);disp('-----------------------------------------------------');disp('各节点电压大小（标幺值）为(节点号从小到大排列)');disp(U);disp('各节点电压相位角（角度）为(节点号从小到大排列)');A=a*180/pi;disp(A);disp('-----------------------------------------------------');disp('各节点视在功率（标幺值）为(节点号从小到大排列)');S=P+Q*i;disp(S);%=============================各条支路功率损耗和总损耗========================= ZSH=0;DS=zeros(nl,1);for h=1:nlp=B1(h,1);q=B1(h,2);DS(h)=S(p)-S(q);ZSH=ZSH+DS(h);DDS(h)=DS(h)*SB;ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(h)),' (MVA) 标么值：',num2str(DS(h))];disp(ZF);enddisp('-----------------------------------------------------');disp(['总损耗为：ZSH=',num2str(ZSH*SB),' (MVA) 标么值：',num2str(ZSH)]);%==============================结果输出到原excel========================== result0=U;%电压result1=A;%相位result2=P;%节点有功result3=Q;%节点无功result4=real(DS);%线路有功损耗result5=imag(DS);%线路无功损耗result6=real(ZSH);%系统总有功损耗result7=imag(ZSH);%系统总无功损耗[filename1, pathname1] = uiputfile('*.xls', 'put the result into the excel with your powersystem parameters ');%从外部excel导入电力系统潮流计算相关参数tryif filename1 ~= 0xlswrite([pathname1,filename1],result0 , 'sheet1', 'J3');xlswrite([pathname1,filename1],result1 , 'sheet1', 'K3');xlswrite([pathname1,filename1],result2 , 'sheet1', 'L3');xlswrite([pathname1,filename1],result3 , 'sheet1', 'M3');xlswrite([pathname1,filename1],result4 , 'sheet1', 'N3');xlswrite([pathname1,filename1],result5 , 'sheet1', 'O3');xlswrite([pathname1,filename1],result6 , 'sheet1', 'P3');xlswrite([pathname1,filename1],result7 , 'sheet1', 'Q3');endcatch%进行出错处理errmsg = lasterr;errordlg(errmsg,'Save as Error');rethrow(lasterror);end%==============================打开excel查看计算结果========================== winopen([pathname1,filename1]);% % %－－－－－－－－－－－－－OK,至此潮流计算已经全部完成－－－－－－－－－－－－－－－－% % %－－－－－－－－－－－－－O(∩_∩)O哈！－－－－－－－－－－－－－－－－% %本程序的功能是用牛顿——拉夫逊法进行潮流计算% %本程序要求用户按照一定的格式将电力系统的参数制成excel表格，系统运行时将从excel中加载这些参数，随后后即可进行潮流计算% %为了方便运算，用户再给系统节点进行编号时，请按照先PQ节点，再PV节点，最后平衡节点的顺序从小到大编号% %电力系统潮流计算excel格式——支路参数：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳;5、支路的变比K:1；6、支路首端处于K侧为1，1侧为0% %电力系统潮流计算excel格式——节点参数：1、节点号；2、电压大小；3、相位角；4、发电机有功；5、发电机无功；6、负载有功；7、负载无功；8、节点类型导纳矩阵 Y=6.3110 -20.4022i -3.5587 +11.3879i -2.7523 + 9.1743i 0 0-3.5587 +11.3879i 8.5587 -31.0093i -5.0000 +15.0000i 0 + 4.9889i-2.7523 + 9.1743i -5.0000 +15.0000i 7.7523 -28.7757i 0 0 + 4.9889i0 0 + 4.9889i 0 0 - 5.2493i0 0 0 + 4.9889i 0 0 - 5.2493iJacobi =-20.5622 11.3879 9.1743 0 -6.3110 3.5587 2.752311.3879 -31.3768 15.0000 4.9889 3.5587 -8.5587 5.00009.1743 15.0000 -29.1632 0 2.7523 5.0000 -7.75230 4.9889 0 -4.9889 0 0 06.3110 -3.5587 -2.7523 0 -20.2422 11.3879 9.1743-3.5587 8.5587 -5.0000 0 11.3879 -30.6418 15.0000-2.7523 -5.0000 7.7523 0 9.1743 15.0000 -28.3882迭代次数k=1节点电压大小（标幺值）1.00361.02971.03271.00001.0000节点电压相位角（弧度）-0.0900-0.0577-0.06120.0425迭代过程中的有功失配功率（标幺值）0.0193-0.0007-0.0140迭代过程中的无功失配功率（标幺值）-0.0148-0.0387-0.0270Jacobi =-21.0649 11.6431 9.4218 0 -5.5312 4.0561 3.1247 11.8809 -32.9618 15.9694 5.1115 3.2952 -9.0811 5.2601 9.5859 15.9316 -30.6597 0 2.5776 5.3736 -8.2678 0 5.1115 0 -5.1115 0 -0.5140 0 7.1808 -4.0561 -3.1247 0 -20.0304 11.6431 9.4218 -3.2952 9.0693 -5.2601 -0.5140 11.8809 -32.7992 15.9694 -2.5776 -5.3736 8.2664 0 9.5859 15.9316 -30.7138迭代次数k=2节点电压大小（标幺值）0.99371.02071.02401.00001.0000节点电压相位角（弧度）-0.0904-0.0587-0.06200.0397迭代过程中的有功失配功率（标幺值）0.0019-0.0010-0.0008-0.0002迭代过程中的无功失配功率（标幺值）0.0046-0.0041Jacobi =-20.6812 11.4294 9.2518 0 -5.4239 3.9739 3.0648 11.6585 -32.4210 15.6951 5.0675 3.2410 -8.9173 5.1741 9.4110 15.6605 -30.1704 0 2.5340 5.2777 -8.1299 0 5.0675 0 -5.0675 0 -0.5002 0 7.0387 -3.9739 -3.0648 0 -19.6079 11.4294 9.2518 -3.2410 8.9154 -5.1741 -0.5002 11.6585 -32.1891 15.6951 -2.5340 -5.2777 8.1284 0 9.4110 15.6605 -30.1786迭代次数k=3节点电压大小（标幺值）0.99351.02031.02371.00001.0000节点电压相位角（弧度）-0.0905-0.0587-0.06200.0397迭代过程中的有功失配功率（标幺值）1.0e-004 *0.6037-0.2358-0.3406-0.0364迭代过程中的无功失配功率（标幺值）1.0e-003 *0.1808-0.1108-0.1555-20.6709 11.4238 9.2471 0 -5.4240 3.9719 3.0632 11.6527 -32.4023 15.6839 5.0657 3.2395 -8.9101 5.1705 9.4062 15.6495 -30.1528 0 2.5328 5.2739 -8.1234 0 5.0657 0 -5.0657 0 -0.5000 0 7.0351 -3.9719 -3.0632 0 -19.6066 11.4238 9.2471 -3.2395 8.9101 -5.1705 -0.5000 11.6527 -32.1625 15.6839 -2.5328 -5.2739 8.1233 0 9.4062 15.6495 -30.1531迭代次数k=4节点电压大小（标幺值）0.99351.02031.02361.00001.0000节点电压相位角（弧度）-0.0905-0.0587-0.06200.0397迭代过程中的有功失配功率（标幺值）1.0e-005 *0.1199-0.0264-0.0861-0.0075迭代过程中的无功失配功率（标幺值）1.0e-005 *0.3521-0.1563-0.3785计算结果总的迭代次数k=4-----------------------------------------------------各节点电压大小（标幺值）为(节点号从小到大排列)1.02031.02361.00001.0000各节点电压相位角（角度）为(节点号从小到大排列)-5.1825-3.3648-3.55382.2723-----------------------------------------------------各节点视在功率（标幺值）为(节点号从小到大排列)-0.8055 - 0.5320i0.0000 - 0.1200i0.0000 + 0.0000i0.5000 + 0.1837iDS(1,2)=-80.5501-41.2005i (MVA) 标么值：-0.8055-0.41201iDS(1,3)=-80.5502-53.2007i (MVA) 标么值：-0.8055-0.53201iDS(2,3)=-5.97263e-005-12.0002i (MVA) 标么值：-5.9726e-007-0.12iDS(4,2)=50+30.3696i (MVA) 标么值：0.5+0.3037iDS(5,3)=-8.6117e-005-0.00037855i (MVA) 标么值：-8.6117e-007-3.7855e-006i -----------------------------------------------------总损耗为：ZSH=-111.1005-76.03228i (MVA) 标么值：-1.111-0.76032i。